########## # File @(#) SfAExpn (SFA PACKAGE) Mercredi 23 juillet 1997 15:51:21 ########## interface(verboseproc=0); # # SfAExpand(sfa) converts subterms in sfa that are on polynomial alphabet # expressions into terms on one single alphabet. # # The following formulae are used (some of them should be optimized later): # # p[I](-A) = (-1)^size(I) p[I](A) # e[I](-A) = (-1)^|I| h[I](A) # h[I](-A) = (-1)^|I| e[I](A) # s[I](-A) = (-1)^|I| s[I~](A) # m[I](-A) = TopA(m[I](A)) # # p[I](AB) = p[I](A) p[I](B) # e[I](AB) = (-1)^|I| product{i=1..size(I)} # sum{|J|=I.i} (-1)^|J| ToeA(m[J](A)) e[J](B) # h[I](AB) = product{i=1..size(I)} # sum{|J|=I.i} TohA(m[J](A)) h[J](B) # s[I](AB) = sum{|J|=|I|} s[I*J](A) s[J](B) where * = SfInternal # m[I](AB) = TopA(m[I](AB)) # # For any k <> 0, 1 : # p[I](kA) = k^(size(I)) p[I](A) # e[I](kA) = (-1)^|I| product{i=1..size(I)} # sum{|J|=I.i} (-1)^|J| m[J](k) h[J](A) # h[I](kA) = product{i=1..size(I)} # sum{|J|=I.i} m[J](k) h[J](A) # s[I](kA) = sum{|J|=|I|} s[I*J](k) s[J](A) where * = SfInternal # m[I](kA) = TopA(m[I](kA)) # # p[I](A+B) = product{i=1..size(I)} p[I.i](A) + p[I.i](B) # e[I](A+B) = product{i=1..size(I)} # sum{j=0..I.i} e[I.i-j](A) e[j](B) # h[I](A+B) = product{i=1..size(I)} # sum{j=0..I.i} h[I.i-j](A) h[j](B) # s[I](A+B) = sum{J in I} s[J](B) SfDiff(s[J], s[I])(A) # m[I](A+B) = sum{I' U I"=I} m[I'](A) m[I"](B) # # For any k <> 0, 1 : # p[I](k) = k^size(I) # e[I](k) = product{i=1..size} # k(k-1)(k-2)..(k-I.i+1)/I.i! # h[I](k) = product{i=1..size} # k(k+1)(k+2)..(k+I.i-1)/I.i! # m[I](k) = product{i=1..size} # k(k-1)(k-2)..(k-n+1)/(a1!a2!..an!) with I=1^a1..n^an # s[I](k) = product{i=1..size} ld[i]*(i-dp+k) / # product(hooks(I)) with (dp, ld) = Part2Diagonal(I) # # VP230797 VP120997 VP151097 VP120298 VP020398 VP030398 VP100299 VP110299 `SFA/SfAExpand` := proc(sfa) local part, # indexing partition... alphab, # alphabet... base, # a basis... lp, # a list of partitions... l, # a list of lists... i, # var for loop and seq... j, # var for loop... k, # var for seq... v, # var... ra, # remaining alphabet... t; # a table... global SFABases; # known bases... option remember; if type(sfa, `+`) then map(`SFA/SfAExpand`, sfa); elif type(sfa, `*`) then t := `SFA/MonomialA`(sfa, SFABases); for i in map(op, [indices(t)]) do if (i <> 'coeff') then for j in map(op, [indices(t[i])]) do if type(t[i][j], `*`) then # VP020398 t['coeff'] := t['coeff'] * eval(cat(`SFA/To`,i,`A`)(expand(map(`SFA/SfAExpand`, t[i][j](j))))); else t['coeff'] := t['coeff'] * (`SFA/SfAExpand`(t[i][j](j))); fi; od; fi; od; expand(t['coeff']); elif type(sfa, `^`) then base := op(0, op(0, op(1, sfa))); if (base='p' or base ='e' or base='h') then # To_A is much faster than SfACollect... `SFA/SfAExpand`(eval(cat(`SFA/To`, base, `A`)(sfa))); elif (type(op(2, sfa), posint) and (base='s' or base='m')) then `SFA/SfAExpand`(eval(cat(`SYMF/To`, base)(op(0, op(1, sfa))^op(2,sfa))) (op(1,op(1, sfa)))); else sfa; fi; elif (type(sfa, function) and (`TYP/IspA`(sfa) or `TYP/IseA`(sfa) or `TYP/IshA`(sfa) or `TYP/IssA`(sfa) or `TYP/IsmA`(sfa))) then part := [op(op(0, sfa))]; base := op(0, op(0, sfa)); alphab := op(1, sfa); # First, we check the partition... if (part=[]) then if (base='p') then # can not handle such a case... ERROR(`p[](...) can not be expanded!`); else 1; fi; elif type(alphab, `^`) then # Alphabet is a power... if type(op(2, alphab), posint) then if (base='p') then `SFA/SfAExpand`(p[op(part)](op(1, alphab))) * `SFA/SfAExpand`(p[op(part)](op(1, alphab)^(op(2, alphab)-1))); elif (base='e') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); v := v, convert([ seq((-1)^(convert(lp[j], `+`))* `SFA/SfAExpand`(`SFA/ToeA`( m[op(lp[j])](op(1, alphab)))) * `SFA/SfAExpand`(e[op(lp[j])](op(1, alphab)^(op(2, alphab)-1))), j = 1..nops(lp))], `+`); od; (-1)^(convert(part, `+`))*`SFA/ToeA`(expand(convert([v], `*`))); elif (base='h') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); v := v, convert([ seq( `SFA/SfAExpand`(`SFA/TohA`(m[op(lp[j])](op(1, alphab)))) * `SFA/SfAExpand`(h[op(lp[j])](op(1, alphab)^(op(2, alphab)-1))), j = 1..nops(lp))], `+`); od; `SFA/TohA`(expand(convert([v], `*`))); elif (base='s') then lp := `PART/ListPart`(convert(part, `+`)); `SFA/TosA`(convert([ # VP110299 seq(`SFA/SfAExpand`(expand( # VP110299 `SFA/Sf2SfA`( `SYMF/SfInternal`(s[op(part)], s[op(lp[i])], 's'),op(1, alphab)) * s[op(lp[i])](op(1, alphab)^(op(2, alphab)-1)))), i=1..nops(lp))], `+`)); elif (base='m') then `SFA/SfAExpand`(`SFA/TopA`(sfa)); fi; elif type(op(2,alphab),negint) or (type(op(2,alphab),`*`) and sign(op(2,alphab))=-1) then # VP110299 # This is a fraction: we put the result on the p-basis... if (base='p') then `SFA/SfAExpand`(p[op(part)](numer(alphab)))/ convert([seq(subs(map((x,y)->x=x^y, SfAVars(),part[i]), denom(alphab)), i=1..nops(part))],`*`); else `SFA/SfAExpand`(`SFA/TopA`(sfa, [alphab])); fi; elif `TYP/IsConstant`(alphab) then # The power is a constant... if (base=`p`) then alphab^nops(part); elif (base=`e`) then convert( [seq(convert([seq(alphab-j, j=0..part[i]-1)], `*`)/factorial(part[i]), i=1..nops(part))], `*`); elif (base=`h`) then convert( [seq(convert([seq(alphab+j, j=0..part[i]-1)], `*`)/factorial(part[i]), i=1..nops(part))], `*`); elif (base=`s`) then convert([seq(seq(alphab-i+j ,j=1..part[i]), i=1..nops(part))], `*`)/ convert(`PART/Part2ListHook`(part), `*`); elif (base=`m`) then v := `PART/Part2Exp`(part); convert([seq(alphab-i, i=0..nops(v)-1)], `*`) / convert(map(factorial, v), `*`); fi; else sfa; fi; elif type(alphab, `*`) then if type(op(1, alphab), negative) then # Negative alphabet... if (base='p') then (-1)^(nops(part)) * `SFA/SfAExpand`(p[op(part)](-alphab)); elif (base='e') then (-1)^(convert(part, `+`)) * `SFA/SfAExpand`(h[op(part)](-alphab)); elif (base='h') then (-1)^(convert(part, `+`)) * `SFA/SfAExpand`(e[op(part)](-alphab)); elif (base='s') then (-1)^(convert(part, `+`)) * `SFA/SfAExpand`(s[op(`PART/Part2Conjugate`(part))](-alphab)); elif (base='m') then `SFA/SfAExpand`(`SFA/TopA`(sfa)); fi; elif type(op(1, alphab), `^`) then # Product with a power as first term... ra := convert([op(2..nops(alphab), alphab)], `*`); if type(op(2, op(1, alphab)), posint) then if (base='p') then `SFA/SfAExpand`(p[op(part)](op(1, op(1, alphab)))) * `SFA/SfAExpand`(p[op(part)]( op(1, op(1, alphab))^(op(2, op(1, alphab))-1)*ra)); elif (base='e') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); v := v, convert([ seq((-1)^(convert(lp[j], `+`))* `SFA/SfAExpand`(`SFA/ToeA`( m[op(lp[j])](op(1, op(1, alphab))))) * `SFA/SfAExpand`(e[op(lp[j])]( op(1, op(1, alphab))^(op(2, op(1, alphab))-1)*ra)), j = 1..nops(lp))], `+`); od; (-1)^(convert(part, `+`))*`SFA/ToeA`(expand(convert([v], `*`))); elif (base='h') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); v := v, convert([ seq( `SFA/SfAExpand`(`SFA/TohA`(m[op(lp[j])]( op(1, op(1, alphab))))) * `SFA/SfAExpand`(h[op(lp[j])]( op(1, op(1, alphab))^(op(2, op(1, alphab))-1)*ra)), j = 1..nops(lp))], `+`); od; `SFA/TohA`(expand(convert([v], `*`))); elif (base='s') then lp := `PART/ListPart`(convert(part, `+`)); `SFA/TosA`(convert([ # VP110299 seq(`SFA/SfAExpand`(expand( # VP110299 `SFA/Sf2SfA`(`SYMF/SfInternal`(s[op(part)], s[op(lp[i])], 's'), op(1, op(1, alphab))) * s[op(lp[i])](op(1, op(1, alphab))^(op(2, op(1,alphab))-1)*ra))), i=1..nops(lp))], `+`)); elif (base='m') then `SFA/SfAExpand`(`SFA/TopA`(sfa)); fi; elif type(op(2,op(1,alphab)),negint) or (type(op(2,op(1,alphab)),`*`) and sign(op(2,op(1,alphab)))=-1) then # VP110299 # This is a fraction: we put the result on the p-basis... if (base='p') then `SFA/SfAExpand`(p[op(part)](numer(op(1,alphab))))/ convert([seq(subs(map((x,y)->x=x^y, SfAVars(),part[i]), denom(op(1,alphab))), i=1..nops(part))],`*`) * `SFA/SfAExpand`(p[op(part)](op(2..nops(alphab),alphab))); else `SFA/SfAExpand`(`SFA/TopA`(sfa, [alphab])); fi; else sfa; fi; elif `TYP/IsArgAlphabet/simple`(op(1, alphab)) or `TYP/IsVariable`(op(1,alphab)) then # Product of two alphabets... if (base='p') then `SFA/SfAExpand`(p[op(part)](op(1, alphab))) * `SFA/SfAExpand`(p[op(part)] (convert([op(2..nops(alphab), alphab)], `*`))); elif (base='e') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); v := v, convert([ seq((-1)^(convert(lp[j], `+`))* `SFA/SfAExpand`(`SFA/ToeA`( m[op(lp[j])](op(1, alphab)))) * `SFA/SfAExpand`(e[op(lp[j])](convert([op(2..nops(alphab), alphab)], `*`))), j = 1..nops(lp))], `+`); od; (-1)^(convert(part, `+`))*`SFA/ToeA`(expand(convert([v], `*`))); elif (base='h') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); v := v, convert([ seq( `SFA/SfAExpand`(`SFA/TohA`(m[op(lp[j])](op(1, alphab)))) * `SFA/SfAExpand`(h[op(lp[j])](convert([op(2..nops(alphab), alphab)], `*`))), j = 1..nops(lp))], `+`); od; `SFA/TohA`(expand(convert([v], `*`))); elif (base='s') then lp := `PART/ListPart`(convert(part, `+`)); `SFA/TosA`(convert([ # VP110299 seq(`SFA/SfAExpand`(expand( # VP110299 `SFA/Sf2SfA`(`SYMF/SfInternal`(s[op(part)], s[op(lp[i])], 's'), op(1, alphab)) * s[op(lp[i])](convert([op(2..nops(alphab), alphab)], `*`)))), i=1..nops(lp))], `+`)); elif (base='m') then `SFA/SfAExpand`(`SFA/TopA`(sfa)); fi; # VP110299 `TYP/IsConstant` matches any sum... I've added this just before. # 1st term of product is a sum => expand... elif type(op(1, alphab), `+`) then `SFA/SfAExpand`(base[op(part)](expand(alphab))); elif `TYP/IsConstant`(op(1, alphab)) then # Product of an alphabet by a numeric/formal constant if (base='p') then op(1, alphab)^(nops(part)) * `SFA/SfAExpand`(p[op(part)] (convert([op(2..nops(alphab), alphab)], `*`))); elif (base='e') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); l := [seq(`PART/Part2Exp`(lp[j]), j=1..nops(lp))]; v := v, convert([ seq((-1)^(part[i])* convert([seq(op(1, alphab)-k, k=0..convert(l[j], `+`)-1)], `*`) / convert(map(factorial, l[j]), `*`) * `SFA/SfAExpand`(e[op(lp[j])](convert([op(2..nops(alphab), alphab)], `*`))), j = 1..nops(lp))], `+`); od; (-1)^(convert(part, `+`))*`SFA/ToeA`(expand(convert([v], `*`))); elif (base='h') then v := NULL; for i from 1 to nops(part) do lp := `PART/ListPart`(part[i]); l := [seq(`PART/Part2Exp`(lp[j]), j=1..nops(lp))]; v := v, convert([seq( # of multinomials... convert([seq(op(1, alphab)-k, k=0..convert(l[j], `+`)-1)], `*`) / convert(map(factorial, l[j]), `*`) * `SFA/SfAExpand`(h[op(lp[j])] (convert([op(2..nops(alphab), alphab)], `*`))), j=1..nops(lp))], `+`); od; `SFA/TohA`(expand(convert([v], `*`))); elif (base='s') then v := convert(part, `+`); lp :=`PART/ListPart`(v); expand( convert([ seq(subs(seq(p.j=op(1, alphab), j=1..v), `SYMF/SfInternal`(s[op(part)], s[op(lp[i])], 's')) * `SFA/SfAExpand`(s[op(lp[i])] (convert([op(2..nops(alphab), alphab)], `*`))), i=1..nops(lp))], `+`)); elif (base='m') then # Nothing better for the moment... `SFA/SfAExpand`(`SFA/TopA`(sfa)); fi; elif (type(op(1,alphab), function) and (`TYP/IspA`(op(1,alphab)) or `TYP/IseA`(op(1,alphab)) or `TYP/IshA`(op(1,alphab)) or `TYP/IssA`(op(1,alphab)) or `TYP/IsmA`(op(1,alphab)))) then # Alphabet may be a product of symmetric function... if (`TYP/IspA`(op(1,alphab)) and base='p' and `TYP/IsArgAlphabet/simple`(op(1, op(1,alphab)))) then p[op(sort([ seq(op(map((n,m) -> n*m, [op(op(0, op(1,alphab)))],part[i])), i=1..nops(part))], (x,y) -> evalb(x>y)))](op(1, op(1,alphab)))* `SFA/SfAExpand`(p[op(part)] (convert([op(2..nops(alphab), alphab)], `*`))); else `SFA/SfAExpand`(`SFA/TopA`(base[op(part)] (`SFA/TopA`(`SFA/SfAExpand`(alphab))))); fi; else # donna what to do with this... sfa; fi; elif type(alphab, `+`) then # Sum of two alphabets... if (base='p') then `SFA/SfACollect`( expand(convert([ seq( `SFA/SfAExpand`(p[part[i]](op(1, alphab))) + `SFA/SfAExpand`(p[part[i]](convert([op(2..nops(alphab), alphab)], `+`))), i=1..nops(part)) ], `*`)), 'p'); elif (base='e') then `SFA/ToeA`(expand(convert([ seq(convert([ seq( `SFA/SfAExpand`(e[subs(0=NULL,part[i]-j)](op(1, alphab))) * `SFA/SfAExpand`(e[subs(0=NULL,j)] (convert([op(2..nops(alphab), alphab)], `+`))), j=0..part[i]) ], `+`), i=1..nops(part)) ], `*`))); elif (base='h') then `SFA/TohA`(expand(convert([ seq(convert([ seq( `SFA/SfAExpand`(h[subs(0=NULL,part[i]-j)](op(1, alphab))) * `SFA/SfAExpand`(h[subs(0=NULL,j)] (convert([op(2..nops(alphab), alphab)], `+`))), j=0..part[i]) ], `+`), i=1..nops(part)) ], `*`))); elif (base='s') then lp := `PART/ListPartIn`(part); expand( convert([ seq(`SFA/TosA`(`SFA/SfAExpand`( # VP100299 s[op(lp[i])](convert([op(1..nops(alphab)-1, alphab)], `+`))) * `SFA/SfAExpand`((`SYMF/SfDiff`(s[op(lp[i])], s[op(part)])) (op(nops(alphab), alphab)))), i=1..nops(lp))], `+`)); elif (base='m') then lp := `PART/SplitPart`(part, 2); convert([ seq(`SFA/TomA`(expand(`SFA/SfAExpand`( # VP100299 m[op(lp[i][1])](op(1, alphab))) * `SFA/SfAExpand`(m[op(lp[i][2])] (convert([op(2..nops(alphab), alphab)], `+`))))), i=1..nops(lp))], `+`); fi; elif `TYP/IsArgAlphabet/simple`(alphab) then # Alphabet is a single alphabet Ai... sfa; elif `TYP/IsVariable`(alphab) then # Alphabet is a variable... VP120298... if (base=`p`) or (base=`h`) then alphab^(convert(part,`+`)); elif (base=`e`) then if part=[] then 1; elif part=[1] then alphab; else 0; fi; elif (base=`s`) or (base=`m`) then if part=[] then 1; elif nops(part)=1 then alphab^(part[1]); else 0; fi; fi; elif (type(alphab, function) and (`TYP/IspA`(alphab) or `TYP/IseA`(alphab) or `TYP/IshA`(alphab) or `TYP/IssA`(alphab) or `TYP/IsmA`(alphab))) then # Alphabet may be a symmetric function... if (`TYP/IspA`(alphab) and base='p' and `TYP/IsArgAlphabet/simple`(op(1, alphab))) then p[op(sort([ seq(op(map((n,m) -> n*m, [op(op(0, alphab))],part[i])), i=1..nops(part))], (x,y) -> evalb(x>y)))](op(1, alphab)); else `SFA/SfAExpand`(`SFA/TopA`(base[op(part)] (`SFA/TopA`(`SFA/SfAExpand`(alphab))))); fi; else # Alphabet is a numerical/formal constant... if (base=`p`) then alphab^nops(part); elif (base=`e`) then convert( [seq(convert([seq(alphab-j, j=0..part[i]-1)], `*`)/factorial(part[i]), i=1..nops(part))], `*`); elif (base=`h`) then convert( [seq(convert([seq(alphab+j, j=0..part[i]-1)], `*`)/factorial(part[i]), i=1..nops(part))], `*`); elif (base=`s`) then convert([seq(seq(alphab-i+j ,j=1..part[i]), i=1..nops(part))], `*`)/ convert(`PART/Part2ListHook`(part), `*`); elif (base=`m`) then v := `PART/Part2Exp`(part); convert([seq(alphab-i, i=0..nops(v)-1)], `*`) / convert(map(factorial, v), `*`); fi; fi; else sfa; fi; end; # # No special check is made, since sfa must be first decomposed # VP120997 `SFA/SfAExpand/check` := proc(sfa) end; # # VP120997 `SFA/SfAExpand/interface` := proc(sfa) `SFA/SfAExpand`(args); end; # # VP120997 `SFA_PACK/SfAExpand` := proc(sfa) `SFA/SfAExpand/check`(args); `SFA/SfAExpand/interface`(args); end; # save `SFA/SfAExpand`, `SFA/SfAExpand/interface`, `SFA/SfAExpand/check`, `SFA_PACK/SfAExpand`, ``.SFALib.`/SFA/src/SfAExpn.m`; if (not assigned(_NOQUIT)) then quit; fi;