#OK to post homework
#Joseph Koutsoutis, 03-31-2024, Assignment 19

read `C19.txt`:

#1 done

#2 done

#3

P := proc(n) local m, j:
 1 + add(trunc( trunc( n / (1 +  add(trunc((1 + (j-1)! )/j - trunc( (j-1)!/j ) ), j=2..m)  ) )^(1/n) ), m=1..2^n):
end:

# I checked that the first 11 outputs of P matched the first 11 primes.

t := proc(n) local p:
 2 + n*trunc( 1 / ( 1 + add( trunc( (n+2)/p - trunc( (n+1)/p ) ) , p=2..n+1) ) ):
end:

# {seq(isprime(t(i)), i=1..1000)} outputted {true}.

Pi1 := proc(n) local k:
 trunc( add( (sin( Pi * GAMMA(k) / (2*k)  ))^2 , k=1..n) ):
end:

# I tested Pi1(i) for i from 2 to 30 and it correctly counted the number of primes between 1 and i.

#4

# TestSameKey(d,K) confirms that Alice and Bob end up with the same key when using random primes with
# d digits. We test this K times.
TestSameKey := proc(d, K) local i,g,P,Alice,Bob,gab,gba:
 for i from 1 to K do:
  P := FindSGP(d,100):
  while P = FAIL do P := FindSGP(d,100): od:
  g := FindPrimiSGP(P):
  Alice := AliceToBob(P,g):
  Bob := BobToAlice(P,g):
  gab := Alice[2]&^Bob[1] mod P:
  gba := Bob[2]&^Alice[1] mod P:
  if gab <> gba then return(FAIL): fi:
 od: 

 return(SUCCESS):
end:

# I ran TestSameKey(10,100) and TestSameKey(100,100) and both returned SUCCESS


#5

DL := proc(x,P,g) local a,g1:
 g1 := 1:
 for a from 1 to P-1 do:
  g1 := g1 * g mod P:
  if g1 = x then return(a): fi:
 od:
 return(FAIL):
end: 


# When I used an 8 digit prime, it took my computer 3 minutes to solve DL(1,P,g), so I did not try
# P with more digits.


#6


# I followed the pseudocode from the wikipedia article for this implementation
# (https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test)
MillerRabin := proc(n,K) local s,d,a,x,y,i,j,ra:
 if n = 2 or n=3 then return(true): fi:
 if n mod 2 = 0 then return(false): fi:

 s := 0:
 d := n-1:
 while d mod 2 = 0 do:
  s += 1:
  d := d / 2:
 od:

 ra := rand(2..n-2):

 for i from 1 to K do:
  a := ra():
  x := a&^d mod n:
  for j from 1 to s do:  
   y := x^2 mod n:
   if y = 1 and x <> 1 and x <> n-1 then return(false): fi:
   x := y:
  od:
  if y <> 1 then return(false): fi:
 od:

 return(true):
end: