Help:=proc(): print(`pn(n), IsPrim(g,p), FindPrim(p), FindSGP(d,K), FindPrimSGP(P), AliceToBob(P,g), BobToAlice(P,g), FDL(x,g,P)  `):end:

pn:=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:



#IsPrimS(g,p): Is g primitive mod p?
IsPrimS:=proc(g,p) local i:
for i from 1 to p-2 do
 if g^i mod p=1 then
   RETURN(false):
 fi:
od:
true:
end:

FindPrimS:=proc(p) local i:

for i from 2 to p-1 do
 if IsPrimS(i,p) then
   RETURN(i):
 fi:
od:

FAIL:

end:

#IsPrim(g,p): Is g primitive mod p?
IsPrim:=proc(g,p) local i:
for i from 1 to p-2 do
 if g&^i mod p=1 then
   RETURN(false):
 fi:
od:
true:
end:

FindPrim:=proc(p) local i:

for i from 2 to p-1 do
 if IsPrim(i,p) then
   RETURN(i):
 fi:
od:

FAIL:

end:



#IsPrimF(g,p): Is g primitive mod p?
IsPrimF:=proc(g,p) local i,g1:
g1:=g:
for i from 1 to p-2 do
 if g1=1 then
   RETURN(false):
  else
  g1:=g1*g mod p
 fi:
od:
true:
end:
 




FindPrimF:=proc(p) local i:

for i from 2 to p-1 do
 if IsPrimF(i,p) then
   RETURN(i):
 fi:
od:

FAIL:

end:

#FindSGP(d,K): finds a random Sophie Germain prime with d digits using K tries
FindSGP:=proc(d,K) local Q,i,ra:
ra:=rand(10^d..10^(d+1)):

for i from 1 to K do
 Q:=nextprime(ra()):
 if isprime(2*Q+1) then
    RETURN(2*Q+1):
 fi:
od:
FAIL:
end:

#FindPrimSGP(P): Finds a primitive element if P is a Sophie Germain prime
FindPrimSGP:=proc(P) local Q,a:

if not isprime(P) then
   RETURN(FAIL):
fi:
Q:=(P-1)/2:
 if not isprime(Q) then
    RETURN(FAIL):
  fi:

a:=rand(2..P-1)():

if a&^Q mod P=P-1 then
  RETURN(a):
else
 a:=P-a :
 if a&^Q mod P<>P-1 then
   print(`Something terrible happened`):
   RETURN(FAIL):
  else
   RETURN(a):
 fi:
fi:
end:

AliceToBob:=proc(P,g) local a,x:
a:=rand(2..P-2)():
x:=g&^a mod P:
print(`Bob, my number to you is`,x, `of course I won't tell you the expoent (discrete logarithmic base`, g, `mod`, P):
RETURN(x):
end:


BobToAlice:=proc(P,g) local b,x:
b:=rand(2..P-2)():
x:=g&^b mod P:
print(`Bob, my number to you is`,x, `of course I won't tell you the expoent (discrete logarithmic base`, g, `mod`, P):
RETURN(x):
end:




FDL:=proc(x,g,P) local a:
for a from 1 to P-1 do
 if g&^a mod P=x then
    RETURN(a):
 fi:
od:
FAIL:
end: