#OK to post homework
#Lucy Martinez, 04-02-2024, Assignment 20


# Problem 1: Using Maple's graphing command, write a procedure Diag(a,b)
# that draws a picture like Fig. 3 
# (with the points A, B,C, A+B labeled and indicated on the diagram)
# for a general elliptic curve y^2=x^3+a*x+b , for general a, and b,
# and A and B picked using RP(a,b)
Diag:=proc(a,b):

end:




# Problem 2: Write a procedure IntSol(a,b,K)
# that finds all integer solutions [x,y] with -K<=x<=K of the diophantine equation
# y^2=x^3+a*x+b
IntSol:=proc(a,b,K) local S,y1,y2,i:
S:=[]:
for i from -K to K do
 y1:=sqrt(i^3+a*i+b):
 y2:=-sqrt(i^3+a*i+b):
 if type(y1,integer) then
  S:=[op(S),[i,y1]]:
 elif type(y2,integer) then
  S:=[op(S),[i,y2]]:
 fi:
od:
S:
end:



# Problem 3: Write a procedure SolChain(a,b,S1,S2,K)
# that starts out with two solutions S1 and S2 of y^2=x^3+a*x+b
# and repeatedly uses AddE(a,b,L[-2],L[-1]) to get a list of length K+2 of solutions
# in rational numbers. Can you get infinitely solutions this way?
SolChain:=proc(a,b,S1,S2,K) local L,i:
if not (type(S1,list) and type(S2,list) and nops(S1)=2 and nops(S2)=2) then
 print(`Check that S1 and S2 are lists of 2 elements`):
 return(FAIL):
fi:
L:=[S1,S2]:
for i from 1 to K do 
 L:=[op(L),AddE(a,b,L[-2],L[-1])]:
od:
L:
end:



####################################### C20.txt
Help20:=proc():
 print(`RP(a,b),DoesLie(a,b,P),AddE(a,b,P,Q)`):
end:


#Random Point on the elliptic curve y^2=x^3+a*x+b
#RP(a,b): Generates a random point on the elliptic curve y^2=x^3+a*x+b
RP:=proc(a,b) local x:
 x:=rand(0..200)()/100:
 [x,sqrt(x^3+a*x+b)]:
end:


#DoesLie(a,b,P): Given P=(x,y), outputs y^2-(x^3+a*x+b)
# if P lies on the curve, then it returns 0
DoesLie:=proc(a,b,P):
 P[2]^2-P[1]^3-a*P[1]-b:
end:


#AddE(a,b,P,Q): the "addition" (over real elliptic curves) of the point P and Q
# that lie on the elliptic curve   y2=x3+ ax+b ?
AddE:=proc(a,b,P,Q) local y,s,eq,sol,xR,R:
s:=(P[2]-Q[2])/(P[1]-Q[1]):
#equation is y-P[2]=(x-P[1])*s
y:=P[2]+(x-P[1])*s:
sol:={solve(y^2-x^3-a*x-b,x)}:

if not (member(P[1],sol) and member(Q[1],sol) and nops(sol)=3) then 
 print(`Something bad happened`):
 return(FAIL):
fi:
sol:=sol minus {P[1],Q[1]}:
xR:=sol[1]:
[xR,sqrt(xR^3+a*xR+b)]:
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: