# hw21 - Pablo Blanco
# OK to post


CFx:=proc(x,k) local X, cfrac,tr,fr,i:
 if not(type(x,float)) or x<0 then:
   error "The first argument of Cfx should be type float and positive.":
 fi:

 if not type(k,posint) then:
   error "The second argument of Cfx should be type posint":
 fi:

 Digits:=2*k^2:              # I'm not honestly not sure what this should be. This is a guess.
 X:=convert(x,rational):
 cfrac:= Array([0$k]):
 tr:=trunc(X):
 fr:=frac(X):
 cfrac[1]:=tr:

 for i from 2 to k do:
  if fr = 0 then:
   break:
  fi:
  cfrac[i]:=1/fr:
  fr:=frac(cfrac[i]):
  tr:=trunc(cfrac[i]):
  cfrac[i]:=tr:

 od:
 cfrac:=convert(cfrac,list):

 [op(1..i-1, cfrac)]: 
end:

# inputs x such that either x is a posint or x^2 is a posint
# returns [L1,L2] where L1 and L2 are lists and L2 represents the
# periodic part of the continued fraction representation of x.
GuessPCF:=proc(x) local repeat,k,A,i,Pguess,j,conserv,N:
 if type(x,posint) then:
  return([[x],[]]):
 fi:

 kernelopts(assertlevel = 1):
 ASSERT(type(x*x,posint),"Input x must be a quadratic irrationality or a positive integer"):

 ### Assumption: periodic part only occurs starting at the fractional part. Is this true? 
 
 repeat:=false:
 conserv:=false: 

 k:=max(trunc(x),50):               # I chose this arbitrarily.
 A:=CFx(evalf(x,2*k^2),k): # long continued fraction
 N:=nops(A):
 i:=3:

 #if x=sqrt(37) then: return(A): fi: 


 while not repeat do: 
  # first, guess that L2 is a single digit:
  if nops({op(2..N,A)})=1 then:
   return([[A[1]],[A[2]]]):
  fi:

  # program has gone on for too long:
  if i > nops(A) then
   error "the procedure failed to make a guess":
  fi:  

  if A[2]<>A[i] then:
   i++:
   continue:
  fi:  

  Pguess:=[op(2..i-1,A)]:
  
  for j from 1 to N-i do:
   if A[i+j]<>Pguess[(j mod (i-2))+1] then:
    i++:
    conserv:=false:
    break:
   fi:
   conserv:=true:
  od:
  repeat:=conserv:

 od:

 return([[A[1]],Pguess ]):
end:

# 3. OEIS seq
# [seq(nops(GuessPCF(sqrt(i))[2]), i=1..100)];
# [0, 1, 2, 0, 1, 2, 4, 2, 0, 1, 2, 2, 5, 4, 2, 0, 1, 2, 6, 2, 6, 6, 4, 2, 0, 1, 2, 4, 5, 2, 8, 4, 4, 4, 2, 0, 1, 2, 2, 2, 3, 2, 10, 8, 6, 12, 4, 2, 0, 1, 2, 6, 5, 6, 4, 2, 6, 7, 6, 4, 11, 4, 2, 0, 1, 2, 10, 2, 8, 6, 8, 2, 7, 5, 4, 12, 6, 4, 4, 2, 0, 1, 2, 2, 5, 10, 2, 6, 5, 2, 8, 8, 10, 16, 4, 4, 11, 4, 2, 0]

# 4. e approx:
# 14665106/5394991
# using CFx(E,19); where E:=evalf(exp(1),30);

# Pi approx:
# 5419351/1725033
# using CFx(P,12); where P:=evalf(Pi,50);