#hw3.txt, Feb. 2 2014, Anthony Zaleski


#Packages we're using 
#(including C3.txt from class,which should be put 
#in the same folder as this file):

with(linalg):
with(combinat):
read("C3.txt"):

#Help()

Help:=proc() 
print('Help(),HilDiff(n,d),MinDigits(N,eps,digits),inv1(p),gf(n,q),gfGuess(n,q)');
end:

##########Problem 1: Garvan pp. 60-90##########

#Note: to see the output--which includes plots
#& animation, you must execute this.


#Differential Equations
f:='f': y:=f(x):
dsolve(diff(y,x)=y,y);

#2D Plots
plot(sin(1/x),x=0..Pi);
plot([t*cos(t),t*sin(t),t=0..8*Pi]);

with(plots):
polarplot(cos(10*t),t=0..2*Pi,numpoints=10000);
implicitplot(sin(x)+sin(y)=.5,x=-4*Pi..4*Pi,y=-4*Pi..4*Pi,numpoints=10000);

#3D Plots
x:='x': y:='y':
plot3d(sin(x^2+y^2),x=-Pi..Pi,y=-Pi..Pi);
plot3d([r*sin(t),r*cos(t),t],r=0..1,t=0..3*Pi);

spacecurve([cos(t),sin(t),t],t=0..6*Pi);

implicitplot3d(y^2-x^2=z,x=-2..2,y=-2..2,z=-2..2);

animate3d([r*sin(t+a),r*cos(t+a),t+a],r=0..1,t=0..3*Pi,a=0..6*Pi,title="Rotating helix animation.",frames=40);

#Linear Algebra 
with(linalg):
A:=matrix([[1,2],[3,4]]);
B:=matrix(2,2,[2,1,3,4]);
multiply(A,B);
A+B;
evalm(%);


##########Problem 2: Floating Point Error##########

#Here we record the number of digits d required to get
#evalf(det(Hil(n))) and det(evalf(Hil(n))) to agree.
#The following function outputs the descrepancy
#in terms of n and d.  

HilDiff:=proc(n,d):

Digits:=d:
abs(evalf(det(Hil(n)))-det(evalf(Hil(n)))):

end:

#Next, MinDigits takes a list N of n-values
#and outputs a list of the min number of 
#digits to get HilDiff<eps.  It checks all
#d in the digits list.

MinDigits:=proc(N,eps,digits) local k,l,n,d,L:

k:=nops(N):
l:=nops(digits):
L:=[]:
d:=1:

for n from 1 to k do
	while d<l and HilDiff(N[n],digits[d])>eps do
		d:=d+1:
	od:
	L:=[op(L),d]:
	print(L):
	d:=1:

od:

L:

end:


#Here is what we find:

#MinDigits(L,10^(-100),[seq(i,i=1..10)]);
#[10, 6, 4, 3, 1, 2, 1, 1, 1, 1] #these are the d-values

#It doesn't seem to make much sense.  There may 
#be a problem since keeping a small number of digits 
#misleadingly makes two numbers agree; for example, 
#if we set Digits:=1, then evalf(Pi-3)=0.

#I think we need better criteria for determining
#when the two computed values "agree."

##########Problem 3: Counting Inversions##########

#inv1(p) counts inversions in the list p.

inv1:=proc(p) local n,i,j,N:

n:=nops(p):
i:=1:
j:=1:
N:=0:

for i from 1 to n-1 do
	for j from i+1 to n do
		if p[i]>p[j] then
			N:=N+1:
		fi:
	od:
od:
N:
end:

#gf(n,q) sums q^inv1(p) over all perms p of length n.

gf:=proc(n,q) local P,p:

P:=permute(n):
add(q^inv1(p),p in P):

end:

#If we factor gf(n,x) for n from 1 to 8, we get
#the following:
(*
for n from 1 to 8 do factor(gf(n,x)); od;

1
1+x
(1+x)*(x^2+x+1)
(x^2+1)*(x^2+x+1)*(1+x)^2
(x^2+x+1)*(x^4+x^3+x^2+x+1)*(x^2+1)*(1+x)^2
(x^2-x+1)*(x^4+x^3+x^2+x+1)*(x^2+1)*(x^2+x+1)^2*(1+x)^3
(x^2-x+1)*(x^4+x^3+x^2+x+1)*(x^6+x^5+x^4+x^3+x^2+x+1)*(x^2+1)*(x^2+x+1)^2*(1+x)^3
(x^2-x+1)*(x^4+x^3+x^2+x+1)*(x^4+1)*(x^6+x^5+x^4+x^3+x^2+x+1)*(x^2+1)^2*(x^2+x+1)^2*(1+x)^4
*)

#There seems to be a pattern in the first 3 lines,
#but then things get weird; perhaps we are overfactoring
#and destroying the pattern.  Let's make a function
#gfGuess(n,q) that outputs the expected formula
#and compare this with the actual gf:

gfGuess:=proc(n,q) local i:
mul(sum(q^j,j=0..i),i=0..n-1):
end:


(*
for n from 1 to 8 do expand(gf(n,x)-gfGuess(n,x)); od;
0
0
0
0
0
0
0
0
*)

#It works!  Apparently our guess is indeed correct.
#As for rigorously proving this.....

#.....left to the reader!