Help:=proc(): print(`BonnSieve(n,k), `):
print(` ApplySieve(Se,P,Si), MakeSeP(n,k)`): 
print(` BrunSieve(k,N), FullBrunSieve(k,N) `):
print(` HMPB(x,y) , BrunSC(r,n,a)`):
end:
with(combinat):
with(numtheory):

#{{set}}: A collection of subsests of {1,2, ....n}
#such that Sum((-1)^s*ElementsWithProperties[s], s in the 
#{set} gives a lower or upper bounds (if k is even or odd)
#(either resp. or not)
BonnSieve:=proc(n,k) local i,N:
N:={seq(i,i=1..n)}:
{seq(op(choose(N,i)),i=0..k)}:

end:


#ApplySieve(Se,P,Si): Inputs a universal se Se
#and a list of "properties" described explicitly
#by a list of subsets of Se P=[P1,P2, ...]
#(meaning Pi is the set of objectts having property i
#(and possiblly other properties), and it also inputs
# a sieve Si. Outputs the "estimate" implied by Si
#and the true value of the set of objects that has
#none of the properties
#For example
#ApplySieve({Emilie,Aisha,Daniela},
#[{Emilie,Aisha},{Emilie,Daniela}], {{},{1},{2}})
#should return
#TrueValue=0
#Estimate=3-2-2=-1
ApplySieve:=proc(Se,P,Si)  local Es, TV,josh,kellen,Si1:

if not {} in Si then
 ERROR(`Bad input`):
fi:

TV:=nops(Se minus {seq(op(P[i]),i=1..nops(P))}):

Si1:=Si minus {{}}:


Es:=add((-1)^nops(josh)*nops(`intersect`(seq(P[kellen],
kellen in josh))),josh in Si1)  +  nops(Se):
Es,TV:

end:





#MakeSeP(n,k): inputs pos. integers n and k
#and  outputs the set {1,...,n} followed by
#a list of k random subsets of {1,...n}
#each with a random number of elements
MakeSeP:=proc(n,k) local Se,P,ra,i,car,j:
Se:={seq(i,i=1..n)}:
P:=[]:
ra:=rand(1..n):

for i from 1 to k do
car:=ra():
P:=[op(P),{seq(ra(),j=1..car)}]:
od:

Se,P:

end:


#k=Number of ordered properties
#BrunSieve1(k,N) inputs the number of properties k
#and a list N of weakly decrreasing pos. integers
#N=[N1,N2,N3, Nr] with r<=k
#outputs all subsets of {1, ..., k} written
#in dec. order i1>i2>i3>.. such that
#i1<=N1 and i2<=N2, ...ir<=Nr  
BrunSieve1:=proc(k,N) local i,r,N1,emilie,s,S,T:
option remember:
r:=nops(N):
if r=1 then
 RETURN({seq([i],i=1..min(N[1],k))}):
fi:

N1:=[op(2..r,N)]:

S:={seq(i,i=1..min(N[1],k))}:

T:={}:

for s in S do
 emilie:=BrunSieve1(s-1,N1):
 T:=T union {seq([s,op(e)],e in emilie)}:
od:
T:

end:

#k=Number of ordered properties
#BrunSieve(k,N) inputs the number of properties k
#and a list N of weakly decrreasing pos. integers
#N=[N1,N2,N3, Nr] with r<=k
#outputs all subsets of exactly r 
#elements of {1, ..., k} written as
#sets
#i1<=N1 and i2<=N2, ...ir<=Nr  
BrunSieve:=proc(k,N) local brian,b:
brian:=BrunSieve1(k,N):
{seq(convert(b,set), b in brian)} union {{}}:
end:

#k=Number of ordered properties
#FullBrunSieve(k,N) inputs the number of properties k
#and a list N of weakly decrreasing pos. integers
#N=[N1,N2,N3, Nr] with r<=k of cardinality <=r
#outputs all subsets of {1, ..., k} written as
#sets
#i1<=N1 and i2<=N2, ...ir<=Nr  
FullBrunSieve:=proc(k,N) local r,i:
r:=nops(N):

{seq(op(BrunSieve(  k  ,[op(1..i,N)])),i=1..r)}:

end:


#How many prime numbers between x and y
HMPBslow:=proc(x,y) local i,co:

co:=0:

for i from x to y-1 do
  if isprime(i) then
     co:=co+1:
  fi:
od:

co:

end:

HMPB:=proc(x,y): pi(y-1)-pi(x-1):end:



#BrunSC(r,n,a): inputs a pos. integer r and a real
#number a>1 and outputs the Brun Staircase
#before Eq. (15) in Brun's paper as form of lists
#N=[sigma1+....+sigman, ...., sigman]
BrunSC:=proc(r,n,a) local pr,sigma,i,N:
pr:=ithprime(r):

for i from 1 to n do
#sigma[i]: number of prime numbers between
#pr^(1/a^(i-1)) and pr^(1/a^i)
sigma[i]:= HMPB( trunc(pr^(1/a^(i))),
trunc(pr^(1/a^(i-1)))):
od:

[seq(sigma[i],i=1..n)];


#Nb:=[sigma[n],sigma[n],
#sigma[n-1]+sigma[n],
#sigma[n-1]+sigma[n],
#sigma[n-2]+sigma[n-1]+sigma[n],


N:=[seq(add(sigma[n-j],j=0..i-1), i=1..n)] :

N:=[seq(N[n-i],i=0..n-1)]:


[N[1],seq(op([N[i],N[i]]),i=2..n)]:


end: