#OK to post homework
#Jingze Li(RUID:195000313),04/21/19
#############READ ME###########
#the homework 23 consists of 2 parts, corresponding to 1 to 2 problems of hw23.

################
Help:=proc() print(`Fc(n,K),Tc(n,K)`): end:


#######  1 ##########

with(LinearAlgebra):

# another algorithm to compute the n_th fibonacci number.
#using matrix iteration. 
# input integer n and K, output the n th fibonacci number mod K.


Fc:=proc(n,K)
local A,B,C:
 
A:=Matrix([[1],[1]]):
if n = 1 or n=2 then
  1:
else
B:= Matrix([[1,1],[1,0]]):
C:=B^(n-2).A mod K:
C[1]:
fi:
end:

# a test  like  time(Fc(10^6 ,2003)) and time(F(10^6,2003))   shows that this algorithm is faster than the best algorithm we write in class. However, when n is larger than 10^10, it takes very long time. we are unable to compute the 10^20 th fibonacci number

# a possible way to reduce the calculation effort in the matrix power, is to diagonalize the iteration matrix, so the computation only concerns the power of the eigenvalues of the iteration matrix.

# this will make calculating 10^20 th fibonacci number accessible.
# however, since the eigenvalues (1+-sqrt(5))/2 are irrational numbers, raising them to n_th power has to use the evalf function, I am afraid that the truncation error can be large.




# algorithm to compute the n_th tribonacci number. Similar to Fc(n,K)
#using matrix iteration. 
# input integer n and K, output the n th tribonacci number mod K.



Tc := proc(n,K)
local    A,B,C:

A:=Matrix([[1],[1],[0]]):
if n=1 then
 0:
elif n=2 or n=3 then 
 1:
else
B:=Matrix([[1,1,1],[1,0,0],[0,1,0]]):
C:=B^(n-3).A mod K:
C[1]:
fi:
end:


# still, this algorithm is not fast enough to compute Tc(10^20,2003). A possible way, like the previous one, is to diagonalize the iteration matrix to  simplify the matrix power camputation.








#C23.txt, April 18, 2019, Computing derivatives and gradients of complicated functions
Help:=proc(): print(`Box(n), IsGG(v) , Fsss(n) , Fss(n) , Fs(n), F(n,K) , EvalC(L,x,x0)`): 
print(`RP(x,d,K), RC(x,d,K,n) , EvalCD(L,x,x0) `):
end:

#Box(n): the set of 0-1 vectors of lenght n
Box:=proc(n) local S,s:
if n=0 then
  RETURN({[]}):
else 
 S:=Box(n-1):
 
RETURN({seq([op(s),0], s in S),seq([op(s),1], s in S)}):
fi:
end:

#IsGG(v): Does v contain 11
IsGG:=proc(v) local i:

for i from 1 to nops(v)-1 do
  if v[i]=1 and v[i+1]=1 then
    RETURN(false):
  fi:
od:
true:
end:
 
#Fsss(n):the n-th Fibonacci number using the extremely stupid way
#One of the Combinatorial definition of fibonacci number, the number
#of 0-1 vectors of length n avoiding 11
Fsss:=proc(n) local S,G,co,s:
S:=Box(n):
co:=0:

for s in S do
 if IsGG(s) then
  co:=co+1:
 fi:
od:
co:
end:


#Fss(n): the n-th Fibonacci number using the very stupid way
Fss:=proc(n):
if n=0 or n=1 then
 n:
else
 Fss(n-1)+Fss(n-2):
fi:
end:

#Fs(n): the n-th Fibonacci number using the usual way mod K
Fs:=proc(n,K)
option remember:
if n=0 or n=1 then
 n:
else
 Fs(n-1,K)+Fs(n-2,K) mod K:
fi:
end:

#F(n,K): the n-th Fibonacci number the clever way mod K
F:=proc(n,K) local a,b,c,i:

a:=0:
b:=1:

for i from 1 to n do
c:=a+b mod K:
a:=b:
b:=c:
od:
a:
end:

#RP(x,d,K): a random poly in x of degree k with coefficients from 1 to K do
RP:=proc(x,d,K) local ra,i:
ra:=rand(1..K):

add(ra()*x^i,i=0..d):

end:
#RC(x,d,K,n): a random chain of length n of poly in x of degree k with coefficients from 1 to K do
RC:=proc(x,d,K,n) local i:
[seq(RP(x,d,K),i=1..n)]:
end:

 #[x0->x1->x2-> ... ->xn] 

  #EvalC(L,x,x0): inputs a listL of expressions in the variable x and a number x0
   #outputs L[n](L[n-1](.... L[1](x0))))))))) mod K:
  #x0->L[1](x0)->L[2](L[1](x0))-....
 EvalC:=proc(L,x,x0) local a,i:
  a:=x0:
   for i from 1 to nops(L) do
     a:=subs(x=a,L[i]):
   od:
a:
end:


 #EvalCD(L,x,x0): inputs a listL of expressions in the variable x and a number x0
   #outputs L[n](L[n-1](.... L[1](x0))))))))) followed by the derivative of the composition
   #at x0:
  #x0->L[1](x0)->L[2](L[1](x0))-....
  #f(g(x))'= f'(g(x))*g'(x) . (F[i](Previous(x)))'(x0)= F[i]'(Previoius(x)|x=x0)* (Previous'(x)|x=x0)
 EvalCD:=proc(L,x,x0) local a,b,i:
  a:=x0:
   b:=1:
   for i from 1 to nops(L) do
      b:=b*subs(x=a,diff(L[i],x)) :
     a:=subs(x=a,L[i]):
     
   od:
a,b:
end: