#OK to post homework
#Lucy Martinez, 03-29-2024, Assignment 19


# Problem 3: Write Maple procedures pn, t(n), &Pi(n), 
# implementing these sequences in Delahaye's article. 
# Verify that they give what they are supposed to give
p:=proc(n) local m,j:
1+add(trunc(trunc(n/(1+add(trunc((1+factorial(j-1))/(j)-trunc((factorial(j-1))/j)),j=2..m)))^(1/n)),m=1..2^n):
end:

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

#NOTE: I opted to use Pii as the name of the procedure
#because Pi is already used in maple
Pii:=proc(n) local k:
trunc(evalf(add((sin((Pi*GAMMA(k))/(2*k)))^2,k=1..n))):
end:



# Problem 4: By using AliceToBob(P,g) and BobToAlice(P,g)
# verify that they share the same key (do it several times)

# ANSWER:
# Example 1: 
# P:=FindSGP(100,100); returns
# P:=13250091537996868200834998692213101306045759725382970424773733953253615521821389079749404909205668319
# g:=FindPrimiSGP(P); returns
# g:=9043800052123305937356306990256411776523304274133508461011332637380011700925384864697123139902474788
# Now, Alice will give Bob ga[2] since ga[1]=a which she keeps
# ga:=AliceToBob(P,g); returns
# ga:=[6541881294504954572786947171376621609818466061896192442184091280542543710515721947261836974799079152, 2031440483608008988346654896243144433629467354213204465477287789589559418867518097003256321473527943]
# Now, Bob will give Alice gb[2] since gb[1]=b which he keeps
# gb:=BobToAlice(P,g); returns
# gb:=[7522441590978620880893419751386251149443436428411974078159249215796575437546791860017740801107395561, 4705914402937975048550299940176604553821061344021116083912864245868236685538815963222249722854477806]
# Now, both Alice and Bob compute the key, which should be the same, as follows:
# gab:=gb[2]&^ga[1] mod P; returns (Alice's computation)
# gab:=1571805319623314060249738059903814667385385616153436691353330470791336863327003148695183362849716386
# gba:=ga[2]&^gb[1] mod P; returns (Bob's computation)
# gba:=1571805319623314060249738059903814667385385616153436691353330470791336863327003148695183362849716386
# We get gab=gba as desired

# Example 2:
# P:=FindSGP(100,100); returns
# P:=13993046450972565375929551708921266219588689830176473802936473193050637988976920397363096398694483687
# g:=FindPrimiSGP(P); returns
# g:=1711605205788439359072986176744405570664898819075050832633732455945492038230632184221539368471095209
# Now, Alice will give Bob ga[2] since ga[1]=a which she keeps
# ga:=AliceToBob(P,g); returns 
# ga:=[7794528659439181996092828013571416778294961632679448806307202940700031859451794972001177741091149552, 7933343622929937709536702308705810417852277316996330861825142122930811305488306406500102214438764955]
# Now, Bob will give Alice gb[2] since gb[1]=b which he keeps
# gb:=BobToAlice(P,g); returns 
# gb:=[6333885204701113560065113308548343109558375136526428225407153913387140983681584276255299129470234867, 1789766198158123076287929574678683768189418204009155349122717003410249005922678225349121001881655752]
# Now, both Alice and Bob compute the key, which should be the same, as follows:
# gab:=gb[2]&^ga[1] mod P; returns (Alice's computation)
# gab:=9554652247776631127407151032276101743702917773060341022382309111060530256978077687888467754077157696
# gba:=ga[2]&^gb[1] mod P; returns (Bob's computation)
# gba:=9554652247776631127407151032276101743702917773060341022382309111060530256978077687888467754077157696
# We get gab=gba as desired


# Example 3:
# P:=FindSGP(100,100); returns
# P:=16703792627864825516179569640251736122004592528443010329856009934897620727789536383679397645485967587
# g:=FindPrimiSGP(P); returns
# g:=6551867880735065519114735912219346700449062697443244670823082887315915950661571805985330482366350607
# Now, Alice will give Bob ga[2] since ga[1]=a which she keeps
# ga:=AliceToBob(P,g); returns 
# ga:=[11897750380378764325361090652889768837234187258779465065621374972234768406651615574640208285930409023, 4568787213786607245807927254816041423120356070745402116484399743614152662327367306485032324148417803]
# Now, Bob will give Alice gb[2] since gb[1]=b which he keeps
# gb:=BobToAlice(P,g); returns 
# gb:=[3078260254309357437241595120566333347437312535482735114495176519979382165105862368386633611488421178, 10932821497378206366865213594916466414488288635149986995461650726816492037682871054556312999176073199]
# Now, both Alice and Bob compute the key, which should be the same, as follows:
# gab:=gb[2]&^ga[1] mod P; returns (Alice's computation)
# gab:=16456399113877713256410748697007810282575665722107586842726496847673970070218067294310418193193285380
# gba:=ga[2]&^gb[1] mod P; returns (Bob's computation)
# gba:=16456399113877713256410748697007810282575665722107586842726496847673970070218067294310418193193285380
# We get gab=gba as desired



# Example 4:
# P:=FindSGP(100,100); returns
# P:=18397249603150204757623885856611765866207997948749406520061727308198074452235710600972849233370434227
# g:=FindPrimiSGP(P); returns
# g:=14877567241021167545124669647260239681142218442075671768858476437670738612586647597928801755378721150
# Now, Alice will give Bob ga[2] since ga[1]=a which she keeps
# ga:=AliceToBob(P,g); returns 
# ga:=[16587863401495229481358187989162337011872375598275701804900085066433492102665039820328991012668722110, 6951968791448420648923096326145562310966219820756545056814246006435608734942010483263186967992779087]
# Now, Bob will give Alice gb[2] since gb[1]=b which he keeps
# gb:=BobToAlice(P,g); returns 
# gb:=[8617587520528738605293690327768082835256253583613960139005237597744159863667476566562564936144284951, 14239193327092053662881136795190967888948781416649453840023891697007787529039090403907077861909730604]
# Now, both Alice and Bob compute the key, which should be the same, as follows:
# gab:=gb[2]&^ga[1] mod P; returns (Alice's computation)
# gab:=10393617316297711250144220684147107279594876019022049695845054180778260631906693312024287767786374904
# gba:=ga[2]&^gb[1] mod P; returns (Bob's computation)
# gba:=10393617316297711250144220684147107279594876019022049695845054180778260631906693312024287767786374904
# We get gab=gba as desired






# Problem 5: Write a procedure DL(x,P,g) that inputs an integer x between 1 and P-1
# and outputs an integer 'a' between 1 and P-1 such that ga=x
# (or returns FAIL if none exits).
# Convince yourself that for large P it is very slow
DL:=proc(x,P,g) local a:
if not (1<x and x<P-1) then
 return(FAIL):
fi:

for a from 2 to P-2 do
 if g*a=x then
  return(a):
 fi:
od:
FAIL:
end:



####################################### C19.txt
Help19:=proc():
 print(`FindRD(d),IsPrim(g,P),FindPrimi(P),FindSGP(d,K),FindPrimiSGP(P)`):
 print(`AliceToBob(P,g),BobToAlice(P,g)`):
end:


#FindRD(d): finds a random prime with d digits
FindRP:=proc(d) local ra:
nextprime(rand(10^(d-1)..10^(d)-1)()):
end:


#IsPrim(g,P): Is g primitive mod P? i.e. are all {g,g^2,...,g^(p-1)} different if
# g^i mod p = 1 doesn't happen until i=p-1
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:


#FindPrimi(P): Given a prime P, finds a primitive g by brute force
FindPrimi:=proc(P) local g:
for g from 2 to P-2 do
 if IsPrim(g,P) then
  return(g):
 fi:
od:
FAIL:
end:


#Finds Sophie Germain Primes
#FindSGP(d,K): finds a random Sophie Germain prime such that (Q-1)/2 has d digits
#trying K times or return FAIL
FindSGP:=proc(d,K) local Q,i,ra:
ra:=rand(10^(d-1)..10^d-1):
for i from 1 to K do
 #this is the candidate prime Q
 Q:=nextprime(ra()):
 if isprime(2*Q+1) then 
  return(2*Q+1):
 fi:
od:
FAIL:
end:


#FindPrimiSGP(P): finds a random primitive element mod P
# if P is a Sophie Germain prime
FindPrimiSGP:=proc(P) local Q,a:
Q:=(P-1)/2:
if not (isprime(P) and isprime(Q)) then
 return(FAIL):
fi:

#a is candidate
a:=rand(2..P-2)():

if a&^Q mod P =P-1 then
 return(a):
fi:
end:


#AliceToBob(P,g): alice picks her secret key 'a' from 2 to P-2
#and does not tell anyone (not even Bob) and sends Bob: g^a mod P
AliceToBob:=proc(P,g) local a:
a:=rand(3..P-3)():
[a,g&^a mod P]:
end:


#BobToAlice(P,g): alice picks her secret key 'b' from 2 to P-2
#and does not tell anyone (not even Alice) and sends Alice: g^b mod P
BobToAlice:=proc(P,g) local b:
b:=rand(3..P-3)():
[b,g&^b mod P]:
end:


####################################### C18.txt
Help18:=proc():
 print(`RW(k), BinToIn(L),InToBin(n,k),BinFun(F,L),Feistel(LR,F),revFeistel(LR,F)`):
end:


#q=2
#RW(k): a random binary word of length k
RW:=proc(k) local ra,i:
ra:=rand(0..1):
[seq(ra(),i=1..k)]:
end:


#BinToIn(L): inputs a binary word and outputs the integer whose binary
# representation it is
#convert from binary to base 10
BinToIn:=proc(L) local k,i:
k:=nops(L):
#ex: [1,0,1]=1*2^2+0*2^1+1*2^0
add(L[i]*2^(k-i),i=1..k):
end:


#InToBin(n,k): inputs a non-neg integer n<2^k into its binary representation with
#0's in front if necessary so that it is of length k
InToBin:=proc(n,k) local i:
#option remember:
if (n>=2^k or n<0) or not (type(n,integer)) then
 return(FAIL):
fi:

if k=1 then
 if n=0 then
  return([0]):
 else
  return([1]):
 fi:
fi:

if n<2^(k-1) then
 return([0,op(InToBin(n,k-1))]):
else
#here we need to remove 2^(k-1)
 return([1,op(InToBin(n-2^(k-1),k-1))]):
fi:
end:


#BinFun(F,L): converts a function on the integers (mod 2^k) to a function
# on binary words
BinFun:=proc(F,L) local k:
k:=nops(L):
InToBin(F(BinToIn(L)) mod 2^k,k):
end:

#Feistel(LR,F): The Feistel transform that takes [L,R]->[R+F(L) mod 2, L]
#For example Feistel([1,1,0,1],n->n^5+n);
Feistel:=proc(LR,F) local k,L,R:
k:=nops(LR):
if k mod 2<>0 then
 return(FAIL):
fi:
L:=[op(1..k/2,LR)]:
R:=[op(k/2+1..k,LR)]:

[op(R + BinFun(F,L) mod 2) , op(L) ]:
end:


#revFeistel(LR,F): The reverse Feistel transform that takes [L,R]->[R,L+F(R) mod 2]
#For example revFeistel([0,0,1,0],n->n^5+n);
revFeistel:=proc(LR,F) local k,L,R:
k:=nops(LR):
if k mod 2<>0 then
 return(FAIL):
fi:
L:=[op(1..k/2,LR)]:
R:=[op(k/2+1..k,LR)]:

[op(R),op(L + BinFun(F,R) mod 2)]:
end:

####################################### C17.txt
Help17:=proc():
 print(`Pols(x,d),PolsE(x,d),IsIr(P),IsPr(P,x),WisW(d,x)`):
end:


#Pols(x,d): The set of all polynomials of degree <=d in x mod 2
# that have 1 (constant term) in it
Pols:=proc(x,d) local S,s:
option remember:
if d=0 then
 return({1}):
fi:

S:=Pols(x,d-1):
S union {seq(s+x^d, s in S)} mod 2:
end:

#PolsE(x,d): The set of polynomials of exactly degree d
# i.e. we only take polynomials that are degree d
PolsE:=proc(x,d) local S,s:
S:=Pols(x,d-1):
{seq(s+x^d, s in S)}:
end:



#IsIr(P): Is P irreducible mod 2?
IsIr:=proc(P) option remember:
evalb(Factor(P) mod 2 = P mod 2):
end:


#IsPr(P,x): is the polynomial P in x of degree d primitive?
IsPr:=proc(P,x) local d,m,i:
d:=degree(P,x):
m:=2^d-1:
for i from 1 to m-1 do
 if rem(x^i,P,x) mod 2=1 then
  return false:
 fi:
od:

if rem(x^m,P,x) mod 2<>1 then
 print(`Something bad happened, Cocks' article is wrong`):
 return(FAIL):
fi:
true:
end:


#WisW(d,x): inputs a pos. integer d and outputs the list of sets
# of polynomials of degree exactly d such that
# (i) neither
# (ii) irreducible but not primitive
# (iii) primitive but not irreducible
# (iv) both 
WisW:=proc(d,x) local S,s,Si,Sp,Sip,Sne:
#Si:=set of polys of degree d (mod 2) that are irreducible but not primitive
#Sp:=set of polys of degree d (mod 2) that are NOT irreducible but primitive
#Sip:=set of polys of degree d (mod 2) that are both irreducible and primitive
#Sne:=set of polys of degree d (mod 2) that are neither irreducible nor primitive
S:=PolsE(x,d):
Si:={}:
Sp:={}:
Sip:={}:
Sne:={}:

for s in S do
 if IsIr(s) and not IsPr(s,x) then
  Si:=Si union {s}:
 elif IsPr(s,x) and not IsIr(s) then
  Sp:=Sp union {s}:
 elif IsIr(s) and IsPr(s,x) then
  Sip:=Sip union {s}:
 else
  Sne:=Sne union {s}:
 fi:  
od:
[Sne,Si,Sp,Sip]:

end:


##################################################### from class:
Help16:=proc():
 print(`SR(INI,L,n),IsPer(L,t),FindPer(L)`):
end:

#SR=Shift Register from Clifford Cocks summary
#SR(INI,L,n): inputs a 0,1 list of length L[nops(L)],
# a list of increasing positive integers L, and 
# outputs the sequence of 0s and 1s generated by them of length n
# outputs the list of length n generated by the recurrence
# x[t]=x[t-L[1]]+..+x[t-L[nops(L)]]
SR:=proc(INI,L,n) local r,i,M,ng:
r:=nops(L):
if nops(INI)<>L[-1] then
 return(FAIL):
fi:

if not (convert(INI,set)={0,1} or convert(INI,set)={1}) then
 return(FAIL):
fi:

if not (type(L,list) and {seq(type(L[i],posint),i=1..nops(L))}={true} and sort(L)=L) then
 return(FAIL):
fi:

if not type(n,posint) then
 return(FAIL)
fi:

M:=INI:
while nops(M)<n do
 ng:=add(M[-L[i]],i=1..r) mod 2:
 M:=[op(M),ng]:
od:
end:

#IsPer(L,t): Is L periodic of period t?
IsPer:=proc(L,t) local i:
if {seq(L[i]-L[i+t],i=1..nops(L)-t)}={0} then
 true:
else
 false:
fi:
end:


#FindPer(L): finds the period of the list L
FindPer:=proc(L) local t:
for t from 1 to trunc(nops(L)/2) while not IsPer(L,t) do od:

if t=trunc(nops(L)/2)+1 then
RETURN(FAIL):
else
 RETURN(t):
fi:

end: