/* Sparse multinomial representation This file implements sparse multivariate polynomial representation, along with a few algorithms. - MakeMultiNomial(expression) converts an expression to internal sparse multinomial format. NormalForm converts back. New proposal for simplify: 1) IsRational a in a+b or a*b etcetera: combine into one rational. 2) if root expression is rational, call Simplify(a)/Simplify(b) or even better: by dividing out the pp of the two (best is gcd). 3) for normal expression: make multinomial, and return normalized content*pp */ /* We need to include this one because we're extending NormalForm */ /*TODO remove this dependency! */ Use("univar.rep/code.ys"); MultiSimp(_expr) <-- [ Local(vars); vars:=MultiExpressionList(expr); MultiSimp2(MM(expr,vars)); ]; 10 # MultiSimp2(_a / _b) <-- [ Local(c1,c2,gcd,cmn); gcd:=MultiGcd(a,b); c1:=MultiContent(a); c2:=MultiContent(b); gcd:=Gcd(c1[2][1][2],c2[2][1][2]); c1[2][1][2] := c1[2][1][2]/gcd; c2[2][1][2] := c2[2][1][2]/gcd; cmn:=Min(c1[2][1][1],c2[2][1][1]); c1[2][1][1] := c1[2][1][1] - cmn; c2[2][1][1] := c2[2][1][1] - cmn; (NormalForm(c1)/NormalForm(c2)) *(NormalForm(MultiPrimitivePart(a))/NormalForm(MultiPrimitivePart(b))); ]; 20 # MultiSimp2(expr_IsMulti) <-- [ NormalForm(MultiContent(expr))*NormalForm(MultiPrimitivePart(expr)); ]; 30 # MultiSimp2(_expr) <-- expr; MultiExpressionList(_expr) <-- VarList(expr,"IsMultiExpression"); 10 # IsMultiExpression(_x + _y) <-- False; 10 # IsMultiExpression(_x - _y) <-- False; 10 # IsMultiExpression( - _y) <-- False; 10 # IsMultiExpression(_x * _y) <-- False; 10 # IsMultiExpression(_x / _y) <-- False; 10 # IsMultiExpression(_x ^ y_IsInteger) <-- False; 11 # IsMultiExpression(_x ^ _y)_(IsInteger(Simplify(y))) <-- False; 10 # IsMultiExpression(x_IsConstant) <-- False; 100 # IsMultiExpression(_x) <-- True; 10 # IsMulti(MultiNomial(vars_IsList,terms_IsList)) <-- True; 20 # IsMulti(_anything) <-- False; RuleBase("MultiNomial",{vars,terms}); 10 # NormalForm(x_IsMulti/y_IsMulti) <-- NormalForm(x)/NormalForm(y); 20 # NormalForm(MultiNomial(vars_IsList,terms_IsList)) <-- [ Local(result); result:=0; ForEach(t,terms) [ result := result + t[2] * Factorize(vars^(t[1])); ]; result; ]; MultiContent(MultiNomial(vars_IsList,terms_IsList)) <-- [ Local(least,list); Set(list,Transpose(terms)); Set(least, Min(list[1])); MultiNomial(vars,{{least,Gcd(list[2])}}); ]; MultiPrimitivePart(_multi) <-- [ Local(cont); cont:=MultiContent(multi); cont := MultiNomial(cont[1],{{-cont[2][1][1],1/cont[2][1][2]}}); MultiNomialMultiply(multi, cont); ]; MM(_expr) <-- MM(expr,MultiExpressionList(expr)); MM(_expr,_vars) <-- MakeMultiNomial(expr,vars); MakeMultiNomial(_expr) <-- MakeMultiNomial(ExpandBrackets(expr),MultiExpressionList(expr)); LocalSymbols(a,vars) [ 10 # MakeMultiNomial(a_IsConstant,vars_IsList) <-- MultiNomial(vars,{{FillList(0,Length(vars)),a}}); ]; LocalSymbols(a,vars,pow) [ 20 # MultiSingleFactor(_vars,_a,_pow) <-- [ //Echo({vars,a,pow}); Local(term); term:={FillList(0,Length(vars)),1}; term[1][Find(vars,a)] := pow; MultiNomial(vars,{term}); ]; ]; LocalSymbols(a,vars) [ 20 # MakeMultiNomial(_a,vars_IsList)_(Contains(vars,a)) <-- MultiSingleFactor(vars,a,1); ]; LocalSymbols(x,y,vars) [ 30 # MakeMultiNomial(_x + _y,vars_IsList) <-- MakeMultiNomial(x,vars) + MakeMultiNomial(y,vars); ]; LocalSymbols(x,y,vars) [ 30 # MakeMultiNomial(_x * _y,vars_IsList) <-- MakeMultiNomial(x,vars) * MakeMultiNomial(y,vars); ]; 10 # MultiRemoveGcd(x_IsMulti/y_IsMulti) <-- [ Local(gcd); Set(gcd,MultiGcd(x,y)); Set(x,MultiDivide(x,{gcd})[1][1]); Set(y,MultiDivide(y,{gcd})[1][1]); x/y; ]; 20 # MultiRemoveGcd(_x) <-- x; LocalSymbols(x,y,vars,gcd,a,a) [ 30 # MakeMultiNomial(_x / _y,vars_IsList) <-- [ MultiRemoveGcd(MakeMultiNomial(x,vars)/MakeMultiNomial(y,vars)); ]; ]; MultiNomial(_vars,_x) + MultiNomial(_vars,_y) <-- MultiNomialAdd(MultiNomial(vars,x), MultiNomial(vars,y)); MultiNomial(_vars,_x) * MultiNomial(_vars,_y) <-- MultiNomialMultiply(MultiNomial(vars,x), MultiNomial(vars,y)); MultiNomial(_vars,_x) - MultiNomial(_vars,_y) <-- MultiNomialAdd(MultiNomial(vars,x), MultiNomialNegate(MultiNomial(vars,y))); - MultiNomial(_vars,_y) <-- MultiNomialNegate(MultiNomial(vars,y)); - MultiNomial(_vars,_x) <-- MultiNomialNegate(MultiNomial(vars,x)); x_IsMulti + (y_IsMulti/z_IsMulti) <-- ((x*z+y)/z); (y_IsMulti/z_IsMulti) + x_IsMulti <-- ((x*z+y)/z); (y_IsMulti/z_IsMulti) + (x_IsMulti/w_IsMulti) <-- ((y*w+x*z)/(z*w)); (y_IsMulti/z_IsMulti) - (x_IsMulti/w_IsMulti) <-- ((y*w-x*z)/(z*w)); (y_IsMulti/z_IsMulti) * (x_IsMulti/w_IsMulti) <-- ((y*x)/(z*w)); (y_IsMulti/z_IsMulti) / (x_IsMulti/w_IsMulti) <-- ((y*w)/(z*x)); x_IsMulti - (y_IsMulti/z_IsMulti) <-- ((x*z-y)/z); (y_IsMulti/z_IsMulti) - x_IsMulti <-- ((y-x*z)/z); (a_IsMulti/(c_IsMulti/b_IsMulti)) <-- ((a*b)/c); ((a_IsMulti/c_IsMulti)/b_IsMulti) <-- (a/(b*c)); ((a_IsMulti/b_IsMulti) * c_IsMulti) <-- ((a*c)/b); (a_IsMulti * (c_IsMulti/b_IsMulti)) <-- ((a*c)/b); - ((a_IsMulti)/(b_IsMulti)) <-- (-a)/b; LocalSymbols(x,vars) [ 30 # MakeMultiNomial(- _x,vars_IsList) <-- [ -MakeMultiNomial(x,vars); ]; ]; LocalSymbols(x,y,vars) [ 30 # MakeMultiNomial(_x - _y,vars_IsList) <-- MakeMultiNomial(x,vars) - MakeMultiNomial(y,vars); ]; LocalSymbols(x,n,vars) [ 30 # MakeMultiNomial(_x ^ n_IsInteger,vars_IsList)_(Contains(vars,x)) <-- MultiSingleFactor(vars,x,n); ]; LocalSymbols(x,n,vars) [ 40 # MakeMultiNomial(_x ^ n_IsPositiveInteger,vars_IsList) <-- [ Local(mult,result); Set(mult,MakeMultiNomial(x,vars)); Set(result,MakeMultiNomial(1,vars)); While(n>0) [ If(n&1 != 0, Set(result, MultiNomialMultiply(result,mult))); Set(n,n>>1); Set(mult,MultiNomialMultiply(mult,mult)); ]; result; ]; 50 # MakeMultiNomial(_x ^ (_n),vars_IsList)_(Contains(vars,x)) <-- [ //Echo({"ENTER!"}); Set(n,Simplify(n)); If(IsInteger(n), MultiSingleFactor(vars,x,n), MultiSingleFactor(vars,x^n,1) ); ]; ]; LocalSymbols(x,vars) [ 100 # MakeMultiNomial(x_IsFunction,vars_IsList) <-- [ //Echo({"Adding unknown",x}); MultiSingleFactor(vars,x,1); ]; ]; MultiNomialAdd(MultiNomial(_vars,_terms1),MultiNomial(_vars,_terms2)) <-- [ ForEach(t,terms2) MultiAddTerm(terms1,t); MultiNomial(vars,terms1); ]; MultiTermLess({_deg1,_fact1},{_deg2,_fact2}) <-- [ Local(deg); deg := deg1-deg2; While(deg != {} And Head(deg) = 0) [ deg := Tail(deg);]; ((deg = {}) And (fact1-fact2 < 0)) Or ((deg != {}) And (deg[1] < 0)); ]; 10 # MultiAddTerm(_terms,_t)_(Assoc(t[1],terms) = Empty) <-- [ Local(i,nr); nr:=Length(terms); i:=1; While (i<=nr And MultiTermLess(t,terms[i])) i++; DestructiveInsert(terms,i,t); /* DestructiveAppend(terms,t); */ ]; 20 # MultiAddTerm(_terms,_t) <-- [ Local(as); as := Assoc(t[1],terms); as[2] := as[2] + t[2]; ]; MultiNomialMultiply(MultiNomial(_vars,_terms1),MultiNomial(_vars,_terms2)) <-- [ Local(result,t); result := {}; ForEach(t1,terms1) ForEach(t2,terms2) [ t := FlatCopy(t1); t[1] := t[1] + t2[1]; t[2] := t[2] * t2[2]; MultiAddTerm(result,t); ]; MultiNomial(vars,result); ]; MultiNomialMultiply( MultiNomial(_vars,_terms1)/MultiNomial(_vars,_terms2), MultiNomial(_vars,_terms3)/MultiNomial(_vars,_terms4)) <-- [ MultiNomialMultiply(MultiNomial(vars,terms1),MultiNomial(vars,terms3))/ MultiNomialMultiply(MultiNomial(vars,terms2),MultiNomial(vars,terms4)); ]; MultiNomialMultiply( MultiNomial(_vars,_terms1)/MultiNomial(_vars,_terms2), MultiNomial(_vars,_terms3)) <-- [ MultiNomialMultiply(MultiNomial(vars,terms1),MultiNomial(vars,terms3))/ MultiNomial(vars,terms2); ]; MultiNomialMultiply( MultiNomial(_vars,_terms3), MultiNomial(_vars,_terms1)/MultiNomial(_vars,_terms2)) <-- [ MultiNomialMultiply(MultiNomial(vars,terms1),MultiNomial(vars,terms3))/ MultiNomial(vars,terms2); ]; 10 # MultiNomialMultiply(_a,_b) <-- [ Echo({"ERROR!",a,b}); Echo({"ERROR!",Type(a),Type(b)}); ]; MultiNomialNegate(MultiNomial(_vars,_terms)) <-- [ ForEach(t,terms) [t[2] := -t[2]; ]; MultiNomial(vars,terms); ]; 10 # MultiDegree(MultiNomial(_vars,{})) <-- FillList(-Infinity,Length(vars)); 20 # MultiDegree(MultiNomial(_vars,_terms)) <-- Head(Head(terms)); 10 # MultiLeadingCoef(MultiNomial(_vars,{})) <-- 0; 20 # MultiLeadingCoef(MultiNomial(_vars,_terms)) <-- (Head(terms)[2]); 10 # MultiLeadingMono(MultiNomial(_vars,{})) <-- 0; 20 # MultiLeadingMono(MultiNomial(_vars,_terms)) <-- Factorize(vars^(Head(Head(terms)))); 20 # MultiLeadingTerm(_m) <-- MultiLeadingCoef(m) * MultiLeadingMono(m); MultiVars(MultiNomial(_vars,_terms)) <-- vars; 10 # MultiLT(MultiNomial(_vars,{})) <-- MultiNomial(vars,{{FillList(0,Length(vars)),0}}); 20 # MultiLT(MultiNomial(_vars,_terms)) <-- MultiNomial(vars,{terms[1]}); 10 # MultiLM(MultiNomial(_vars,{})) <-- FillList(0,Length(vars)); 20 # MultiLM(MultiNomial(_vars,_terms)) <-- terms[1][1]; 10 # MultiLC(MultiNomial(_vars,{})) <-- 0; 20 # MultiLC(MultiNomial(_vars,_terms)) <-- terms[1][2]; 10 # MultiZero(MultiNomial(_vars,{})) <-- True; 20 # MultiZero(MultiNomial(_vars,_terms)) <-- False; DropZeroLC(MultiNomial(_vars,_terms)) <-- MultiNomial(vars,DropZeroLC(terms)); 10 # DropZeroLC({}) <-- {}; 20 # DropZeroLC(terms_IsList) <-- DropZeroLC(Head(terms),Tail(terms)); 30 # DropZeroLC({_pows,0},terms_IsList) <-- DropZeroLC(terms); 40 # DropZeroLC(_pows,terms_IsList) <-- pows:terms; /* MultiDivide : input f - a multivariate polynomial g[1 .. n] - a list of polynomials to divide by output {q[1 .. n],r} such that f = q[1]*g[1] + ... + q[n]*g[n] + r Basically quotient and remainder after division by a group of polynomials. */ 20 # MultiDivide(_f,g_IsList) <-- [ Local(i,v,q,r,nr); v:=MultiExpressionList(f+Sum(g)); f:=MakeMultiNomial(f,v); nr := Length(g); For(i:=1,i<=nr,i++) [ g[i] := MakeMultiNomial(g[i],v); ]; {q,r}:=MultiDivide(f,g); q:=MapSingle("NormalForm",q); r:=NormalForm(r); {q,r}; ]; 10 # MultiDivide(f_IsMulti,g_IsList) <-- [ Local(i,nr,q,r,p,v,finished); Set(nr, Length(g)); Set(v, MultiVars(f)); Set(q, FillList(0,nr)); Set(r, 0); Set(p, f); Set(finished,MultiZero(p)); Local(plt,glt); While (Not finished) [ //Set(g,LocalSymbols(a,b)(Subst(a,b)g)); //Echo({"p = ",p}); //hier //Echo({""}); //hier Set(plt, MultiLT(p)); For(i:=1,i<=nr,i++) [ Set(glt, MultiLT(g[i])); if (MultiLM(glt) = MultiLM(plt) Or MultiTermLess({MultiLM(glt),1}, {MultiLM(plt),1})) if (Select({{n},n<0},MultiLM(plt)-MultiLM(glt)) = {}) [ Local(ff); Set(ff, MultiNomial(v,{{MultiLM(plt)-MultiLM(glt),MultiLC(plt)/MultiLC(glt)}})); q[i] := q[i] + ff; Set(p, p - ff*g[i]); Set(i,nr+2); ]; ]; //Echo({"p3 = ",p}); //hier //Echo({"r = ",r}); If (i = nr+1, [ Set(r, r + LocalSymbols(a,b)(Subst(a,b)plt)); Set(p, p - LocalSymbols(a,b)(Subst(a,b)plt)); ]); //Echo({"r = ",r}); //Echo({"p4 = ",p}); //hier Set(p,DropZeroLC(p)); // Set(p,DropZeroLC(LocalSymbols(a,b)(Subst(a,b)p))); //Echo({"p5 = ",p}); //hier Set(finished,MultiZero(p)); ]; {q,r}; ]; //TODO optimize this! keeps on converting to and from internal format! 10 # MultiGcd( 0,_g) <-- 1; 10 # MultiGcd(_f, 0) <-- 1; 20 # MultiGcd(_f,_g) <-- [ Local(v); v:=MultiExpressionList(f+g); //hier NormalForm(MultiGcd(MakeMultiNomial(f,v),MakeMultiNomial(g,v))); ]; 5 # MultiGcd(f_IsMulti,g_IsMulti)_(MultiTermLess({MultiLM(f),1},{MultiLM(g),1})) <-- [ MultiGcd(g,f); ]; 5 # MultiGcd(MultiNomial(_vars,_terms),g_IsMulti)_(MultiLM(MultiNomial(vars,terms)) = MultiLM(g)) <-- MultiNomial(vars,{{FillList(0,Length(vars)),1}}); 5 # MultiGcd(MultiNomial(_vars,_terms),g_IsMulti)_(Select({{n},n<0},MultiLM(MultiNomial(vars,terms))-MultiLM(g)) != {}) <-- MultiNomial(vars,{{FillList(0,Length(vars)),1}}); 5 # MultiGcd(MultiNomial(_vars,_terms),g_IsMulti)_(NormalForm(g) = 0) <-- MultiNomial(vars,{{FillList(0,Length(vars)),1}}); 10 # MultiGcd(f_IsMulti,g_IsMulti) <-- [ LocalSymbols(a) [ Set(f,Subst(a,a)f); Set(g,Subst(a,a)g); ]; Local(new); While(g != 0) [ Set(new, MultiDivide(f,{g})[2]); Set(f,g); Set(g,new); ]; MultiPrimitivePart(f); ]; MultiDivTerm(MultiNomial(_vars,{_term1}),MultiNomial(_vars,{_term2})) <-- [ MultiNomial(vars,{{term1[1]-term2[1],term1[2]/term2[2]}}); ]; MultiS(_g,_h,MultiNomial(_vars,_terms)) <-- [ Local(gamma); gamma :=Max(MultiDegree(g),MultiDegree(h)); Local(result); result := MultiDivTerm(MultiNomial(vars,{{gamma,1}}),MultiLT(g))*g - MultiDivTerm(MultiNomial(vars,{{gamma,1}}),MultiLT(h))*h; result; ]; /* Groebner : Calculate the Groebner basis of a set of polynomials. Nice example of its power is In> TableForm(Groebner({x*(y-1),y*(x-1)})) x*y-x x*y-y y-x y^2-y In> Factor(y^2-y) Out> y*(y-1); From which you can see that x = y, and x^2 = x so x is 0 or 1. */ Groebner(f_IsList) <-- [ Local(vars,i,j,S,nr,r); nr:=Length(f); vars:=VarList(f); For(i:=1,i<=nr,i++) [ f[i] := MakeMultiNomial(f[i],vars); ]; S:={}; For(i:=1,i ",f}); Local(ms); ms:=MultiSimp(f); Echo({"Out> ",ms}); ]; //try(D(x,4)Exp((-x^2)/2)); //try(1/(1+1/x)); //try(1/(1+1/(1+1/x))); tryfull():= [ exp:=1/x; While(True) [ try(exp); exp:=1/(1+exp); ]; ]; //try(1/(1+1/(1+1/(1+1/x)))); //try(tst); //(-2 m^3 y)/M^2-(m^2 y)/M+m y+1/(1/(2 y m+(-6 y m^2)/M+(4 m^3 y)/M^2)+1/((-6 m^2 y)/M+(4 m^3 y)/M^2+2 y m)) */