#OK to Post

#Mingjia Yang 3/12/2016

#HW 14


###########
#Problem 1
###########

#RandElem(x,MaxC,Maxr) inputs a list of variables x, a large positive integer MaxC, and a
#not-too-large (in applications) positive integer Maxr, and first picks random 1 ≤ i ≠ j ≤ nops(x), 
#then a random r between 1 and Maxr, then a random C between -MaxC and MaxC, and outputs the 
#pair of transformations [Elem(x,C,i,j,r), Elem(x,-C,i,j,r)]

RandElem:=proc(x,MaxC,Maxr) local i,j,r,C:

i:=rand(1..nops(x))():
j:=(rand(1..(nops(x)-1))()+i) mod nops(x) :
if j=0 then 
j:=nops(x):
fi:
r:=rand(Maxr)():
C:=rand(MaxC)():
[Elem(x,C,i,j,r),Elem(x,-C,i,j,r)]:
end:

###########
#Problem 2
###########

#RandEpair(x,MaxC,Maxr,T) uses RandElem(x,MaxC,Maxr) T times and outputs a pair of transformations,
#inverses of each other. 

RandEpair:=proc(x,MaxC,Maxr,T) local P,i,a,L,V:
P:=[x,x]:
L:=[seq(RandElem(x,MaxC,Maxr),i=1..T)]:


for i from 1 to T do
P[1]:=Comp(P[1],L[i][1],x):
P[2]:=Comp(P[2],L[T-i+1][2],x):
od:
print(P):
end:

#It's verified that P[1] and P[2] are inverse to each other.


###########
#Problem 3
###########

#Invk(P[1],x,N) seemed to give different inverses than P[2]...

###########
#Problem 4
###########

#Vecs(n,k)  inputs positive integers n and k and outputs the set of 0-1 vectors (given as lists) 
#of length n that add-up to k

Vecs:=proc(n,k) local L1,L2:

if k=0 then 
RETURN({[seq(0,i=1..n)]}):
fi:

if n=k then 
RETURN({[seq(1,i=1..n)]}):
fi:

if n<k then
RETURN({}):
fi:

L1:=convert(Vecs(n-1,k-1),list):
L2:=convert(Vecs(n-1,k),list):

{seq([1,op(L1[i])],i=1..nops(L1))} union {seq([0,op(L2[i])],i=1..nops(L2))}:


end:

#ZeroOneMatrices(RowSums,ColumnSums) inputs lists of non-negative integers 
#RowSums, ColumnSums, of lengtha m, and n, say (i.e. nops(RowSums)=m, 
#nops(ColumnSums)=n) and such that they have the same sum (k, say, if not it 
#returns FAIL), and outputs the set of m by n 0-1 matrices (given as lists-of-lists) 
#whose row sums are given by the list RowSums and whose column-sums are 
#given by the list ColumnSums.

ZeroOneMatrices1:=proc(RowSums,ColumnSums) local m,n,k,L,U,a,i,j,r,U1,b:

m:=nops(RowSums):
n:=nops(ColumnSums):
k:=sum(RowSums[i],i=1..m):

if m=1 then 
RETURN ({[[seq(ColumnSums[i],i=1..n)]]}):
fi:


L:=convert(Vecs(n,RowSums[m]),list):  
U:={}:

for i from 1 to nops(L) do

if i=nops(L)+1 then
RETURN(U):
fi:

a:=convert(ZeroOneMatrices1(RowSums[1..m-1],[seq(ColumnSums[j]-L[i][j],j=1..n)]),list): 
U:=U union {seq([op(a[j]),L[i]],j=1..nops(a))}:
od:
U:=convert(U,list):
end:


ZeroOneMatrices:=proc(RowSums,ColumnSums) local k,m,n,a1:
a1:={1}:
m:=nops(RowSums):
n:=nops(ColumnSums):
k:=sum(RowSums[i],i=1..m):

if k<>sum(ColumnSums[i],i=1..n) then
RETURN({}):
fi:

Filter(ZeroOneMatrices1(RowSums,ColumnSums)):


end:



Filter:=proc(L) local a,n,i,j,b,L1:
a:={}:
n:=nops(L[1][1]):
for i from 1 to nops(L) do
for j from 1 to n do 
if L[i][1][j]<0 or L[i][1][j]>1 then 
a:=a union{i}:
fi:
od:
od:
b:={seq(i,i=1..nops(L))} minus a:
L1:={}:
for i from 1 to nops(L) do
if member(i,b) then
L1:=L1 union {L[i]}:
fi:
od:
L1:

end:



#ZeroOneMatrices([1,1,1,1],[1,1,1,1]) is 

#{[[0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]], [[0, 0, 0, 1], [0, 0, 1, 0], [1, 0, 0, 0], [0, 1, 0, 0]], [[0, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [1, 0, 0, 0]], [[0, 0, 0, 1], [0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0]], [[0, 0, 0, 1], [1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0]], [[0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]], [[0, 0, 1, 0], [0, 0, 0, 1], [0, 1, 0, 0], [1, 0, 0, 0]], [[0, 0, 1, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0]], [[0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1], [1, 0, 0, 0]], [[0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1]], [[0, 0, 1, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 1, 0, 0]], [[0, 0, 1, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1]], [[0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [1, 0, 0, 0]], [[0, 1, 0, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 0, 1, 0]], [[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [1, 0, 0, 0]], [[0, 1, 0, 0], [0, 0, 1, 0], [1, 0, 0, 0], [0, 0, 0, 1]], [[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]], [[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]], [[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]], [[1, 0, 0, 0], [0, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0]], [[1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [0, 1, 0, 0]], [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]], [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]], [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]}



#seq(nops(ZeroOneMatrices([1$i],[1$i])),i=1..6) returns 1, 2, 6, 24, 120, 720;  It's in the OEIS.

#seq(nops(ZeroOneMatrices([2$i],[2$i])),i=1..5) returns 0, 1, 6, 90, 2040; It's in the OEIS.

#seq(nops(ZeroOneMatrices([3$i],[3$i])),i=1..5) returns 0,0,1,24,2040; It's in the OEIS.

#seq(nops(ZeroOneMatrices([2$i],[i,i])),i=1..5) returns 1,1,1,1,1; It's in the OEIS.

#seq(nops(ZeroOneMatrices([3$i],[i,i,i])),i=1..5) returns 1,1,1,1,1; It's in the OEIS.




##########################################
#Old Stuff


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

if r=0 and C=0 then
RETURN(FAIL):
fi:

[op(1..i-1,P),P[i]+C*P[j]^r, op(i+1..nops(P),P)]:
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:

####
#


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