/* todos for univariates: - Factorize */ /* we need to use "predicates" because we are adding to IsZero */ Use("predicates.ys"); RuleBase("NormalForm",{expression}); Rule("NormalForm",1,1000,True) expression; 0 # NormalForm(UniVariate(_var,_first,_coefs)) <-- ExpandUniVariate(var,first,coefs); /*TODO remove Rule("NormalForm",1,0,Type(expression) = "UniVariate") ExpandUniVariate(expression[1],expression[2],expression[3]); */ Function("ExpandUniVariate",{var,first,coefs}) [ Local(result,i); result:=0; For(i:=Length(coefs),i>0,i--) [ Local(term); term:=NormalForm(coefs[i])*var^(first+i-1); result:=result+term; ]; result; ]; 10 # IsUniVar(UniVariate(_var,_first,_coefs)) <-- True; 20 # IsUniVar(_anything) <-- False; /*TODO remove? Function("IsUniVar",{expr}) Type(expr) = "UniVariate"; */ RuleBase("UniVariate",{var,first,coefs}); Rule("UniVariate",3,10,Length(coefs)>0 And coefs[1]=0) UniVariate(var,first+1,Tail(coefs)); Rule("UniVariate",3,1000,IsComplex(var) Or IsList(var)) ExpandUniVariate(var,first,coefs); 20 # IsZero(UniVariate(_var,_first,_coefs)) <-- IsZeroVector(coefs); RuleBase("Degree",{expr}); Rule("Degree",1,0, IsUniVar(expr)) [ Local(i,min,max); min:=expr[2]; max:=min+Length(expr[3]); i:=max; While(i >= min And IsZero(Coef(expr,i))) i--; i; ]; 10 # Degree(_poly) <-- Degree(MakeUni(poly)); 10 # Degree(_poly,_var) <-- Degree(MakeUni(poly,var)); 500 # UniVariate(_var,_f1,_c1) + UniVariate(_var,_f2,_c2) <-- [ Local(from,result); Local(curl,curr,left,right); Set(curl, f1); Set(curr, f2); Set(left, c1); Set(right, c2); Set(result, {}); Set(from, Min(curl,curr)); While(MathAnd(LessThan(curl,curr),left != {})) [ DestructiveAppend(result,Head(left)); Set(left,Tail(left)); Set(curl,MathAdd(curl,1)); ]; While(LessThan(curl,curr)) [ DestructiveAppend(result,0); Set(curl,MathAdd(curl,1)); ]; While(MathAnd(LessThan(curr,curl), right != {})) [ DestructiveAppend(result,Head(right)); Set(right,Tail(right)); Set(curr,MathAdd(curr,1)); ]; While(LessThan(curr,curl)) [ DestructiveAppend(result,0); Set(curr,MathAdd(curr,1)); ]; While(MathAnd(left != {}, right != {})) [ DestructiveAppend(result,Head(left)+Head(right)); Set(left, Tail(left)); Set(right, Tail(right)); ]; While(left != {}) [ DestructiveAppend(result,Head(left)); Set(left, Tail(left)); ]; While(right != {}) [ DestructiveAppend(result,Head(right)); Set(right, Tail(right)); ]; UniVariate(var,from,result); ]; 200 # UniVariate(_var,_first,_coefs) + a_IsNumber <-- UniVariate(var,first,coefs) + UniVariate(var,0,{a}); 200 # a_IsNumber + UniVariate(_var,_first,_coefs) <-- UniVariate(var,first,coefs) + UniVariate(var,0,{a}); /*TODO remove? 200 # aLeft_IsUniVar + aRight_IsNumber <-- aRight+aLeft; 200 # aLeft_IsNumber + aRight_IsUniVar <-- UniVariate(aRight[1],0,{aLeft})+aRight; */ 200 # - UniVariate(_var,_first,_coefs) <-- UniVariate(var,first,-coefs); /*TODO remove? 200 # - (aLeft_IsUniVar) <-- Apply("UniVariate",{aLeft[1],aLeft[2],-(aLeft[3])}); */ 200 # aLeft_IsUniVar - aRight_IsUniVar <-- [ Local(from,result); Local(curl,curr,left,right); curl:=aLeft[2]; curr:=aRight[2]; left:=aLeft[3]; right:=aRight[3]; result:={}; from:=Min(curl,curr); While(curl0) [ Set(b,b*b); If (MathNot(Equals(BitAnd(m,n), 0)),Set(b,b*a)); Set(m, ShiftRight(m,1)); ]; b; ]; 200 # aLeft_IsUniVar ^ aRight_IsPositiveInteger <-- RepeatedSquaresMultiply(aLeft,aRight); /*TODO this can be made twice as fast!*/ 200 # (aLeft_IsUniVar * _aRight)_(Not(Contains(VarList(aRight),aLeft[1]))) <-- aRight*aLeft; 200 # (_factor * UniVariate(_var,_first,_coefs))_(Not(Contains(VarList(factor),var))) <-- UniVariate(var,first,coefs*factor); 200 # (UniVariate(_var,_first,_coefs)/_factor)_(Not(Contains(VarList(factor),var))) <-- UniVariate(var,first,coefs/factor); MaxUniOrder:=3000; Function("SetOrder",{order}) MaxUniOrder:=order; ShiftUniVar(UniVariate(_var,_first,_coefs),_fact,_shift) <-- UniVariate(var,first+shift,fact*coefs); 200 # UniVariate(_var,_f1,_c1) * UniVariate(_var,_f2,_c2) <-- [ Local(i,j,n,shifted,result); Set(result,MakeUni(0,var)); Set(n,Length(c1)); For(i:=1,i<=n,i++) [ Set(result,result+ShiftUniVar(UniVariate(var,f2,c2),MathNth(c1,i),f1+i-1)); ]; result; ]; /*TODO remove? Function("ShiftUniVar",{uni,fact,shift}) [ Apply("UniVariate",{uni[1],uni[2]+shift,fact*(uni[3])}); ]; 200 # (aLeft_IsUniVar * aRight_IsUniVar)_(aLeft[1] = aRight[1]) <-- [ Local(i,j,n,shifted,result); Set(result,MakeUni(0,aLeft[1])); Set(n,Length(aLeft[3])); For(i:=1,i<=n,i++) [ result:=result+ShiftUniVar(aRight,aLeft[3][i],aLeft[2]+i-1); ]; result; ]; */ 5 # Coef(uv_IsUniVar,order_IsList) <-- [ Local(result); result:={}; ForEach(item,order) [ DestructiveAppend(result,Coef(uv,item)); ]; result; ]; 10 # Coef(uv_IsUniVar,_order)_(order=uv[2]+Length(uv[3])) <-- 0; 20 # Coef(uv_IsUniVar,_order) <-- uv[3][(order-uv[2])+1]; 30 # Coef(uv_CanBeUni,_order) <-- Coef(MakeUni(uv),order); Function("Coef",{expression,var,order}) NormalForm(Coef(MakeUni(expression,var),order)); 10 # LeadingCoef(uv_IsUniVar) <-- Coef(uv,Degree(uv)); 20 # LeadingCoef(uv_CanBeUni) <-- [ Local(uvi); uvi:=MakeUni(uv); Coef(uvi,Degree(uvi)); ]; 10 # LeadingCoef(uv_CanBeUni,_var) <-- [ Local(uvi); uvi:=MakeUni(uv,var); Coef(uvi,var,Degree(uvi)); ]; Function("UniTaylor",{taylorfunction,taylorvariable,taylorat,taylororder}) [ Local(n,result,dif,polf); result:={}; [ MacroLocal(taylorvariable); MacroSet(taylorvariable,taylorat); DestructiveAppend(result,Eval(taylorfunction)); ]; dif:=taylorfunction; polf:=(taylorvariable-taylorat); For(n:=1,n<=taylororder,n++) [ dif:= Deriv(taylorvariable) dif; MacroLocal(taylorvariable); MacroSet(taylorvariable,taylorat); DestructiveAppend(result,(Eval(dif)/n!)); ]; UniVariate(taylorvariable,0,result); ]; Function("MakeUni",{expression}) MakeUni(expression,VarList(expression)); /* Convert normal form to univariate expression */ RuleBase("MakeUni",{expression,var}); 1 # MakeUni(_expr,{}) <-- UniVariate(dummyvar,0,{expression}); 2 # MakeUni(_expr,var_IsList) <-- [ Local(result,item); result:=expression; ForEach(item,var) [ result:=MakeUni(result,item); ]; result; ]; 10 # MakeUni(UniVariate(_var,_first,_coefs),_var) <-- UniVariate(var,first,coefs); 20 # MakeUni(UniVariate(_v,_first,_coefs),_var) <-- [ Local(reslist,item); reslist:={}; ForEach(item,expression[3]) [ If(IsFreeOf(item,var), DestructiveAppend(reslist,item), DestructiveAppend(reslist,MakeUni(item,var)) ); ]; UniVariate(expression[1],expression[2],reslist); ]; /*TODO remove? Rule("MakeUni",2,10,Type(expression) = "UniVariate") [ Local(reslist,item); reslist:={}; ForEach(item,expression[3]) [ If(IsFreeOf(item,var), DestructiveAppend(reslist,item), DestructiveAppend(reslist,MakeUni(item,var)) ); ]; Apply("UniVariate",{expression[1],expression[2],reslist}); ]; */ 20 # MakeUni(_expression,_var)_(IsFreeOf(expression,var)) <-- UniVariate(var,0,{expression}); 30 # MakeUni(_var,_var) <-- UniVariate(var,1,{1}); 30 # MakeUni(_a + _b,_var) <-- MakeUni(a,var) + MakeUni(b,var); 30 # MakeUni(_a - _b,_var) <-- MakeUni(a,var) - MakeUni(b,var); 30 # MakeUni( - _b,_var) <-- - MakeUni(b,var); 30 # MakeUni(_a * _b,_var) <-- MakeUni(a,var) * MakeUni(b,var); 30 # MakeUni(_a ^ n_IsInteger,_var) <-- MakeUni(a,var) ^ n; 30 # MakeUni(_a / _b,_var)_(IsFreeOf(b,var)) <-- MakeUni(a,var) * (1/b); /*TODO remove? Rule("MakeUni",2,20,IsFreeOf(expression,var)) UniVariate(var,0,{expression}); Rule("MakeUni",2,30,expression=var) UniVariate(var,1,{1}); Rule("MakeUni",2,30,Type(expression) = "+") MakeUni(expression[1],var)+MakeUni(expression[2],var); Rule("MakeUni",2,30,Type(expression) = "*") MakeUni(expression[1],var)*MakeUni(expression[2],var); Rule("MakeUni",2,30,Type(expression) = "^" And IsInteger(expression[2])) MakeUni(expression[1],var)^expression[2]; Rule("MakeUni",2,30,Type(expression) = "-" And NrArgs(expression) = 1) -(MakeUni(expression[1],var)); Rule("MakeUni",2,30,Type(expression) = "/" And Not(Contains(VarList(expression[2]),var))) MakeUni(expression[1],var)*(1/expression[2]); Rule("MakeUni",2,30,Type(expression) = "-" And NrArgs(expression) = 2) MakeUni(expression[1],var)-MakeUni(expression[2],var); */ 0 # Div(n_IsUniVar,m_IsUniVar)_(Degree(n) < Degree(m)) <-- 0; 0 # Mod(n_IsUniVar,m_IsUniVar)_(Degree(n) < Degree(m)) <-- n; 1 # Div(n_IsUniVar,m_IsUniVar)_ (n[1] = m[1] And Degree(n) >= Degree(m)) <-- [ UniVariate(n[1],0, UniDivide(Concat(ZeroVector(n[2]),n[3]), Concat(ZeroVector(m[2]),m[3]))[1]); ]; 1 # Mod(n_IsUniVar,m_IsUniVar)_ (n[1] = m[1] And Degree(n) >= Degree(m)) <-- [ UniVariate(n[1],0, UniDivide(Concat(ZeroVector(n[2]),n[3]), Concat(ZeroVector(m[2]),m[3]))[2]); ]; /* division algo: (for zero-base univariates:) */ Function("UniDivide",{u,v}) [ Local(m,n,q,r,k,j); m := Length(u)-1; n := Length(v)-1; While (m>0 And IsZero(u[m+1])) m--; While (n>0 And IsZero(v[n+1])) n--; q := ZeroVector(m-n+1); r := FlatCopy(u); /* (m should be >= n) */ For(k:=m-n,k>=0,k--) [ q[k+1] := r[n+k+1]/v[n+1]; For (j:=n+k-1,j>=k,j--) [ r[j+1] := r[j+1] - q[k+1]*v[j-k+1]; ]; ]; Local(end); end:=Length(r); While (end>n) [ DestructiveDelete(r,end); end:=end-1; ]; {q,r}; ]; DropEndZeroes(list):= [ Local(end); end:=Length(list); While(list[end] = 0) [ DestructiveDelete(list,end); end:=end-1; ]; ]; Function("UniGcd",{u,v}) [ Local(l,div,mod,m); DropEndZeroes(u); DropEndZeroes(v); /* If(Length(v)>Length(u), [ Locap(swap); swap:=u; u:=v; v:=swap; ] ); If(Length(u)=Length(v) And v[Length(v)] > u[Length(u)], [ Locap(swap); swap:=u; u:=v; v:=swap; ] ); */ l:=UniDivide(u,v); div:=l[1]; mod:=l[2]; DropEndZeroes(mod); m := Length(mod); /* Echo({"v,mod = ",v,mod}); */ /* If(m <= 1, */ If(m = 0, v, /* v/v[Length(v)], */ UniGcd(v,mod)); ]; 0 # Gcd(n_IsUniVar,m_IsUniVar)_ (n[1] = m[1] And Degree(n) < Degree(m)) <-- Gcd(m,n); 1 # Gcd(nn_IsUniVar,mm_IsUniVar)_ (nn[1] = mm[1] And Degree(nn) >= Degree(mm)) <-- [ UniVariate(nn[1],0, UniGcd(Concat(ZeroVector(nn[2]),nn[3]), Concat(ZeroVector(mm[2]),mm[3]))); ]; RuleBase("PSolve",{uni}); Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 1) -Coef(uni,0)/Coef(uni,1); Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 2) [ Local(a,b,c,d); c:=Coef(uni,0); b:=Coef(uni,1); a:=Coef(uni,2); d:=b*b-4*a*c; {(-b+Sqrt(d))/(2*a),(-b-Sqrt(d))/(2*a)}; ]; Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 3 ) [ Local(p,q,r,w,ww,a,b); Local(coef0,coef1,coef3,adjust); /* Get coefficients for a new polynomial, such that the coefficient of degree 2 is zero: Take f(x)=a0+a1*x+a2*x^2+a3*x^3 and substitute x = x' + adjust This gives g(x) = b0+b1*x+b2*x^2+b3*x^3 where b3 = a3; b2 = 0 => adjust = (-a2)/(3*a3); b1 = 2*a2*adjust+3*a3*adjust^2+a1; b0 = a2*adjust^2+a3*adjust^3+adjust*a1+a0; After solving g(x') = 0, return x = x' + adjust. */ adjust := (-Coef(uni,2))/(3*Coef(uni,3)); coef3 := Coef(uni,3); coef1 := 2*Coef(uni,2)*adjust+3*Coef(uni,3)*adjust^2+Coef(uni,1); coef0 := Coef(uni,2)*adjust^2+Coef(uni,3)*adjust^3+ adjust*Coef(uni,1)+Coef(uni,0); p:=coef3; q:=coef1/p; r:=coef0/p; w:=Complex(-1/2,Sqrt(3/4)); ww:=Complex(-1/2,-Sqrt(3/4)); /* Equation is xxx + qx + r = 0 */ /* Let x = a + b a^3 + b^3 + 3(aab + bba) + q(a + b) + r = 0 a^3 + b^3 + (3ab+q)x + r = 0 Let 3ab+q = 0. This is permissible, for we can still find a+b == x a^3 + b^3 = -r (ab)^3 = -q^3/27 So a^3 and b^3 are the roots of t^2 + rt - q^3/27 = 0 Let a^3 = -r/2 + Sqrt(q^3/27+ rr/4) b^3 = -r/2 - Sqrt(q^3/27+ rr/4) Therefore there are three values for each of a and b. Clearly if ab = -q/3 is true then (wa)(wwb) == (wb)(wwa) == -q/3 */ a:=(-r/2 + Sqrt(q^3/27+ r*r/4))^(1/3); b:=(-r/2 - Sqrt(q^3/27+ r*r/4))^(1/3); {a+b+adjust,w*a+ww*b+adjust,ww*a+w*b+adjust}; ]; Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 4 ) [ Local(coef4,a1,a2,a3,a4,y,y1,z); coef4:=Coef(uni,4); a1:=Coef(uni,3)/coef4; a2:=Coef(uni,2)/coef4; a3:=Coef(uni,1)/coef4; a4:=Coef(uni,0)/coef4; y1:=Head(PSolve(y^3-a2*y^2+(a1*a3-4*a4)*y+(4*a2*a4-a3^2-a1^2*a4),y)); Concat(PSolve(z^2+1/2*(a1+Sqrt(a1^2-4*a2+4*y1))*z+1/2*(y1-Sqrt(y1^2-4*a4)),z), PSolve(z^2+1/2*(a1-Sqrt(a1^2-4*a2+4*y1))*z+1/2*(y1+Sqrt(y1^2-4*a4)),z)); ]; Function("PSolve",{uni,var}) [ PSolve(MakeUni(uni,var)); ]; /* Generate a random polynomial */ RandomPoly(_var,_degree,_coefmin,_coefmax) <-- NormalForm(UniVariate(var,0,RandomIntegerVector(degree+1,coefmin,coefmax))); /* CanBeUni returns whether the function can be converted to a * univariate, with respect to a variable. */ Function("CanBeUni",{expression}) CanBeUni(expression,VarList(expression)); /* Accepting an expression as being convertable to univariate */ /* Dealing wiht a list of variables. The poly should be expandable * to each of these variables (smells like tail recursion) */ 10 # CanBeUni(_expression,{}) <-- True; 20 # CanBeUni(_expression,var_IsList) <-- CanBeUni(expression,Head(var)) And CanBeUni(expression,Tail(var)); /* Atom can always be a polynom to any variable */ 30 # CanBeUni(expression_IsAtom,_var) <-- True; 35 # CanBeUni(_expression,_var)_IsFreeOf(expression,var) <-- True; /* Other patterns supported. */ 40 # CanBeUni(_x + _y,_var) <-- CanBeUni(x,var) And CanBeUni(y,var); 40 # CanBeUni(_x - _y,_var) <-- CanBeUni(x,var) And CanBeUni(y,var); 40 # CanBeUni( + _y,_var) <-- CanBeUni(y,var); 40 # CanBeUni( - _y,_var) <-- CanBeUni(y,var); 40 # CanBeUni(_x * _y,_var) <-- CanBeUni(x,var) And CanBeUni(y,var); 40 # CanBeUni(_x / _y,_var) <-- CanBeUni(x,var) And IsFreeOf(y,var); /* Special case again: raising powers */ 40 # CanBeUni(_x ^ y_IsInteger,_var)_CanBeUni(x,var) <-- True; 41 # CanBeUni(_x ^ _y,_var)_(IsFreeOf(x,var) And IsFreeOf(y,var)) <-- True; 50 # CanBeUni(UniVariate(_var,_first,_coefs),_var) <-- True; 1000 # CanBeUni(_f,_var)_(Not(IsFreeOf(f,var))) <-- False; 1001 # CanBeUni(_f,_var) <-- True; 10 # Content(UniVariate(_var,_first,_coefs)) <-- Gcd(coefs)*var^first; 20 # Content(poly_CanBeUni) <-- NormalForm(Content(MakeUni(poly))); 10 # PrimitivePart(UniVariate(_var,_first,_coefs)) <-- UniVariate(var,0,coefs/Gcd(coefs)); 20 # PrimitivePart(poly_CanBeUni) <-- NormalForm(PrimitivePart(MakeUni(poly))); 10 # Monic(UniVariate(_var,_first,_coefs)) <-- [ DropEndZeroes(coefs); UniVariate(var,first,coefs/coefs[Length(coefs)]); ]; 20 # Monic(poly_CanBeUni) <-- NormalForm(Monic(MakeUni(poly))); 30 # Monic(_poly,_var)_CanBeUni(poly,var) <-- NormalForm(Monic(MakeUni(poly,var))); 10 # BigOh(UniVariate(_var,_first,_coefs),_var,_degree) <-- [ While(first+Length(coefs)>=(degree+1) And Length(coefs)>0) DestructiveDelete(coefs,Length(coefs)); UniVariate(var,first,coefs); ]; 20 # BigOh(_uv,_var,_degree)_CanBeUni(uv,var) <-- NormalForm(BigOh(MakeUni(uv,var),var,degree)); Horner(_e,_v) <-- [ Local(uni,coefs,result); uni := MakeUni(e,v); coefs:=DestructiveReverse(uni[3]); result:=0; While(coefs != {}) [ result := result*v; result := result+Head(coefs); coefs := Tail(coefs); ]; result:=result*v^uni[2]; result; ]; DivPoly(_A,_B,_var,_deg) <-- [ Local(a,b,c,i,j,denom); b:=MakeUni(B,var); denom:=Coef(b,0); if (denom = 0) [ Local(f); f:=Content(b); b:=PrimitivePart(b); A:=Simplify(A/f); denom:=Coef(b,0); ]; a:=MakeUni(A,var); c:=FillList(0,deg+1); For(i:=0,i<=deg,i++) [ Local(sum,j); sum:=0; For(j:=0,j