Help:=proc(): print(`RandM(m,n,N),IterPolygon(n,N), IterPolygon1(P,N), RandomGoodPoly(n), A(n), Nor(v), EVC(M) , Champs(L) `): end:


with(linalg):

#######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:



IterPolygon1:=proc(P,N) local i,P1:
P1:=P:
for i from 1 to N do 
  P1:=AvePolygon(P1):
od:
plot(P1):
end:






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


#From Edna Jones' hw24
#A(n): Inputs a positive integer n, and outputs the n by n matrix (given as a list of lists), corresponding to the "averaging process" [x1,x2, ..., xn] ->[(x1+x2)/2, (x2+x3)/2, ...., (xn+x1)/2]
A:=proc(n) local j:
    if n=1 then
        return [[1]];
    else
        return [seq([0$(j-1), 1/2, 1/2, 0$(n-j-1)], j=1..(n-1)), [1/2, 0$(n-2), 1/2]];
    fi:
end:
#End from Edna Jones' hw24

#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:

Nor:=proc(v) local i,a:
a:=evalf(add(abs(v[i])^2,i=1..nops(v))):
if a<10^(-Digits+2) then
  RETURN(FAIL):
fi:
[seq(v[i]/a,i=1..nops(v))]:
end:


#EVC(M): the list of [eigenvalue,eigenvector] in weakly decreasing order of the absolute value of the eigenvalues
EVC:=proc(M)  local i,L:
L:=[eigenvectors(M)]:


if {seq(L[i][2],i=1..nops(L))}<>{1} then
  print(`Some eigenspaces are not one-dimensional`):
  RETURN(FAIL):
fi:

[seq([L[i][1],Nor(  convert(L[i][3][1],list) )],i=1..nops(L))]:

end:



#Champs(L): given a list of pairs outputs the set indices of pairs with the largest first component  according to absolute value
Champs:=proc(L) local cha,rec,i:
if nops(L)=0 then
 RETURN(FAIL):
fi:

cha:={1}:
rec:=abs(L[1][1]):

for i from  2 to nops(L) do
 if abs(L[i][1])>rec then
   cha:={i}:
   rec:= abs(L[i][1]):
  elif abs(L[i][1])=rec then
    cha:=cha union {i}:
fi:
od:
cha:
end:


##