;; Calculator for GNU Emacs, part II [calc-poly.el];; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.;; Written by Dave Gillespie, daveg@synaptics.com.;; This file is part of GNU Emacs.;; GNU Emacs is distributed in the hope that it will be useful,;; but WITHOUT ANY WARRANTY. No author or distributor;; accepts responsibility to anyone for the consequences of using it;; or for whether it serves any particular purpose or works at all,;; unless he says so in writing. Refer to the GNU Emacs General Public;; License for full details.;; Everyone is granted permission to copy, modify and redistribute;; GNU Emacs, but only under the conditions described in the;; GNU Emacs General Public License. A copy of this license is;; supposed to have been given to you along with GNU Emacs so you;; can know your rights and responsibilities. It should be in a;; file named COPYING. Among other things, the copyright notice;; and this notice must be preserved on all copies.;; This file is autoloaded from calc-ext.el.(require'calc-ext)(require'calc-macs)(defuncalc-Need-calc-poly()nil)(defuncalcFunc-pcont(expr&optionalvar)(cond((Math-primpexpr)(cond((Math-zeropexpr)1)((Math-messy-integerpexpr)(math-truncexpr))((Math-objectpexpr)expr)((or(equalexprvar)(notvar))1)(texpr)))((eq(carexpr)'*)(math-mul(calcFunc-pcont(nth1expr)var)(calcFunc-pcont(nth2expr)var)))((eq(carexpr)'/)(math-div(calcFunc-pcont(nth1expr)var)(calcFunc-pcont(nth2expr)var)))((and(eq(carexpr)'^)(Math-natnump(nth2expr)))(math-pow(calcFunc-pcont(nth1expr)var)(nth2expr)))((memq(carexpr)'(negpolar))(calcFunc-pcont(nth1expr)var))((conspvar)(let((p(math-is-polynomialexprvar)))(ifp(let((lead(nth(1-(lengthp))p))(cont(math-poly-gcd-listp)))(if(math-guess-if-neglead)(math-negcont)cont))1)))((memq(carexpr)'(+-cplxsdev))(let((cont(calcFunc-pcont(nth1expr)var)))(if(eqcont1)1(let((c2(calcFunc-pcont(nth2expr)var)))(if(and(math-negpcont)(if(eq(carexpr)'-)(math-pospc2)(math-negpc2)))(math-neg(math-poly-gcdcontc2))(math-poly-gcdcontc2))))))(varexpr)(t1)))(defuncalcFunc-pprim(expr&optionalvar)(let((cont(calcFunc-pcontexprvar)))(if(math-equal-intcont1)expr(math-poly-div-exactexprcontvar))))(defunmath-div-poly-const(exprc)(cond((memq(car-safeexpr)'(+-))(list(carexpr)(math-div-poly-const(nth1expr)c)(math-div-poly-const(nth2expr)c)))(t(math-divexprc))))(defuncalcFunc-pdeg(expr&optionalvar)(if(Math-zeropexpr)'(neg(varinfvar-inf))(ifvar(or(math-polynomial-pexprvar)(math-reject-argexpr"Expected a polynomial"))(math-poly-degreeexpr))))(defunmath-poly-degree(expr)(cond((Math-primpexpr)(if(eq(car-safeexpr)'var)10))((eq(carexpr)'neg)(math-poly-degree(nth1expr)))((eq(carexpr)'*)(+(math-poly-degree(nth1expr))(math-poly-degree(nth2expr))))((eq(carexpr)'/)(-(math-poly-degree(nth1expr))(math-poly-degree(nth2expr))))((and(eq(carexpr)'^)(natnump(nth2expr)))(*(math-poly-degree(nth1expr))(nth2expr)))((memq(carexpr)'(+-))(max(math-poly-degree(nth1expr))(math-poly-degree(nth2expr))))(t1)))(defuncalcFunc-plead(exprvar)(cond((eq(car-safeexpr)'*)(math-mul(calcFunc-plead(nth1expr)var)(calcFunc-plead(nth2expr)var)))((eq(car-safeexpr)'/)(math-div(calcFunc-plead(nth1expr)var)(calcFunc-plead(nth2expr)var)))((and(eq(car-safeexpr)'^)(math-natnump(nth2expr)))(math-pow(calcFunc-plead(nth1expr)var)(nth2expr)))((Math-primpexpr)(if(equalexprvar)1expr))(t(let((p(math-is-polynomialexprvar)))(if(cdrp)(nth(1-(lengthp))p)1)))));;; Polynomial quotient, remainder, and GCD.;;; Originally by Ove Ewerlid (ewerlid@mizar.DoCS.UU.SE).;;; Modifications and simplifications by daveg.(setqmath-poly-modulus1);;; Return gcd of two polynomials(defuncalcFunc-pgcd(pnpd)(if(math-any-floatspn)(math-reject-argpn"Coefficients must be rational"))(if(math-any-floatspd)(math-reject-argpd"Coefficients must be rational"))(let((calc-prefer-fract)(math-poly-modulus(math-poly-moduluspnpd)))(math-poly-gcdpnpd)));;; Return only quotient to top of stack (nil if zero)(defuncalcFunc-pdiv(pnpd&optionalbase)(let*((calc-prefer-fract)(math-poly-modulus(math-poly-moduluspnpd))(res(math-poly-divpnpdbase)))(setqcalc-poly-div-remainder(cdrres))(carres)));;; Return only remainder to top of stack(defuncalcFunc-prem(pnpd&optionalbase)(let((calc-prefer-fract)(math-poly-modulus(math-poly-moduluspnpd)))(cdr(math-poly-divpnpdbase))))(defuncalcFunc-pdivrem(pnpd&optionalbase)(let*((calc-prefer-fract)(math-poly-modulus(math-poly-moduluspnpd))(res(math-poly-divpnpdbase)))(list'vec(carres)(cdrres))))(defuncalcFunc-pdivide(pnpd&optionalbase)(let*((calc-prefer-fract)(math-poly-modulus(math-poly-moduluspnpd))(res(math-poly-divpnpdbase)))(math-add(carres)(math-div(cdrres)pd))));;; Multiply two terms, expanding out products of sums.(defunmath-mul-thru(lhsrhs)(if(memq(car-safelhs)'(+-))(list(carlhs)(math-mul-thru(nth1lhs)rhs)(math-mul-thru(nth2lhs)rhs))(if(memq(car-saferhs)'(+-))(list(carrhs)(math-mul-thrulhs(nth1rhs))(math-mul-thrulhs(nth2rhs)))(math-mullhsrhs))))(defunmath-div-thru(numden)(if(memq(car-safenum)'(+-))(list(carnum)(math-div-thru(nth1num)den)(math-div-thru(nth2num)den))(math-divnumden)));;; Sort the terms of a sum into canonical order.(defunmath-sort-terms(expr)(if(memq(car-safeexpr)'(+-))(math-list-to-sum(sort(math-sum-to-listexpr)(function(lambda(ab)(math-beforep(cara)(carb))))))expr))(defunmath-list-to-sum(lst)(if(cdrlst)(list(if(cdr(carlst))'-'+)(math-list-to-sum(cdrlst))(car(carlst)))(if(cdr(carlst))(math-neg(car(carlst)))(car(carlst)))))(defunmath-sum-to-list(tree&optionalneg)(cond((eq(car-safetree)'+)(nconc(math-sum-to-list(nth1tree)neg)(math-sum-to-list(nth2tree)neg)))((eq(car-safetree)'-)(nconc(math-sum-to-list(nth1tree)neg)(math-sum-to-list(nth2tree)(notneg))))(t(list(constreeneg)))));;; Check if the polynomial coefficients are modulo forms.(defunmath-poly-modulus(expr&optionalexpr2)(or(math-poly-modulus-recexpr)(andexpr2(math-poly-modulus-recexpr2))1))(defunmath-poly-modulus-rec(expr)(if(and(eq(car-safeexpr)'mod)(Math-natnump(nth2expr)))(list'mod1(nth2expr))(and(memq(car-safeexpr)'(+-*/))(or(math-poly-modulus-rec(nth1expr))(math-poly-modulus-rec(nth2expr))))));;; Divide two polynomials. Return (quotient . remainder).(defunmath-poly-div(uv&optionalmath-poly-div-base)(ifmath-poly-div-base(math-do-poly-divuv)(math-do-poly-div(calcFunc-expandu)(calcFunc-expandv))))(setqmath-poly-div-basenil)(defunmath-poly-div-exact(uv&optionalbase)(let((res(math-poly-divuvbase)))(if(eq(cdrres)0)(carres)(math-reject-arg(list'vecuv)"Argument is not a polynomial"))))(defunmath-do-poly-div(uv)(cond((math-constpu)(if(math-constpv)(cons(math-divuv)0)(cons0u)))((math-constpv)(cons(if(eqv1)u(if(memq(car-safeu)'(+-))(math-add-or-sub(math-poly-div-exact(nth1u)v)(math-poly-div-exact(nth2u)v)nil(eq(caru)'-))(math-divuv)))0))((Math-equaluv)(consmath-poly-modulus0))((and(math-atomic-factorpu)(math-atomic-factorpv))(cons(math-simplify(math-divuv))0))(t(let((base(ormath-poly-div-base(math-poly-div-baseuv)))vpupres)(if(or(nullbase)(null(setqvp(math-is-polynomialvbasenil'gen))))(cons0u)(setqup(math-is-polynomialubasenil'gen)res(math-poly-div-coefsupvp))(cons(math-build-polynomial-expr(carres)base)(math-build-polynomial-expr(cdrres)base)))))))(defunmath-poly-div-rec(uv)(cond((math-constpu)(math-divuv))((math-constpv)(if(eqv1)u(if(memq(car-safeu)'(+-))(math-add-or-sub(math-poly-div-rec(nth1u)v)(math-poly-div-rec(nth2u)v)nil(eq(caru)'-))(math-divuv))))((Math-equaluv)math-poly-modulus)((and(math-atomic-factorpu)(math-atomic-factorpv))(math-simplify(math-divuv)))(math-poly-div-base(math-divuv))(t(let((base(math-poly-div-baseuv))vpupres)(if(or(nullbase)(null(setqvp(math-is-polynomialvbasenil'gen))))(math-divuv)(setqup(math-is-polynomialubasenil'gen)res(math-poly-div-coefsupvp))(math-add(math-build-polynomial-expr(carres)base)(math-div(math-build-polynomial-expr(cdrres)base)v)))))));;; Divide two polynomials in coefficient-list form. Return (quot . rem).(defunmath-poly-div-coefs(uv)(cond((nullv)(math-reject-argnil"Division by zero"))((<(lengthu)(lengthv))(consnilu))((cdru)(let((qnil)(urev(reverseu))(vrev(reversev)))(while(let((qk(math-poly-div-rec(math-simplify(carurev))(carvrev)))(upurev)(vpvrev))(if(orq(not(math-zeropqk)))(setqq(consqkq)))(while(setqup(cdrup)vp(cdrvp))(setcarup(math-sub(carup)(math-mul-thruqk(carvp)))))(setqurev(cdrurev))up))(while(andurev(Math-zerop(carurev)))(setqurev(cdrurev)))(consq(nreverse(mapcar'math-simplifyurev)))))(t(cons(list(math-poly-div-rec(caru)(carv)))nil))));;; Perform a pseudo-division of polynomials. (See Knuth section 4.6.1.);;; This returns only the remainder from the pseudo-division.(defunmath-poly-pseudo-div(uv)(cond((nullv)nil)((<(lengthu)(lengthv))u)((or(cdru)(cdrv))(let((urev(reverseu))(vrev(reversev))up)(while(let((vpvrev))(setqupurev)(while(setqup(cdrup)vp(cdrvp))(setcarup(math-sub(math-mul-thru(carvrev)(carup))(math-mul-thru(carurev)(carvp)))))(setqurev(cdrurev))up)(whileup(setcarup(math-mul-thru(carvrev)(carup)))(setqup(cdrup))))(while(andurev(Math-zerop(carurev)))(setqurev(cdrurev)))(nreverse(mapcar'math-simplifyurev))))(tnil)));;; Compute the GCD of two multivariate polynomials.(defunmath-poly-gcd(uv)(cond((Math-equaluv)u)((math-constpu)(if(Math-zeropu)v(calcFunc-gcdu(calcFunc-pcontv))))((math-constpv)(if(Math-zeropv)v(calcFunc-gcdv(calcFunc-pcontu))))(t(let((base(math-poly-gcd-baseuv)))(ifbase(math-simplify(calcFunc-expand(math-build-polynomial-expr(math-poly-gcd-coefs(math-is-polynomialubasenil'gen)(math-is-polynomialvbasenil'gen))base)))(calcFunc-gcd(calcFunc-pcontu)(calcFunc-pcontu)))))))(defunmath-poly-div-list(lsta)(if(eqa1)lst(if(eqa-1)(math-mul-listlsta)(mapcar(function(lambda(x)(math-poly-div-exactxa)))lst))))(defunmath-mul-list(lsta)(if(eqa1)lst(if(eqa-1)(mapcar'math-neglst)(and(not(eqa0))(mapcar(function(lambda(x)(math-mulxa)))lst)))));;; Run GCD on all elements in a list.(defunmath-poly-gcd-list(lst)(if(or(memq1lst)(memq-1lst))(math-poly-gcd-frac-listlst)(let((gcd(carlst)))(while(and(setqlst(cdrlst))(not(eqgcd1)))(or(eq(carlst)0)(setqgcd(math-poly-gcdgcd(carlst)))))(iflst(setqlst(math-poly-gcd-frac-listlst)))gcd)))(defunmath-poly-gcd-frac-list(lst)(while(andlst(not(eq(car-safe(carlst))'frac)))(setqlst(cdrlst)))(iflst(let((denom(nth2(carlst))))(while(setqlst(cdrlst))(if(eq(car-safe(carlst))'frac)(setqdenom(calcFunc-lcmdenom(nth2(carlst))))))(list'frac1denom))1));;; Compute the GCD of two monovariate polynomial lists.;;; Knuth section 4.6.1, algorithm C.(defunmath-poly-gcd-coefs(uv)(let((d(math-poly-gcd(math-poly-gcd-listu)(math-poly-gcd-listv)))(g1)(h1)(z0)hhrdeltaghd)(while(anduv(Math-zerop(caru))(Math-zerop(carv)))(setqu(cdru)v(cdrv)z(1+z)))(or(eqd1)(setqu(math-poly-div-listud)v(math-poly-div-listvd)))(while(progn(setqdelta(-(lengthu)(lengthv)))(if(<delta0)(setqruuvvrdelta(-delta)))(setqr(math-poly-pseudo-divuv))(cdrr))(setquvv(math-poly-div-listr(math-mulg(math-powhdelta)))g(nth(1-(lengthu))u)h(if(<=delta1)(math-mul(math-powgdelta)(math-powh(-1delta)))(math-poly-div-exact(math-powgdelta)(math-powh(1-delta))))))(setqv(ifr(listd)(math-mul-list(math-poly-div-listv(math-poly-gcd-listv))d)))(if(math-guess-if-neg(nth(1-(lengthv))v))(setqv(math-mul-listv-1)))(while(>=(setqz(1-z))0)(setqv(cons0v)))v));;; Return true if is a factor containing no sums or quotients.(defunmath-atomic-factorp(expr)(cond((eq(car-safeexpr)'*)(and(math-atomic-factorp(nth1expr))(math-atomic-factorp(nth2expr))))((memq(car-safeexpr)'(+-/))nil)((memq(car-safeexpr)'(^neg))(math-atomic-factorp(nth1expr)))(tt)));;; Find a suitable base for dividing a by b.;;; The base must exist in both expressions.;;; The degree in the numerator must be higher or equal than the;;; degree in the denominator.;;; If the above conditions are not met the quotient is just a remainder.;;; Return nil if this is the case.(defunmath-poly-div-base(ab)(let(a-baseb-base)(and(setqa-base(math-total-polynomial-basea))(setqb-base(math-total-polynomial-baseb))(catch'return(whilea-base(let((maybe(assoc(car(cara-base))b-base)))(ifmaybe(if(>=(nth1(cara-base))(nth1maybe))(throw'return(car(cara-base))))))(setqa-base(cdra-base)))))));;; Same as above but for gcd algorithm.;;; Here there is no requirement that degree(a) > degree(b).;;; Take the base that has the highest degree considering both a and b.;;; ("a^20+b^21+x^3+a+b", "a+b^2+x^5+a^22+b^10") --> (a 22)(defunmath-poly-gcd-base(ab)(let(a-baseb-base)(and(setqa-base(math-total-polynomial-basea))(setqb-base(math-total-polynomial-baseb))(catch'return(while(anda-baseb-base)(if(>(nth1(cara-base))(nth1(carb-base)))(if(assoc(car(cara-base))b-base)(throw'return(car(cara-base)))(setqa-base(cdra-base)))(if(assoc(car(carb-base))a-base)(throw'return(car(carb-base)))(setqb-base(cdrb-base)))))))));;; Sort a list of polynomial bases.(defunmath-sort-poly-base-list(lst)(sortlst(function(lambda(ab)(or(>(nth1a)(nth1b))(and(=(nth1a)(nth1b))(math-beforep(cara)(carb))))))));;; Given an expression find all variables that are polynomial bases.;;; Return list in the form '( (var1 degree1) (var2 degree2) ... ).;;; Note dynamic scope of mpb-total-base.(defunmath-total-polynomial-base(expr)(let((mpb-total-basenil))(math-polynomial-baseexpr'math-polynomial-p1)(math-sort-poly-base-listmpb-total-base)))(defunmath-polynomial-p1(subexpr)(or(assocsubexprmpb-total-base)(memq(carsubexpr)'(+-*/neg))(and(eq(carsubexpr)'^)(natnump(nth2subexpr)))(let*((math-poly-base-variablesubexpr)(exponent(math-polynomial-pmpb-top-exprsubexpr)))(ifexponent(setqmpb-total-base(cons(listsubexprexponent)mpb-total-base)))))nil)(defuncalcFunc-factors(expr&optionalvar)(let((math-factored-vars(ifvartnil))(math-to-listt)(calc-prefer-fract))(orvar(setqvar(math-polynomial-baseexpr)))(let((res(math-factor-finish(or(catch'factor(math-factor-expr-tryvar))expr))))(math-simplify(if(math-vectorpres)res(list'vec(list'vecres1)))))))(defuncalcFunc-factor(expr&optionalvar)(let((math-factored-varsnil)(math-to-listnil)(calc-prefer-fract))(math-simplify(math-factor-finish(ifvar(let((math-factored-varst))(or(catch'factor(math-factor-expr-tryvar))expr))(math-factor-exprexpr))))))(defunmath-factor-finish(x)(if(Math-primpx)x(if(eq(carx)'calcFunc-Fac-Prot)(math-factor-finish(nth1x))(cons(carx)(mapcar'math-factor-finish(cdrx))))))(defunmath-factor-protect(x)(if(memq(car-safex)'(+-))(list'calcFunc-Fac-Protx)x))(defunmath-factor-expr(expr)(cond((eqmath-factored-varst)expr)((or(memq(car-safeexpr)'(*/^neg))(assq(car-safeexpr)calc-tweak-eqn-table))(cons(carexpr)(mapcar'math-factor-expr(cdrexpr))))((memq(car-safeexpr)'(+-))(let*((math-factored-varsmath-factored-vars)(y(catch'factor(math-factor-expr-partexpr))))(ify(math-factor-expry)expr)))(texpr)))(defunmath-factor-expr-part(x); uses "expr"(if(memq(car-safex)'(+-*/^neg))(while(setqx(cdrx))(math-factor-expr-part(carx)))(and(not(Math-objvecpx))(not(assocxmath-factored-vars))(>(math-factor-containsexprx)1)(setqmath-factored-vars(cons(listx)math-factored-vars))(math-factor-expr-tryx))))(defunmath-factor-expr-try(x)(if(eq(car-safeexpr)'*)(let((res1(catch'factor(let((expr(nth1expr)))(math-factor-expr-tryx))))(res2(catch'factor(let((expr(nth2expr)))(math-factor-expr-tryx)))))(and(orres1res2)(throw'factor(math-accum-factors(orres1(nth1expr))1(orres2(nth2expr))))))(let*((p(math-is-polynomialexprx30'gen))(math-poly-modulus(math-poly-modulusexpr))res)(and(cdrp)(setqres(math-factor-poly-coefsp))(throw'factorres)))))(defunmath-accum-factors(facpowfacs)(ifmath-to-list(if(math-vectorpfac)(progn(while(setqfac(cdrfac))(setqfacs(math-accum-factors(nth1(carfac))(*pow(nth2(carfac)))facs)))facs)(if(and(eq(car-safefac)'^)(natnump(nth2fac)))(setqpow(*pow(nth2fac))fac(nth1fac)))(if(eqfac1)facs(or(math-vectorpfacs)(setqfacs(if(eqfacs1)'(vec)(list'vec(list'vecfacs1)))))(let((foundfacs))(while(and(setqfound(cdrfound))(not(equalfac(nth1(carfound))))))(iffound(progn(setcar(cdr(cdr(carfound)))(+pow(nth2(carfound))))facs);; Put constant term first.(if(and(cdrfacs)(Math-ratp(nth1(nth1facs))))(cons'vec(cons(nth1facs)(cons(list'vecfacpow)(cdr(cdrfacs)))))(cons'vec(cons(list'vecfacpow)(cdrfacs))))))))(math-mul(math-powfacpow)facs)))(defunmath-factor-poly-coefs(p&optionalsquare-free); uses "x"(let(t1t2)(cond((not(cdrp))(or(carp)0));; Strip off multiples of x.((Math-zerop(carp))(let((z0))(while(andp(Math-zerop(carp)))(setqz(1+z)p(cdrp)))(if(cdrp)(setqp(math-factor-poly-coefspsquare-free))(setqp(math-sort-terms(math-factor-expr(carp)))))(math-accum-factorsxz(math-factor-protectp))));; Factor out content.((and(notsquare-free)(not(eq1(setqt1(math-mul(math-poly-gcd-listp)(if(math-guess-if-neg(nth(1-(lengthp))p))-11))))))(math-accum-factorst11(math-factor-poly-coefs(math-poly-div-listpt1)'cont)));; Check if linear in x.((not(cdr(cdrp)))(math-add(math-factor-protect(math-sort-terms(math-factor-expr(carp))))(math-mulx(math-factor-protect(math-sort-terms(math-factor-expr(nth1p)))))));; If symbolic coefficients, use FactorRules.((let((ppp))(while(andpp(or(Math-ratp(carpp))(and(eq(car(carpp))'mod)(Math-integerp(nth1(carpp)))(Math-integerp(nth2(carpp))))))(setqpp(cdrpp)))pp)(let((res(math-rewrite(list'calcFunc-thecoefsx(cons'vecp))'(varFactorRulesvar-FactorRules))))(or(and(eq(car-saferes)'calcFunc-thefactors)(=(lengthres)3)(math-vectorp(nth2res))(let((facs1)(vec(nth2res)))(while(setqvec(cdrvec))(setqfacs(math-accum-factors(carvec)1facs)))facs))(math-build-polynomial-exprpx))));; Check if rational coefficients (i.e., not modulo a prime).((eqmath-poly-modulus1);; Check if there are any squared terms, or a content not = 1.(if(or(eqsquare-freet)(equal(setqt1(math-poly-gcd-coefsp(setqt2(math-poly-deriv-coefsp))))'(1)));; We now have a square-free polynomial with integer coefs.;; For now, we use a kludgey method that finds linear and;; quadratic terms using floating-point root-finding.(if(setqt1(let((calc-symbolic-modenil))(math-poly-all-rootsnilpt)))(let((roots(cart1))(csign(if(math-negp(nth(1-(lengthp))p))-11))(expr1)(unfac(nth1t1))(scale(nth2t1)))(whileroots(let((coef0(car(carroots)))(coef1(cdr(carroots))))(setqexpr(math-accum-factors(ifcoef1(let((den(math-lcm-denomscoef0coef1)))(setqscale(math-divscaleden))(math-add(math-add(math-mulden(math-powx2))(math-mul(math-mulcoef1den)x))(math-mulcoef0den)))(let((den(math-lcm-denomscoef0)))(setqscale(math-divscaleden))(math-add(math-muldenx)(math-mulcoef0den))))1expr)roots(cdrroots))))(setqexpr(math-accum-factorsexpr1(math-mulcsign(math-build-polynomial-expr(math-mul-list(nth1t1)scale)x)))))(math-build-polynomial-exprpx)); can't factor it.;; Separate out the squared terms (Knuth exercise 4.6.2-34).;; This step also divides out the content of the polynomial.(let*((cabs(math-poly-gcd-listp))(csign(if(math-negp(nth(1-(lengthp))p))-11))(t1s(math-mul-listt1csign))(uunil)(v(car(math-poly-div-coefspt1s)))(w(car(math-poly-div-coefst2t1s))))(while(not(math-poly-zerop(setqt2(math-poly-simplify(math-poly-mixw1(math-poly-deriv-coefsv)-1)))))(setqt1(math-poly-gcd-coefsvt2)uu(const1uu)v(car(math-poly-div-coefsvt1))w(car(math-poly-div-coefst2t1))))(setqt1(lengthuu)t2(math-accum-factors(math-factor-poly-coefsvt)(1+t1)1))(whileuu(setqt2(math-accum-factors(math-factor-poly-coefs(caruu)t)t1t2)t1(1-t1)uu(cdruu)))(math-accum-factors(math-mulcabscsign)1t2))));; Factoring modulo a prime.((and(=(length(setqtemp(math-poly-gcd-coefsp(math-poly-deriv-coefsp))))(lengthp)))(setqp(cartemp))(while(cdrtemp)(setqtemp(nthcdr(nth2math-poly-modulus)temp)p(cons(cartemp)p)))(and(setqtemp(math-factor-poly-coefsp))(math-powtemp(nth2math-poly-modulus))))(t(math-reject-argnil"*Modulo factorization not yet implemented")))))(defunmath-poly-deriv-coefs(p)(let((n1)(dpnil))(while(setqp(cdrp))(setqdp(cons(math-mul(carp)n)dp)n(1+n)))(nreversedp)))(defunmath-factor-contains(xa)(if(equalxa)1(if(memq(car-safex)'(+-*/neg))(let((sum0))(while(setqx(cdrx))(setqsum(+sum(math-factor-contains(carx)a))))sum)(if(and(eq(car-safex)'^)(natnump(nth2x)))(*(math-factor-contains(nth1x)a)(nth2x))0))));;; Merge all quotients and expand/simplify the numerator(defuncalcFunc-nrat(expr)(if(math-any-floatsexpr)(setqexpr(calcFunc-pfracexpr)))(if(or(math-vectorpexpr)(assq(car-safeexpr)calc-tweak-eqn-table))(cons(carexpr)(mapcar'calcFunc-nrat(cdrexpr)))(let*((calc-prefer-fract)(res(math-to-ratpolyexpr))(num(math-simplify(math-sort-terms(calcFunc-expand(carres)))))(den(math-simplify(math-sort-terms(calcFunc-expand(cdrres)))))(g(math-poly-gcdnumden)))(or(eqg1)(let((num2(math-poly-divnumg))(den2(math-poly-divdeng)))(and(eq(cdrnum2)0)(eq(cdrden2)0)(setqnum(carnum2)den(carden2)))))(math-simplify(math-divnumden)))));;; Returns expressions (num . denom).(defunmath-to-ratpoly(expr)(let((res(math-to-ratpoly-recexpr)))(cons(math-simplify(carres))(math-simplify(cdrres)))))(defunmath-to-ratpoly-rec(expr)(cond((Math-primpexpr)(consexpr1))((memq(carexpr)'(+-))(let((r1(math-to-ratpoly-rec(nth1expr)))(r2(math-to-ratpoly-rec(nth2expr))))(if(equal(cdrr1)(cdrr2))(cons(list(carexpr)(carr1)(carr2))(cdrr1))(if(eq(cdrr1)1)(cons(list(carexpr)(math-mul(carr1)(cdrr2))(carr2))(cdrr2))(if(eq(cdrr2)1)(cons(list(carexpr)(carr1)(math-mul(carr2)(cdrr1)))(cdrr1))(let((g(math-poly-gcd(cdrr1)(cdrr2))))(let((d1(and(not(eqg1))(math-poly-div(cdrr1)g)))(d2(and(not(eqg1))(math-poly-div(math-mul(carr1)(cdrr2))g))))(if(and(eq(cdrd1)0)(eq(cdrd2)0))(cons(list(carexpr)(card2)(math-mul(carr2)(card1)))(math-mul(card1)(cdrr2)))(cons(list(carexpr)(math-mul(carr1)(cdrr2))(math-mul(carr2)(cdrr1)))(math-mul(cdrr1)(cdrr2)))))))))))((eq(carexpr)'*)(let*((r1(math-to-ratpoly-rec(nth1expr)))(r2(math-to-ratpoly-rec(nth2expr)))(g(math-mul(math-poly-gcd(carr1)(cdrr2))(math-poly-gcd(cdrr1)(carr2)))))(if(eqg1)(cons(math-mul(carr1)(carr2))(math-mul(cdrr1)(cdrr2)))(cons(math-poly-div-exact(math-mul(carr1)(carr2))g)(math-poly-div-exact(math-mul(cdrr1)(cdrr2))g)))))((eq(carexpr)'/)(let*((r1(math-to-ratpoly-rec(nth1expr)))(r2(math-to-ratpoly-rec(nth2expr))))(if(and(eq(cdrr1)1)(eq(cdrr2)1))(cons(carr1)(carr2))(let((g(math-mul(math-poly-gcd(carr1)(carr2))(math-poly-gcd(cdrr1)(cdrr2)))))(if(eqg1)(cons(math-mul(carr1)(cdrr2))(math-mul(cdrr1)(carr2)))(cons(math-poly-div-exact(math-mul(carr1)(cdrr2))g)(math-poly-div-exact(math-mul(cdrr1)(carr2))g)))))))((and(eq(carexpr)'^)(integerp(nth2expr)))(let((r1(math-to-ratpoly-rec(nth1expr))))(if(>(nth2expr)0)(cons(math-pow(carr1)(nth2expr))(math-pow(cdrr1)(nth2expr)))(cons(math-pow(cdrr1)(-(nth2expr)))(math-pow(carr1)(-(nth2expr)))))))((eq(carexpr)'neg)(let((r1(math-to-ratpoly-rec(nth1expr))))(cons(math-neg(carr1))(cdrr1))))(t(consexpr1))))(defunmath-ratpoly-p(expr&optionalvar)(cond((equalexprvar)1)((Math-primpexpr)0)((memq(carexpr)'(+-))(let((p1(math-ratpoly-p(nth1expr)var))p2)(andp1(setqp2(math-ratpoly-p(nth2expr)var))(maxp1p2))))((eq(carexpr)'*)(let((p1(math-ratpoly-p(nth1expr)var))p2)(andp1(setqp2(math-ratpoly-p(nth2expr)var))(+p1p2))))((eq(carexpr)'neg)(math-ratpoly-p(nth1expr)var))((eq(carexpr)'/)(let((p1(math-ratpoly-p(nth1expr)var))p2)(andp1(setqp2(math-ratpoly-p(nth2expr)var))(-p1p2))))((and(eq(carexpr)'^)(integerp(nth2expr)))(let((p1(math-ratpoly-p(nth1expr)var)))(andp1(*p1(nth2expr)))))((notvar)1)((math-poly-dependsexprvar)nil)(t0)))(defuncalcFunc-apart(expr&optionalvar)(cond((Math-primpexpr)expr)((eq(carexpr)'+)(math-add(calcFunc-apart(nth1expr)var)(calcFunc-apart(nth2expr)var)))((eq(carexpr)'-)(math-sub(calcFunc-apart(nth1expr)var)(calcFunc-apart(nth2expr)var)))((not(math-ratpoly-pexprvar))(math-reject-argexpr"Expected a rational function"))(t(let*((calc-prefer-fract)(rat(math-to-ratpolyexpr))(num(carrat))(den(cdrrat))(qr(math-poly-divnumden))(q(carqr))(r(cdrqr)))(orvar(setqvar(math-polynomial-baseden)))(math-addq(or(andvar(math-expr-containsdenvar)(math-partial-fractionsrdenvar))(math-divrden)))))))(defunmath-padded-polynomial(exprvardeg)(let((p(math-is-polynomialexprvardeg)))(appendp(make-list(-deg(lengthp))0))))(defunmath-partial-fractions(rdenvar)(let*((fden(calcFunc-factorsdenvar))(tdeg(math-polynomial-pdenvar))(fpfden)(dlistnil)(eqns0)(lznil)(tz(make-list(1-tdeg)0))(calc-matrix-mode'scalar))(and(not(and(=(lengthfden)2)(eq(nth2(nth1fden))1)))(progn(while(setqfp(cdrfp))(let((rpt(nth2(carfp)))(deg(math-polynomial-p(nth1(carfp))var))dnumdvardeg2)(while(>rpt0)(setqdeg2degdnum0)(while(>deg20)(setqdvar(append'(vec)lz'(1)tz)lz(cons0lz)tz(cdrtz)deg2(1-deg2)dnum(math-adddnum(math-muldvar(math-powvardeg2)))dlist(cons(and(=deg2(1-deg))(math-pow(nth1(carfp))rpt))dlist)))(let((fppfden)(mult1))(while(setqfpp(cdrfpp))(or(eqfppfp)(setqmult(math-mulmult(math-pow(nth1(carfpp))(nth2(carfpp)))))))(setqdnum(math-muldnummult)))(setqeqns(math-addeqns(math-muldnum(math-pow(nth1(carfp))(-(nth2(carfp))rpt))))rpt(1-rpt)))))(setqeqns(math-div(cons'vec(math-padded-polynomialrvartdeg))(math-transpose(cons'vec(mapcar(function(lambda(x)(cons'vec(math-padded-polynomialxvartdeg))))(cdreqns))))))(and(math-vectorpeqns)(let((res0)(numnil))(setqeqns(nreverseeqns))(whileeqns(setqnum(cons(careqns)num)eqns(cdreqns))(if(cardlist)(setqnum(math-build-polynomial-expr(nreversenum)var)res(math-addres(math-divnum(cardlist)))numnil))(setqdlist(cdrdlist)))(math-normalizeres)))))))(defunmath-expand-term(expr)(cond((and(eq(car-safeexpr)'*)(memq(car-safe(nth1expr))'(+-)))(math-add-or-sub(list'*(nth1(nth1expr))(nth2expr))(list'*(nth2(nth1expr))(nth2expr))nil(eq(car(nth1expr))'-)))((and(eq(car-safeexpr)'*)(memq(car-safe(nth2expr))'(+-)))(math-add-or-sub(list'*(nth1expr)(nth1(nth2expr)))(list'*(nth1expr)(nth2(nth2expr)))nil(eq(car(nth2expr))'-)))((and(eq(car-safeexpr)'/)(memq(car-safe(nth1expr))'(+-)))(math-add-or-sub(list'/(nth1(nth1expr))(nth2expr))(list'/(nth2(nth1expr))(nth2expr))nil(eq(car(nth1expr))'-)))((and(eq(car-safeexpr)'^)(memq(car-safe(nth1expr))'(+-))(integerp(nth2expr))(if(>(nth2expr)0)(or(and(or(>mmt-many500000)(<mmt-many-500000))(math-expand-power(nth1expr)(nth2expr)nilt))(list'*(nth1expr)(list'^(nth1expr)(1-(nth2expr)))))(if(<(nth2expr)0)(list'/1(list'^(nth1expr)(-(nth2expr))))))))(texpr)))(defuncalcFunc-expand(expr&optionalmany)(math-normalize(math-map-tree'math-expand-termexprmany)))(defunmath-expand-power(xn&optionalvarelse-nil)(or(and(natnumpn)(memq(car-safex)'(+-))(let((termsnil)(ctermsnil))(while(memq(car-safex)'(+-))(setqterms(cons(if(eq(carx)'-)(math-neg(nth2x))(nth2x))terms)x(nth1x)))(setqterms(consxterms))(ifvar(let((pterms))(whilep(or(math-expr-contains(carp)var)(setqterms(delq(carp)terms)cterms(cons(carp)cterms)))(setqp(cdrp)))(ifcterms(setqterms(cons(apply'calcFunc-addcterms)terms)))))(if(=(lengthterms)2)(let((i0)(accum0))(while(<=in)(setqaccum(list'+accum(list'*(calcFunc-chooseni)(list'*(list'^(nth1terms)i)(list'^(carterms)(-ni)))))i(1+i)))accum)(if(=n2)(let((accum0)(p1terms)p2)(whilep1(setqaccum(list'+accum(list'^(carp1)2))p2p1)(while(setqp2(cdrp2))(setqaccum(list'+accum(list'*2(list'*(carp1)(carp2))))))(setqp1(cdrp1)))accum)(if(=n3)(let((accum0)(p1terms)p2p3)(whilep1(setqaccum(list'+accum(list'^(carp1)3))p2p1)(while(setqp2(cdrp2))(setqaccum(list'+(list'+accum(list'*3(list'*(list'^(carp1)2)(carp2))))(list'*3(list'*(carp1)(list'^(carp2)2))))p3p2)(while(setqp3(cdrp3))(setqaccum(list'+accum(list'*6(list'*(carp1)(list'*(carp2)(carp3))))))))(setqp1(cdrp1)))accum))))))(and(notelse-nil)(list'^xn))))(defuncalcFunc-expandpow(xn)(math-normalize(math-expand-powerxn)))