#Version of 12:00noon (check for updates)

#C24.txt: April 16, 2018, Math 640 (Spring 2018), Rutgers University, Dr. Z.
#Homemade linear algebra
#Rather than use Maple's linear algebra package, it is more instructive to write our own
#Today's virtual class will consist in understanding and testing the Maple code below
#we also use linalg[det] (the home-made version sometimes FAILs)
#At the top you ave Yukun Yao's homework hw22.txt, illustrating the amazing Iterated polygon phenomenon
#Try, e.g. IterPolygon(20,N), for N=10,20,30, ...,  several times (of course, the way it is programmed now
#the intial polygon changes every time, you may want to write a version that
#plots successive iterations, and even animate them).  Notice that sometimes it takes
#a fairly long time until it reaches the ellipse (i.e. oval) shape, and that slant is alway close to 45 degrees
#but sometimes pointing left and sometimes ponting right, and the "fatness" (eccentricity) changes
print(`Welcome to virtual class C24`):
print(`For a list of the new procedures, type: Help(); `):
print(`For a list of the procedures from h22.txt (thanks to Yukun Yao)  type: HelpYY(); `):


Help:=proc(): print(` MtV(M,V),  LCdet(M), RandM(m,n,N), ChM(M,x), ChP(M,x), EV(M)    `): end:

HelpYY:=proc(): print(`IterPolygon(n,N)`): end:

#######From Yukun Yao's homework
#OK to post homework
#Yukun Yao, Apr 14, 2018, Assignment 22

##Problem 2## For the random polygon project, get a feel for it, by
#Writing a procedure
#RandomGoodPoly(n)
#that inputs a positive integer and outputs a random polygon with n vertices, whose centroid is the origin.
#Hint: generate two random vectors, x and y, and substract from them their averages, also call these new vectors x and y to make their sum to zero and let the polygon be [[x1,y1], ..., [xn,yn],[x1,y1]]

#Write a procedure
#AvePolygon(P)
#that inputs a polygon, (given by its list of vertices), and outputs the polygon whose vertices are the midpoints of the edges.

#Write a procedure
#IterPolygon(n,N)
#that inputs a RandomGoodPoly(n) and a large pos. integer N, and performs AvePolygon(P) N times, and also plots it (use the command plot(P)). Perform IterPolygon(20,200) many times and look at the picture.

#To generate a random polygon usually the x and y are rational. So by scaling let x,y be integers and give a range for them. 
RandomGoodPoly:=proc(n) local i,S,Ave,ra,L:
S:=[]:
ra:=rand(0..1000):
for i from 1 to n do
  S:=[op(S), [ra(),ra()]]:
od:
Ave:=convert(S,`+`)/n:
L:=[Ave$n]:
S:=S-L:
S:=[op(S), S[1]]:
end:


AvePolygon:=proc(P) local n,S,i:
n:=nops(P)-1:
S:=[]:
for i from 1 to n do
  S:=[op(S),(P[i]+P[i+1])/2]:
od:
S:=[op(S), S[1]]:
end:


IterPolygon:=proc(n,N) local P,i:
P:=RandomGoodPoly(n):
for i from 1 to N do 
  P:=AvePolygon(P):
od:
plot(P):
end:

#It seems that the picture will converge to an oval.
#######End From Yukun Yao's homework

#start today's class
#Our date structure for vectors will be lists;  For matrices, lists of lists, we will NOT use Maple's linear algebra
#packages EXCEPT the solve command

#MtV(M,V): inputs a matrix M and a column vector V (expressed as a list), outputs the vector MV
#For example, try
#MtV([[1,2],[3,4]],[5,6]); MtV([[a11,a12],[a21,a22]],[b1,b2]);
MtV:=proc(M,V) local i,j:

if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and type(V,list) and
         nops(M)>0 and nops(M[1])=nops(V)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

[seq(add(M[i][j]*V[j],j=1..nops(M[i])),i=1..nops(M))]:

end:

#LCdet(M): the determinant of a matrix using the Rev. Dodgson's method (as nicely proved in
# http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/twotime.html
LCdet:=proc(M) local M11,M12,M21,M22,Mid,i,n:

if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and nops(M[1])=nops(M)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

n:=nops(M):
if n=1 then
  RETURN(M[1][1]):
elif n=2 then
  RETURN(M[1][1]*M[2][2]-M[1][2]*M[2][1]):
else

 M11:=[seq([op(1..n-1,M[i])],i=1..n-1)]:
 M12:=[seq([op(2..n,M[i])],i=1..n-1)]:

 M21:=[seq([op(1..n-1,M[i])],i=2..n)]:
 M22:=[seq([op(2..n,M[i])],i=2..n)]:

 Mid:=[seq([op(2..n-1,M[i])],i=2..n-1)]:

 if CLdet(Mid)=0 then
  print(`Thou should not divide by zero`):
   RETURN(FAIL):
  else
 RETURN(normal((LCdet(M11)*LCdet(M22)-LCdet(M12)*LCdet(M21))/LCdet(Mid))):
 fi:
fi:
 
end:

#RandM(m,n,N): a random m by n matrix with entries from 1 to N
RandM:=proc(m,n,N) local ra,i,j:
ra:=rand(1..N):
[seq([seq(ra(),i=1..n)],j=1..m)]:
end:



#ChM(M,x): the characteristic matrix of the square matrix M, w.r.t the variable x, i.e. M-Ix
ChM:=proc(M,x) local n,i,j:
if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and nops(M[1])=nops(M)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

n:=nops(M):

[seq([seq(M[i][j],j=1..i-1),M[i][i]-x,seq(M[i][j],j=i+1..n)],i=1..n)]:
end:

#ChP(M,x): the characteristic polynomial of the square matrix M, w.r.t the variable x, i.e. M-Ix
ChP:=proc(M,x)
linalg[det](ChM(M,x)):
end:




#EV(M): the list of eigenvalues of the square matrix M, and decreasing order of their absolute values. Exact version
#Warning, sorting is very naive
EV:=proc(M) local P,x,L:

if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and nops(M[1])=nops(M)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:
P:=ChP(M,x):

L:=[solve(P,x)]:


end:


##