#OK to post homework
#Joseph Koutsoutis, 03-10-2024, Assignment 15


read `C14.txt`:
read `C15.txt`:
with(LinearAlgebra):

Help := proc(): print(` qS(x,y,q), qOurPF(R,t,q), qaSr(a,q), EncodeBCH(d,q,x), DecodeBCH(d,q,x)  `): end:

#1 

# First we'll show that lim t->1/y[i] ((1-t*y[i])*A(t)/B(t)) = -y[i]*A(1/y[i])/B'(1/y[i]).
# Since lim t->1/y[i] (1-t*y[i]) = 0 and
# lim t->1/y[i] (B(t)) = 0, we may apply L'Hopital's rule to obtain:

# lim t->1/y[i] ((1-t*y[i])*A(t)/B(t)) = lim t->1/y[i] (((1-t*y[i])*A'(t)-y[i]*A(t))/B'(t))
#                                       = -y[i]*A(1/y[i])/B'(1/y[i]).

# Since A(t)/B(t) = C[1]/(1-t*y[1])+ ... + C[n]/(1-t*y[n]), we have that

# (1-t*y[i])*A(t)/B(t) = (1-t*y[i])*C[1]/(1-t*y[1])+ ... + (1-t*y[i])*C[n]/(1-t*y[n]).

# Hence -y[i]*A(1/y[i])/B'(1/y[i]) = lim t->1/y[i] ((1-t*y[i])*A(t)/B(t))
#                                  = lim t->1/y[i] ((1-t*y[i])*C[1]/(1-t*y[1])+ ... + (1-t*y[i])*C[n]/(1-t*y[n]))
#                                  = C[i].


#2

# To check Ramanujan's example on his seminal paper I ran the following code:

# x := [(8-sqrt(5))/(2*sqrt(5)),-3/5,(18-sqrt(5))/10,-(8+sqrt(5))/(2*sqrt(5)),(18+sqrt(5))/10]:
# y := [-(sqrt(5)+1)/(2),-1,(3-sqrt(5))/2,(sqrt(5)-1)/2,(3+sqrt(5))/2]:
# a := S(x,y):
# simplify(aSr(a)- [x,y]);

# I also changed one line in OurPF to this:

OurPF:=proc(R,t) local n,Top,Bot,zer,Y,i:
Top:=numer(R):
Bot:=denom(R):
zer:=[solve(Bot,t)]:
n:=nops(zer):
Y:=sort([seq(1/zer[i],i=1..nops(zer))], (a,b) -> evalf(a) < evalf(b)): # I changed this line
[[seq(-Y[i]*subs(t=1/Y[i],Top)/subs(t=1/Y[i],diff(Bot,t)),i=1..n)],Y]:
end:

# since Maple appears to sort expressions with radicals in an unexpected way.
# For example, sort([sqrt(5)/2-1/2,sqrt(5)/2+3/2]) outputs [sqrt(5)/2+3/2,sqrt(5)/2-1/2].


#3

qS := proc(x,y,q) local n,i,r:
if not (type(x,list) and type(y,list) and nops(x)=nops(y)) then
 RETURN(FAIL):
fi:
n:=nops(x):
[ seq(add(x[i]*y[i]^r,i=1..n) mod q,r=0..2*n-1)]:
end:


qOurPF:=proc(R,t,q) local n,Top,Bot,zer,Y,i:
Top:=numer(R):
Bot:=denom(R):
zer:=[msolve(Bot,q)]:
n:=nops(zer):
Y:=sort([seq(subs(zer[i], 1/t) mod q,i=1..nops(zer))]):
[[seq(-Y[i]*subs(t=1/Y[i],Top)/subs(t=1/Y[i],diff(Bot,t)) mod q,i=1..n)],Y]:
end:


qaSr:=proc(a,q) local n,A,B,i,Top,Bot,f,eq,var,t:
if not( type(a,list) and nops(a) mod 2=0) then
 RETURN(FAIL):
fi:
n:=nops(a)/2:

Top:=add(A[i+1]*t^i,i=0..n-1):

Bot:=1+add(B[i]*t^i,i=1..n):

f:=expand(Top-Bot*add(a[i]*t^(i-1),i=1..2*n)) mod q:
eq:={seq(coeff(f,t,i),i=0..2*n-1)}:
var:={seq(A[i],i=1..n),seq(B[i],i=1..n)}:
var:=msolve(eq,q):
qOurPF(normal(subs(var,Top)/subs(var,Bot)) mod q,t, q):
end:

#4

EncodeBCH := proc(d,q,x) local S,eqs,i,j,n,y:
 n := nops(x) + d - 1:
 eqs := {seq(add(i^j * x[i], i=1..nops(x)) + add(i^j * y[i], i=(nops(x)+1)..n), j=0..(d-2))}:
 S := msolve(eqs, q):
 subs(S, [op(x), seq(y[i], i=(nops(x)+1)..n)]):
end:



#5


UpperTriangular := proc(mat117,aug117,q) local t,curr_row,i,j,k,tmp,M,aug:
 M := mat117:
 aug := aug117:
 t := nops(M):
 curr_row := 1:
 for j from 1 to t do:
  for i from curr_row to t do:
   if M[i][j] <> 0 then:
    tmp := M[i]:
    M[i] := M[curr_row]:
    M[curr_row] := tmp:

    tmp := aug[i]:
    aug[i] := aug[curr_row]:
    aug[curr_row] := tmp:
    
    for k from curr_row+1 to t do:
     aug[k] := aug[k] - aug[curr_row] * M[k][j] / M[curr_row][j] mod q:
     M[k] := M[k] - M[curr_row] * M[k][j] / M[curr_row][j] mod q:
    od:
    curr_row += 1:
   fi:
  od:
 od:
 M,aug:
end:


# DecodeBCH should be able to correct up to t errors. If more than t errors have occurred, this function
# will either return FAIL or some valid member of the code.
DecodeBCH := proc(d,q,x) local t,a,i,j,r,n,e,matrix117,aug117,eqs117,A,B,locator,evaluator,theta,errs,X,Y,y,zer:
 t := (d-1) / 2: 
 n := nops(x):
 a := [seq(add(i^r * x[i], i=1..n) mod q, r=0..(2*t-1))]:
 
 matrix117 := [seq([seq(a[t-i+j], i=1..t)], j=1..t)]:
 aug117 := [seq(-a[t+j] mod q, j=1..t)]:
 matrix117, aug117 := UpperTriangular(matrix117, aug117, q):
 e := Rank(matrix117):
 eqs117 := {seq(aug117[i] = add(matrix117[i][j] * B[j], j=1..e), i=1..e)}:

 locator := 1 + subs(msolve(eqs117, q), add(B[i] * theta^i, i=1..e)):

 for i from 1 to e do:
  A[i] := add(a[j] * coeff(locator, theta, i-j), j=1..i) mod q:
 od:
 evaluator := add(A[i] * theta^(i-1), i=1..e):

 zer := [msolve(locator, q)]:
 for y in zer do:
  if not type(subs(y, theta),integer) or y mod q = 0 then:
   return(FAIL):
  fi:
 od:

 Y := sort([seq(subs(zer[i], 1/theta) mod q,i=1..nops(zer))]):

 for y in Y do:
  if y > n then:
   return(FAIL):
  fi:
 od:

 X := [seq(-Y[i]*subs(theta=1/Y[i] mod q,evaluator)/subs(theta=1/Y[i] mod q,diff(locator,theta)),i=1..nops(Y))]:

 for y in X do:
  if denom(y) mod q = 0 then:
   return(FAIL):
  fi:
 od:

 X := X mod q:
 y := x + add(-[0$(Y[i]-1), X[i], 0$(n-Y[i])], i=1..nops(X)) mod q:
 if EncodeBCH(d,q,y[1..n+1-d]) = y then:
  return(y):
 else:
  return(FAIL):
 fi:
end: 




# TestDecodeBCH(n,d,q,K) tests DecodeBCH by creating K vectors x of size n, adding up to t errors, and
# checks to see if DecodeBCH(d,q,x + errors) = x.
# This function also calls DecodeBCH using these K vectors with up to d errors added just to
# make sure our code above does not fail in an unexpected way.
TestDecodeBCH:=proc(n, d, q, K) local i,j,x,y,rand_noise,rand_index,t:
 t := (d-1)/2:
 rand_noise := rand(1..q-1):
 rand_index := rand(1..n):
 
 for i from 1 to K do:
  x := EncodeBCH(d,q,RV(q,n+1-d)):
  y := x:
  for j from 1 to t do:
   y[rand_index()] += rand_noise():
  od:
  y := y mod q:
  
  if x <> DecodeBCH(d,q,y) then:
   return(FAIL):
  fi:

  y := x:
  for j from 1 to d do:
   y[rand_index()] += rand_noise():
  od:
  y := y mod q:
  DecodeBCH(d,q,y):
 od:

 return(SUCCESS):
end: