#March 21, 2016, finishing Jacobian conjecture

Help:=proc(): print(` ExtractC(P,x) , Vecs1(d,r), Vecs(d,r), GenPol(x,d,a) `): end:

##old stuff
#C14.txt: March 7, 2016, JC (cont'd)
Help14:=proc(): print(` ChopkT(P,x,i), Jac(P,x), RP(x,n,K) `): 
print(` Elem(P,C,i,j,r) `):
end:

with(linalg):

###Code from hw13.txt (thanks to Emily Kukura)#OK to post

Help13:=proc() print(`Hw13:`,`Chop3(P,x,y,z,i)`,`Inv3(P,x,y,z,N)` ,`Comp(P,Q,x)`,`Chopk(P,x,i)`,`Invk(P,x,N)`): end:



#Comp(P,Q,x): Inputs polynomial transformations 
#P and Q in the list of variables x=[x1,...,xn]
#outputs its composition P(Q(x))
Comp:=proc(P,Q,x) local i,n:
n:=nops(x):
if nops(Q)<>n then
  RETURN(FAIL):
fi:

expand(subs({seq(x[i]=Q[i],i=1..n)},P)):

end:


#3. 

#Chopk(P,x,i): inputs a polynomial P in the list of variables x and outputs
#the all the terms of total degree <=i
Chopk:=proc(P,x,i) local n,j,Q,P2,x2,xNew,k:

n:=nops(x):
x2:=x:

Q:=0:

if n>1 then
	for j from 0 to i do:
		P2:=coeff(P,x2[1],j):
		xNew:= [seq(x2[k+1], k = 1 .. n-1)]: #new variables, omitting x2[1]:
		Q:=Q+Chopk(P2,xNew,i-j)*x2[1]^j:
	od:
else #if n=1
	for j from 0 to i do:
		P2:=coeff(P,x2[1],j):
		Q:=Q+P2*x2[1]^j:
	od:
fi:

Q:=expand(Q):

end:


#Invk(P,x,N): inputs a list, P, of length, k, say, a list x of variables (of the same length), say x=[x1,x2,.., xk], and the members of P are each polynomials in the variables [x1,..., xk], and outputs a list of length k that consists of the beginning (up to degree N) of the k components of the inverse transformation

Invk:=proc(P,x,N) local n,L,Q,F,Qa,X,var,A,K,i,j,setVars,Eqns,S,C,C2,f,q,qa,l,
v,subsEqs,subsEqs2,varSubsEqs,varSubsEqsSet, eqn, se, vse:

f:='f': q:='q': qa:='qa': l:='l': v:='v':
n:=nops(x):
F:=[seq(f[i], i = 1 .. n)]:
Q:=[seq(q[i], i = 1 .. n)]:
Qa:=[seq(qa[i], i = 1 .. n)]:
L:=[seq(l[i], i = 1 .. n)]:
X:=[seq(v[i], i = 1 .. n)]:
A:=[]:

for i from 1 to n do:
	F[i]:=P[i]:
od:

for i from 1 to n do:
	if Chopk(F[i],x,0)<>0 then
		RETURN(FAIL):
	fi:
od:

for i from 1 to n do:
	L[i]:=Chopk(F[i],x,1):
	Q[i]:=F[i]-L[i]:
od:

Eqns:=[seq(eqn[i], i = 1 .. n)]:
for i from 1 to n do:
	Eqns[i]:=X[i]=L[i]+Qa[i]:
od:

setVars:=convert(x,set):

var:=solve(Eqns,setVars):

for i from 1 to n do:
	var:=subs(Qa[i]=Q[i],var):
od:

for i from 1 to n do:
	A:=[op(A),subs(var,x[i])]:
od:

K:=[0$n]:

subsEqs := [seq(se[i], i = 1 .. n)]:

for i from 1 to N do
	for j from 1 to n do:
		subsEqs[j] := x[j] = K[j]:
	od:
	subsEqs2 := convert(subsEqs, set):

	K:=subs(subsEqs2,A):

	for j from 1 to n do
		K[j]:=Chopk(K[j],X,i):
	od:

od:

varSubsEqs:=[seq(vse[i],i=1..n)]:
for j from 1 to n do
	varSubsEqs[j]:=X[j]=x[j]:
od:	

varSubsEqsSet:=convert(varSubsEqs,set):
S:=subs(varSubsEqsSet,K):

C:=Comp(S,P,x):
C2 := [seq(Chopk(C[i], x, N), i = 1 .. n)]:
if C2<>x then
	RETURN(FAIL):
else 
	print(`Comp(P,Q,x) routine has verified that this is the correct inverse transformation`):
fi:

S:

end:

###End of Emily's code


#ChopkT(P,x,i): inputs a polynomial transformation
#P  (given as a list of polyniamials in the list of variables x and outputs
#the list of all the chopped terms of total degree <=i
ChopkT:=proc(P,x,i) local j:
[ seq(Chopk(P[j],x,i),j=1..nops(P))]:
end:

Jac:=proc(P,x) local i,j:

det( [ seq([seq( diff(P[i],x[j]), j=1..nops(x))] , i=1..nops(P))]):

end:

#RP1(x,n,K): inputs a list of variables x and a pos. integer n
#outputs a random polynomial of degree n in the x, whose coefficients
#are from -K to K
RP1:=proc(x,n,K) local ra,k,xk,i,xB:
ra:=rand(-K..K):

k:=nops(x):

xk:=x[k]:

if k=1 then
 RETURN(add(ra()*xk^i,i=0..n)):
fi:

xB:=x[1..k-1]:

expand(add(RP1(xB,n-i,K)*xk^i, i=0..n)):

end:



#RP(x,n,K): inputs a list of variables x and a pos. integer n
#outputs a random polynomial of degree n in the x, whose coefficients
#are from -K to K WITHOUT a constant term
RP:=proc(x,n,K) local P:
P:=RP1(x,n,K):
P-Chopk(P,x,0):
end:

#RPT(x,n,K): a random transformation, of degree n in the list of variables x,
# from nops(x)-dim space to itself
RPT:=proc(x,n,K) local i:


[seq(RP(x,n,K),i=1..nops(x))]:

end:

#Elem(P,C,i,j,r): inputs a poly. trans. P, a constant C, and i,j (from 1 to nops(P)
#outputs the trans. obtained by replacing P[i] by P[i]+C*P[j]^r
#it also its inverse
Elem:=proc(P,C,i,j,r) local k:

if not (i<>j and 1<=i and i<=nops(P) and 1<=j and j<=nops(P)) then
 RETURN(FAIL):
fi:

[[op(1..i-1,P),P[i]+C*P[j]^r, op(i+1..nops(P),P)],
[op(1..i-1,P),P[i]-C*P[j]^r, op(i+1..nops(P),P)]]:
end:
####end C13.txt

###end old stuff

#ExtractC(P,x): inputs a polynomial P in the list of variables x, outputs
#the set of coefficients
ExtractC:=proc(P,x) local n,i,x1:

n:=nops(x):

if n=1 then
 RETURN({seq(coeff(P,x[1],i),i=0..degree(P,x[1]))}):
fi:

x1:=x[2..n]:

{seq( op(ExtractC(coeff(P,x[1],i),x1))   ,i=0..degree(P,x[1]))}:
end:

#Vecs1(d,r): all vectors of non-neg. integers of length r that  add-up to d
Vecs1:=proc(d,r) local S,d1,S1,v:

option remember:

if r=0 then
  if d=0 then
    RETURN({[]}):
  else
    RETURN({}):
  fi:
fi:

S:={}:

for d1 from 0 to d do
  S1:=Vecs1(d-d1,r-1):

  S:=S union {seq( [op(v),d1], v in S1)}:
od:

S:

end:

#Vecs(d,r): all vectors of length r whose 0<sum <=d
Vecs:=proc(d,r) local d1:
option remember:
{seq(op(Vecs1(d1,r)),d1=1..d)}:
end:



#GenPol(x,d,a): inputs a list of variables x, a degree d, and a symbol a
#outputs a generic polynomial of degree d with coeffices of the form a[.]
#followed by the SET of coefficients
GenPol:=proc(x,d,a) local S,P,coe,v,i:
S:=Vecs(d, nops(x)):

  add(a[op(v)]*mul(x[i]^v[i],i=1..nops(x)), v in S),  
 {seq(a[op(v)], v in S)}  :

end: