#April 2, 2018, ExpMath (RU) C20.txt
Help:=proc(): print(` HS(p) , Rev(L) , HJS( p), SmithTable(N), SmithTableF(N), EvalCFg(A,B) `): end:
#HS(p): all the continued fractions of p/i for i from 2 to (p-1)/2
HS:=proc(p) local i:
{seq(RatToCF(p/i),i=2..(p-1)/2)}:
end:

Rev:=proc(L) local i:
[seq(L[nops(L)+1-i],i=1..nops(L))]:
end:

#HJS(p): finding [a,b] such that p=a^2+b^2 following Henry Smith
HJS:=proc(p) local S,i,s,a,b:

S:=convert(HS(p),list):

for i from 1 to nops(S)  while Rev(S[i])<>S[i] do od:
if i=nops(S)+1 then
  RETURN(FAIL):
else
s:=S[i]:
fi: 

if nops(s) mod 2=1 then
 RETURN(FAIL):
fi:

a:=P([seq(s[i],i=1..nops(s)/2)]): b:= P([seq(s[i],i=1..nops(s)/2 -1)]):

if p<>a^2+b^2 then
 print(`Fire Prof. Smith or fire me for messing up`):
  RETURN(FAIL):
fi:
[a,b]:

end:

#SmithTableF: all [p,[a,b]] where p=a^2+b^2 and p is 1 mod 4 for the first such primes 1 mod 4
SmithTableF:=proc(N) local L,p:
p:=5:
L:=[]:


while nops(L)<N do
  if p mod 4=1 then
   L:=[op(L), [p,Rev(SofS(p)[1])]]:
  fi:
   p:=nextprime(p):
od:
 L: 

end:



#SmithTable: all [p,[a,b]] where p=a^2+b^2 and p is 1 mod 4 for the first such primes 1 mod 4
SmithTable:=proc(N) local L,p:
p:=5:
L:=[]:


while nops(L)<N do
  if p mod 4=1 then
   L:=[op(L), [p,HJS(p)]]:
  fi:
   p:=nextprime(p):
od:
 L: 

end:


#Corrected after class
#EvalCFg(A,B): evaluates the non-simple continued with numerators B and denominators A s given in Chrystal's book p. 512, Eq. (1)
EvalCFg:=proc(A,B) local A1,B1,x:
if not (type(A,list) and nops(A)>0) then
  RETURN(FAIL):
fi:

if not (type(B,list) and nops(B)>0) then
  RETURN(FAIL):
fi:

if nops(A)<>nops(B) then
  RETURN(FAIL):
fi:

if nops(A)=1 then
  RETURN(B[1]/A[1]):
fi:

A1:=[op(2..nops(A),A)]:

B1:=[op(2..nops(B),B)]:

x:=EvalCFg(A1,B1):

normal(B[1]/(A[1]+x)):
end:


#end new stuff


#Start from C19.txt: 
Help19:=proc(): print(` SofS(n) , P(L), Q(L), FibS(n) , wt(s,a) `): end:

#SofS(n): the set of all the pairs [a,b], a<=b,  such that a^2+b^2=n
SofS:=proc(n) local S,a,b:

S:={}:

for a from 1 to trunc(sqrt(n))+1 do
  b:=sqrt(n-a^2):
  if type(b,integer)  and a<=b then
    S:=S union {[a,b]}:
  fi:
od:

S:

end:



#The numerator of the (simple) continued fraction L
P:=proc(L) option remember:

if L=[] then
 1:
elif nops(L)=1 then
  L[1]:
else
expand(L[1]*P([op(2..nops(L),L)])+ P([op(3..nops(L),L)])):

fi:
 end:


Q:=proc(L):
P([op(2..nops(L),L)]):
end:

#FibS(n): the set of vectors in {1,2} that add-up to n
FibS:=proc(n) local S1,S2,s:
option remember:
if n=0 then
 RETURN({[]}):
elif n=1 then
  RETURN({[1]}):
else
 S1:=FibS(n-1):
 S2:=FibS(n-2):
  
 {seq([1,op(s)],s in S1), seq([2,op(s)],s in S2)}:
fi:
end:

#wt(s,a): the weight of the Fibonacci word s, product of a[j] over all intermediate stops
#j that follow immediately j-1
#wt(122212)=a[1]*a[8] wt(22222)=1
wt:=proc(s,a) local i,p:
p:=1:
for i from 1 to nops(s) do
 if s[i]=1 then
  p:=p*a[add(s[j],j=1..i)]:
 fi:
od:

p:

end:

#CheckQI(L,k): Checks the identity P([a1,...., an])=P([a1...,ak])*P([a(k+1)... an])+  ... 
CheckQI:=proc(L,k) local i:
expand(P(L)-P([op(1..k,L)])*P([op(k+1..nops(L),L)])-
P([op(1..k-1,L)])*P([op(k+2..nops(L),L)])):
end:

#end C19.txt




###from C18.txt
#C18.txt March 26, 2018, ExpMath (RU)
Help18:=proc(): print(`RatToCF(x),NuToCF(x,n) , EvalCF(L), SymCF(a,n) `):end:

#RatToCF(x): converts a pos. rational number x into a continued fraction
RatToCF:=proc(x) local a,x1:
 if not (x>0 and type(x,rational)) then
  RETURN(FAIL):
 fi:

a:=trunc(x):

if a=x then
  RETURN([a]):
fi:

x1:=1/(x-a):

[a, op(RatToCF(x1))]:


end:


#NuToCF(x,n): appx. a pos. number x into a continued fraction with <=n terms
NuToCF:=proc(x,n) local a,x1:

 
 if not (x>0 and type(n,integer) and n>=1) then
  RETURN(FAIL):
 fi:


a:=trunc(x):

if n=1 then
  RETURN([a]):
fi:
  
if x=a then
  RETURN([a]):
fi:

x1:=1/(x-a):

[a,op(NuToCF(x1,n-1))]:
end:

#EvalCF(L): evaluates the simple continued 
EvalCF:=proc(L) local L1,x:
if not (type(L,list) and nops(L)>0) then
  RETURN(FAIL):
fi:

if nops(L)=1 then
 RETURN(L[1]):
fi:

L1:=[op(2..nops(L),L)]:

x:=EvalCF(L1):

normal(L[1]+1/x):
end:

#SymCF: the rational function of a[1], ..., a[n] describing EvalCF([a[1], ..., a[n]]):
#the numerator followed by the denominator
SymCF:=proc(a,n) local i,yonah:
yonah:=EvalCF([seq(a[i],i=1..n)]):

numer(yonah),denom(yonah):
end:



####