#=====================================================================
#    Homework for 640:640 Experimental Math
#      Due 1/28/2010
#      David Wilson
#=====================================================================
#   Challenge to improve PrimeList(N)
#=====================================================================


#Challenge( 5 dollars for the best entry): write a faster program for
#PrimeList(N) ; (Note: Brian Nakamura told me that the Sieve of 
#Eratosthenes (SA(N) in our program) is pretty good for the first million
#primes, so try to use it!)

Help:=proc(): print(`SA(n)`):
print(`IsDiv(n,L), OneStep(L), PrimeList(N), PrimeListM(N)`):
print(`PrimeListNew(n), TimeDiff(m,n)`);
end :

#============================================================================
#       New Code
#============================================================================

#PrimeListNew(n): The list of the first n primes
PrimeListNew:=proc(n) local A,N:

A:=SA(n):
for N from nops(A) to n-1 do
A:=OneStep(A):
od:
A:
end:


#SA(n): The list of primes <= n
SA:=proc(n) local i,j:
{seq(i,i=1..n)} minus 
{seq(seq(j*i,i=2..n/j),j=2..trunc(sqrt(n)))} minus {1}:

end:

#TimeDiff(m,n)
TimeDiff:=proc(m,n) local i,s,t,A:

A:=Array(1..n,1..2):

for i from 1 to n do
s:=time():
PrimeList(m^i):
A[i,1]:=time()-s:

t:=time():
PrimeListNew(m^i):
A[i,2]:=time()-t:

od:

A:
end:


#============================================================================
#       Old Code
#============================================================================



#IsDiv(n,L): inputs a pos. integer n and a list
#of pos. integers L, and outputs true if
#n is divisible by any of the members of L

IsDiv:=proc(n,L) local i:
for i from 1 to nops(L) do 
 if n mod L[i]=0 then
  RETURN(true):
 fi:
od:
false:

end:


#OneStep(L): inputs the list of the first primes
#and outputs the same list with one more prime that
#comes right after
OneStep:=proc(L) local i,n:
#n is the largest prime so-far
n:=L[nops(L)]:

 for i from n+1 while IsDiv(i,L) do od:

[op(L),i]:

end:

#The list of the first N prime numbers, our way
PrimeList:=proc(N) local i,L:
L:=[2]:
for i from 2 to N do
 L:=OneStep(L):
od:
L:
end:

#The list of the first N prime numbers, the Maple way
PrimeListM:=proc(N) local i:
 [seq(ithprime(i),i=1..N)]:
end: