#$Id: specific.txt,v 1.2 1996/10/07 14:54:17 tscoffe2 Exp tscoffe2 $ # # homoMap # # Constructs all monomials of degree g: z^a w^b Z^c W^d, a+b+c+d = g # input degree g # output set of coefficient lists [ [1,,,,0]], [1,,,,0], ] homoMap := proc(g) local a,b,c; {seq(seq(seq([1,a,b,c,g-a-b-c,0],c=0..g-a-b),b=0..g-a),a=0..g)}; end: # # h # # Applies group action to a specific coefficient list. # input: [1,,,,0], p and q picked such that (p,q)=1 and q!=+-1(mod p) # typical examples are 2,1 5,2 5,3 Good one: 8,3 # output: a number, if an integer then monomial is invariant under # group action. h := proc(l,p,q) (-l[2]-l[3]*q+l[4]+l[5]*q)/p; end: # # Invmonos # # Applies the group action:, h, to each item in set and removes those that # are not invariant under the action. # input: { [,,,], [,,,], }, p, q takes set input. # output: { [,,,], [,,,], } Invmonos := proc(L,p,q) local l,x; l := NULL; for x in L do if type(h(x,p,q),integer) then l := l,x; fi; od; {l}; end: # # listdiff # # This will take the derivative of every elment in the list with respect # to the list element specified. # input: [ [e,a,b,c,d,u], ], x an integer # output: [ [e,a,b,c,d,u], ] listdiff := proc(L,x) local i,l; l := NULL; for i in L do if i[x] > 0 then l := l,[i[1]*i[x],op(2..x-1,i),i[x]-1,op(x+1..nops(i),i)]; fi; od; [l]; end: # # listprod # # This will take two polynomials in list form and take their product. # no collection is done for I don't know what I will need to collect # on. # input: [ [,,,,,], [,,,,,], ], [ [,,,,,], [,,,,,], ] # output: [ [,,,,,], [,,,,,], ] # This multiplies the first element, adds the next four, and puts the # last element in a list. listprod := proc(L1,L2) local i,j; [seq(seq([ L1[i][1]*L2[j][1], L1[i][2]+L2[j][2], L1[i][3]+L2[j][3], L1[i][4]+L2[j][4], L1[i][5]+L2[j][5], sort([L1[i][6],L2[j][6]]) ],j=1..nops(L2)),i=1..nops(L1))]; end: # # Rightside # # Outputs the right hand side of the equation in list form. # input: equation number n, degree g # output: [ [,,,], [,,,], ] Rightside := proc(n,g) local k; if n = 1 then [seq([1,g-2-k,k,g-k,k,0],k=0..g-2)]; elif n = 2 then [seq([1,g-2-k,k,g-2-k,k+2,0],k=0..g-2)]; elif n = 3 then [seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)]; elif n = 4 then [seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)]; elif n = 5 then [seq([1,g-2-k,k,g-1-k,k+1,0],k=0..g-2)]; elif n = 6 then [seq([1,g-k-1,k,g-2-k,k+1,0],k=0..g-2)]; # New stuff here: elif n = 7 then [seq([1,g-k,k,g-k-2,k,0],k=0..g-2)]; elif n = 8 then [seq([1,g-k-2,k,g-k,k,0],k=0..g-2)]; elif n = 9 then [seq([1,g-1-k,k+1,g-2-k,k,0],k=0..g-2)]; elif n = 10 then [seq([1,g-2-k,k+1,g-1-k,k,0],k=0..g-2)] fi; end: # # CoeffGrab # # Note: the input to this should be collected on the monomials first. # # This procedure will take in two polynomials and a choice monomial # Each monomial in the first polynomial is then taken and a corresponding # monomial is constructed that would multiply to the choice monomial. This # monomial is then checked for in the second polynomial. This is continued # until one is found in the second polynomial. It is then multiplied # and added to a queue. This queue is then output as a polynomial # # This will now also take in the Rightside monomials and subtract them # at the end from the current coefficients, if necessary. CoeffGrab := proc(A,B,M,R) local m,n,N,G,L,queue,this,j; L := NULL; for m in A do # printf(`m = %a\n`,m); N := [1,M[2]-m[2],M[3]-m[3],M[4]-m[4],M[5]-m[5],0]; G := NULL; for this in B do if [op(2..5,this)] = [op(2..5,N)] then G := G,this; fi; od; G := [G]; queue := NULL; for n in G do queue := queue,[n[1]*m[1],sort([n[6],m[6]])]; od; if queue <> NULL then L := L,queue; fi; od; for j from 1 to nops(R) do if op(2..5,M) = op(2..5,R[j]) then L := L,[-R[j][1],[0,0]]; fi; od; [L]; end: # # RightPoly # # This constructs the righthand side in list form. With coefficients. # It also changes the sign on the constant term so that they # can just be appended to the list of monomials. # # fzfz fwfw fzfZ fwfW fzfw fZfw fZfZ fWfW fZfW fzfW # RightPoly := proc(n,g) local k; if n = 1 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-k-2)!*k!),g-k-2,k,g-k,k,0],k=0..g-2)]; elif n = 2 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-k-2)!*k!),g-k-2,k,g-k-2,k+2,0],k=0..g-2)]; elif n = 3 then [[1/2,0,g-1,0,g-1,0],seq([((g-2)!/((g-2-k)!*k!))*((g-1)/(2*(g-1-k))+3*g/(4*g+8)-1/4),g-1-k,k,g-1-k,k,0],k=0..g-2)]; elif n = 4 then [seq([1/2*(g-2)!/((g-2-k)!*k!),g-1-k,k,g-1-k,k,0],k=0..g-2), seq([(3*g/(4*g+8)+1/4)*(g-2)!/((g-2-k)!*k!),g-2-k,k+1,g-2-k,k+1,0],k=0..g-2)]; elif n = 5 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-2-k)!*k!),g-2-k,k,g-1-k,k+1,0],k=0..g-2)]; elif n = 6 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-2-k)!*k!),g-1-k,k,g-2-k,k+1,0],k=0..g-2)]; #New stuff here: elif n = 7 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-k-2)!*k!),g-k,k,g-k-2,k,0],k=0..g-2)]; elif n = 8 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-k-2)!*k!),g-k-2,k,g-k,k,0],k=0..g-2)]; elif n = 9 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-2-k)!*k!),g-1-k,k+1,g-2-k,k,0],k=0..g-2)]; elif n = 10 then [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-2-k)!*k!),g-2-k,k+1,g-1-k,k,0],k=0..g-2)] fi; end: # # RightSide # # This proc outpus the polynomials on the right hand side. Without # coefficients. # input: equation number n, degree g # output: [ [,,,], [,,,], ] #RightSide := proc(n,g) #local k; #if n = 1 then #[seq([1,g-2-k,k,g-k,k,0],k=0..g-2)]; #elif n = 2 then #[seq([1,g-2-k,k,g-2-k,k+2,0],k=0..g-2)]; #elif n = 3 then #[seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)]; #elif n = 4 then #[seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)]; #elif n = 5 then #[seq([1,g-2-k,k,g-1-k,k+1,0],k=0..g-2)]; #elif n = 6 then #[seq([1,g-k-1,k,g-2-k,k+1,0],k=0..g-2)]; # ## New stuff here: #elif n = 7 then # [seq([1,g-k,k,g-k-2,k,0],k=0..g-2)]; #elif n = 8 then # [seq([1,g-k-2,k,g-k,k,0],k=0..g-2)]; #elif n = 9 then # [seq([1,g-1-k,k+1,g-2-k,k,0],k=0..g-2)]; #elif n = 10 then # [seq([1,g-2-k,k+1,g-1-k,k,0],k=0..g-2)] #fi; #end: # # CollectMonomials # # This proc should take in a polynomial in list form and collect on all # the monomials by adding together the constant coefficients. This will # only occur if a,b,c,d, and the u term are the same. # #CollectMonomials := proc(A); #local m,B; # B := NULL; # while nops(A) > 0 do # m := A[1]; # A := [op(2..nops(A),A)]; # for i from 1 to nops(A) do # if [op(2..6,A[i])] = [op(2..6,m)] then # m := [m[1]+A[i][1],op(2..6,m)]; # # fi; # od; # od; # #end: varmax := 1000: # ********************************************************** # # basic # basic := proc(g,p,q) local homo,n,harmonics,nakedMono,polyhold,a,b,c,d,capN,i,coeff1,coeff2, mono1,mono2,f,fz,fw,fZ,fW,arbitrary,fzfz,fzfZ,fwfw,fwfW,fzfw,fwfZ, cl,ch,RHS,st,T,r,nakedMono_conj,Check,prodcts,j; global varmax; homo := homoMap(g); n := Invmonos(homo,p,q); print(`Constructing harmonics`); harmonics := NULL; arbitrary := 0; #varmax := 1000; while nops(n) > 0 do #print(nops(n)); arbitrary := arbitrary+1; nakedMono := n[1]; #print(`NakedMono`,nakedMono); n := n minus { nakedMono }; if member(ListConjugate(nakedMono), n, 'r') then nakedMono_conj := n[r]; n := n minus { nakedMono_conj }; else nakedMono_conj := NULL; fi; # nakedMono here: polyhold := [op(1..5,nakedMono),arbitrary]; a := nakedMono[2]; b := nakedMono[3]; c := nakedMono[4]; d := nakedMono[5]; capN := max(min(a,c),min(b,d)); for i from 0 to capN do #print(capN,i); coeff1:=(-1)^(i+1)*product((a-k)*(c-k)/((b+k+1)*(d+k+1)),k=0..i); #print(`Coeff1`,coeff1); mono1:=[coeff1,a-i-1,b+i+1,c-i-1,d+i+1,arbitrary]; coeff2:=(-1)^(i+1)*product((b-k)*(d-k)/((a+k+1)*(c+k+1)),k=0..i); #print(`Coeff2`,coeff2); mono2:=[coeff2,a+i+1,b-i-1,c+i+1,d-i-1,arbitrary]; #print(`Mono1: `,mono1,` Mono2: `,mono2); if coeff1 <> 0 then polyhold := polyhold,mono1; fi; if coeff2 <> 0 then polyhold := polyhold,mono2; fi; n := n minus { [1,a-i-1,b+i+1,c-i-1,d+i+1,0], [1,a+i+1,b-i-1,c+i+1,d-i-1,0] }; od; #print(`Polyhold`,polyhold); harmonics := harmonics,polyhold; # nakedMono_conj here: if nakedMono_conj <> NULL then nakedMono := nakedMono_conj; polyhold := [op(1..5,nakedMono),evalf(arbitrary/varmax)]; a := nakedMono[2]; b := nakedMono[3]; c := nakedMono[4]; d := nakedMono[5]; capN := max(min(a,c),min(b,d)); for i from 0 to capN do #print(capN,i); coeff1:=(-1)^(i+1)*product((a-k)*(c-k)/((b+k+1)*(d+k+1)),k=0..i); #print(`Coeff1`,coeff1); mono1:=[coeff1,a-i-1,b+i+1,c-i-1,d+i+1,evalf(arbitrary/varmax)]; coeff2:=(-1)^(i+1)*product((b-k)*(d-k)/((a+k+1)*(c+k+1)),k=0..i); #print(`Coeff2`,coeff2); mono2:=[coeff2,a+i+1,b-i-1,c+i+1,d-i-1,evalf(arbitrary/varmax)]; #print(`Mono1: `,mono1,` Mono2: `,mono2); if coeff1 <> 0 then polyhold := polyhold,mono1; fi; if coeff2 <> 0 then polyhold := polyhold,mono2; fi; n := n minus { [1,a-i-1,b+i+1,c-i-1,d+i+1,0], [1,a+i+1,b-i-1,c+i+1,d-i-1 ,0] }; od; #print(`Polyhold`,polyhold); harmonics := harmonics,polyhold; fi; od; f := [harmonics]; print(`Starting differentiation`); fz := listdiff(f,2); fw := listdiff(f,3); fZ := listdiff(f,4); fW := listdiff(f,5); save(f,fz,fw,fZ,fW,`data_diff.m`); print(`Starting product and colleciton.`); cl := NULL; ch := NULL; prodcts := [[fz,fz],[fw,fw],[fz,fZ],[fw,fW],[fz,fw],[fw,fZ],[fZ,fZ],[fW,fW],[fZ,fW],[fz,fW]]; for j from 1 to nops(prodcts) do print(cat(`Polynomial #`,convert(j,string))); st := time(); cl := cl,SetUpEqn(op(prodcts[j]),g,j); print(time()-st); T := `Poly `.convert(j,string).` done!`; save(T,`cl.txt`); od; print(`Finished constructing equations.`); {cl}; end: # # This is the kernel in basic for the coefficient comparison. # SetUpEqn := proc(u,v,g,n) local RHS,Right,cl,ch,i,Check,Check2; RHS := Rightside(n,g); Right := RightPoly(n,g); cl := seq(CoeffGrab(u,v,RHS[i],Right),i=1..nops(RHS)); Check := [op({op(map(x->[x[2..5]],listprod(u,v)))} minus {op(map(x->[x[2..5]],Right))})]; ch := seq(CoeffGrab(u,v,[0,op(Check[i]),0],[[0,0,0,0,0,0]]),i=1..nops(Check)); cl,ch; end: ListConjugate := proc(L) [L[1],L[4],L[5],L[2],L[3],L[6]]; end: # # Note this g only works because N = g+1 for g = 2 and 6. # # I have written this code to accept the list format output of basic and to output # the polynomial format that will be required by solve. # TDerivativeAtZero := proc(L,g) local output,intermed,u,v,k,c,d,N; global varmax; if g = 2 then N := 8 elif g = 6 then N := 6 fi; output := NULL; for u in L do intermed := NULL; for v in u do if v[2][1] <> 0 then if ( type(v[2][2],float) and type(v[2][1],float) ) then c := round(v[2][1]*varmax); d := round(v[2][2]*varmax); intermed := intermed, convert([seq( v[1]*(cat(A,_,convert(k,string),_,convert(d,string),_,'RE') -I*cat(A,_,convert(k,string),_,convert(d,string),_,'IM')) *conjugate(a(k,c)) +v[1]*(cat(A,_,convert(k,string),_,convert(c,string),_,'RE') -I*cat(A,_,convert(k,string),_,convert(c,string),_,'IM')) *conjugate(a(k,d)),k=1..N+1)],`+`); elif type(v[2][1],float) then c := round(v[2][1]*varmax); d := v[2][2]; intermed := intermed, convert([seq( v[1]*(cat(A,_,convert(k,string),_,convert(d,string),_,'RE') +I*cat(A,_,convert(k,string),_,convert(d,string),_,'IM')) *conjugate(a(k,c)) +v[1]*(cat(A,_,convert(k,string),_,convert(c,string),_,'RE') -I*cat(A,_,convert(k,string),_,convert(c,string),_,'IM')) *a(k,d),k=1..N+1)],`+`); elif type(v[2][2],float) then c := v[2][1]; d := round(v[2][2]*varmax); intermed := intermed, convert([seq( v[1]*(cat(A,_,convert(k,string),_,convert(d,string),_,'RE') -I*cat(A,_,convert(k,string),_,convert(d,string),_,'IM')) *a(k,c) +v[1]*(cat(A,_,convert(k,string),_,convert(c,string),_,'RE') +I*cat(A,_,convert(k,string),_,convert(c,string),_,'IM')) *conjugate(a(k,d)),k=1..N+1)],`+`); else c := v[2][1]; d := v[2][2]; intermed := intermed, convert([seq( v[1]*(cat(A,_,convert(k,string),_,convert(d,string),_,'RE') +I*cat(A,_,convert(k,string),_,convert(d,string),_,'IM')) *a(k,c) +v[1]*(cat(A,_,convert(k,string),_,convert(c,string),_,'RE') +I*cat(A,_,convert(k,string),_,convert(c,string),_,'IM')) *a(k,d),k=1..N+1)],`+`); fi; fi; od; output := output,convert([intermed],`+`); od; {output}; end: EndListToPoly := proc(L) local i,j; subs(u_0_0=1,{seq(convert([seq(L[i][j][1]*cat(u,_,convert(L[i][j][2][1],string),_,convert(L[i][j][2][2],string)),j=1..nops(L[i]))],`+`),i=1..nops(L))}); end: BasicListToPoly := proc(L) local i,j,d1,d2,M,N; N := NULL; for i in L do M := NULL; for j in i do if type(j[2][1],float) then d1 := cat(c,convert(round(j[2][1]*1000),string)); else d1 := convert(j[2][1],string); fi; if type(j[2][2],float) then d2 := cat(c,convert(round(j[2][2]*1000),string)); else d2 := convert(j[2][2],string); fi; # M := M,j[1]*``.u.`_`.d1.`_`.d2; M := M,j[1]*cat(u,_,d1,_,d2); od; N := N,convert({M},`+`); od; subs(u_0_0=1,{N}); end: PolyListToPoly := proc(L) local i,j; convert([seq(L[i][1]*z^L[i][2]*w^L[i][3]*zbar^L[i][4]*wbar^L[i][5]*u(L[i][6]),i=1..nops(L))],`+`); end: PolyToList := proc(solu) local Ii,Jj,i,j,d1,d2; Ii := NULL; for i in solu do Jj := NULL; if type(i,`function`) then d1 := op(1,i); d2 := op(2,i); if type(op(1,i),`function`) then d1 := evalf(op(op(1,i))/1000); fi; if type(op(2,i),`function`) then d2 := evalf(op(op(2,i))/1000); fi; Jj := Jj,[1,[d1,d2]]; else for j in i do if has(j,u) then if type(j,`*`) then d1 := op(1,op(2,j)); d2 := op(2,op(2,j)); if type(op(1,op(2,j)),`function`) then d1 := evalf(op(op(1,op(2,j)))/1000); fi; if type(op(2,op(2,j)),`function`) then d2 := evalf(op(op(2,op(2,j)))/1000); fi; Jj := Jj,[op(1,j),[d1,d2]]; else d1 := op(1,j); d2 := op(2,j); if type(op(1,j),`function`) then d1 := evalf(op(op(1,j))/1000); fi; if type(op(2,j),`function`) then d2 := evalf(op(op(2,j))/1000); fi; Jj := Jj,[1,[d1,d2]]; fi; else Jj := Jj,[j,[0,0]]; fi; od; fi; Ii := Ii,[Jj]; od; [Ii]; end: FPolyToList := proc(Ldata) local A,B,C,bs,i,j,x,split,f; B := [0,`z`,`w`,`zbar`,`wbar`,`u`]; f := NULL; for A in Ldata do split := {op(A*z^2*w^2*zbar^2*wbar^2)}; C := 'C'; for i from 2 to 5 do bs := select(has,split,B[i]); if nops(bs) = 1 then C(i) := op(2,bs[1])-2; split := split minus bs; else C(i) := 0; fi; od; bs := select(has,split,B[6]); if nops(bs) = 1 then if type(op(bs[1]),function) then C(6) := evalf(op(op(bs[1]))/1000); else C(6) := op(bs[1]); fi; split := split minus bs; else C(6) := 0; fi; if nops(split) = 1 then C(1) := op(split); else C(1) := 1; fi; f := f,[seq(C(j),j=1..6)]; od; [f]; end: