\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra intrf.spad} \author{Barry Trager, Renaud Rioboo, Manuel Bronstein} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package SUBRESP SubResultantPackage} <>= )abbrev package SUBRESP SubResultantPackage ++ Subresultants ++ Author: Barry Trager, Renaud Rioboo ++ Date Created: 1987 ++ Date Last Updated: August 2000 ++ Description: ++ This package computes the subresultants of two polynomials which is needed ++ for the `Lazard Rioboo' enhancement to Tragers integrations formula ++ For efficiency reasons this has been rewritten to call Lionel Ducos ++ package which is currently the best one. ++ SubResultantPackage(R, UP): Exports == Implementation where R : IntegralDomain UP: UnivariatePolynomialCategory R Z ==> Integer N ==> NonNegativeInteger Exports ==> with subresultantVector: (UP, UP) -> PrimitiveArray UP ++ subresultantVector(p, q) returns \spad{[p0,...,pn]} ++ where pi is the i-th subresultant of p and q. ++ In particular, \spad{p0 = resultant(p, q)}. if R has EuclideanDomain then primitivePart : (UP, R) -> UP ++ primitivePart(p, q) reduces the coefficient of p ++ modulo q, takes the primitive part of the result, ++ and ensures that the leading coefficient of that ++ result is monic. Implementation ==> add Lionel ==> PseudoRemainderSequence(R,UP) if R has EuclideanDomain then primitivePart(p, q) == rec := extendedEuclidean(leadingCoefficient p, q, 1)::Record(coef1:R, coef2:R) unitCanonical primitivePart map((rec.coef1 * #1) rem q, p) subresultantVector(p1, p2) == F : UP -- auxiliary stuff ! res : PrimitiveArray(UP) := new(2+max(degree(p1),degree(p2)), 0) -- -- kind of stupid interface to Lionel's Package !!!!!!!!!!!! -- might have been wiser to rewrite the loop ... -- But I'm too lazy. [rr] -- l := chainSubResultants(p1,p2)$Lionel -- -- this returns the chain of non null subresultants ! -- we must rebuild subresultants from this. -- we really hope Lionel Ducos minded what he wrote -- since we are fully blind ! -- null l => -- Hum it seems that Lionel returns [] when min(|p1|,|p2|) = 0 zero?(degree(p1)) => res.degree(p2) := p2 if degree(p2) > 0 then res.((degree(p2)-1)::NonNegativeInteger) := p1 res.0 := (leadingCoefficient(p1)**(degree p2)) :: UP else -- both are of degree 0 the resultant is 1 according to Loos res.0 := 1 res zero?(degree(p2)) => if degree(p1) > 0 then res.((degree(p1)-1)::NonNegativeInteger) := p2 res.0 := (leadingCoefficient(p2)**(degree p1)) :: UP else -- both are of degree 0 the resultant is 1 according to Loos res.0 := 1 res error "SUBRESP: strange Subresultant chain from PRS" Sn := first(l) -- -- as of Loos definitions last subresultant should not be defective -- l := rest(l) n := degree(Sn) F := Sn null l => error "SUBRESP: strange Subresultant chain from PRS" zero? Sn => error "SUBRESP: strange Subresultant chain from PRS" while (l ~= []) repeat res.(n) := Sn F := first(l) l := rest(l) -- F is potentially defective if degree(F) = n then -- -- F is defective -- null l => error "SUBRESP: strange Subresultant chain from PRS" Sn := first(l) l := rest(l) n := degree(Sn) res.((n-1)::NonNegativeInteger) := F else -- -- F is non defective -- degree(F) < n => error "strange result !" Sn := F n := degree(Sn) -- -- Lionel forgets about p1 if |p1| > |p2| -- forgets about p2 if |p2| > |p1| -- but he reminds p2 if |p1| = |p2| -- a glance at Loos should correct this ! -- res.n := Sn -- -- Loos definition -- if degree(p1) = degree(p2) then res.((degree p1)+1) := p1 else if degree(p1) > degree(p2) then res.(degree p1) := p1 else res.(degree p2) := p2 res @ \section{package MONOTOOL MonomialExtensionTools} <>= )abbrev package MONOTOOL MonomialExtensionTools ++ Tools for handling monomial extensions ++ Author: Manuel Bronstein ++ Date Created: 18 August 1992 ++ Date Last Updated: 3 June 1993 ++ Description: Tools for handling monomial extensions. MonomialExtensionTools(F, UP): Exports == Implementation where F : Field UP: UnivariatePolynomialCategory F RF ==> Fraction UP FR ==> Factored UP Exports ==> with split : (UP, UP -> UP) -> Record(normal:UP, special:UP) ++ split(p, D) returns \spad{[n,s]} such that \spad{p = n s}, ++ all the squarefree factors of n are normal w.r.t. D, ++ and s is special w.r.t. D. ++ D is the derivation to use. splitSquarefree: (UP, UP -> UP) -> Record(normal:FR, special:FR) ++ splitSquarefree(p, D) returns ++ \spad{[n_1 n_2\^2 ... n_m\^m, s_1 s_2\^2 ... s_q\^q]} such that ++ \spad{p = n_1 n_2\^2 ... n_m\^m s_1 s_2\^2 ... s_q\^q}, each ++ \spad{n_i} is normal w.r.t. D and each \spad{s_i} is special ++ w.r.t D. ++ D is the derivation to use. normalDenom: (RF, UP -> UP) -> UP ++ normalDenom(f, D) returns the product of all the normal factors ++ of \spad{denom(f)}. ++ D is the derivation to use. decompose : (RF, UP -> UP) -> Record(poly:UP, normal:RF, special:RF) ++ decompose(f, D) returns \spad{[p,n,s]} such that \spad{f = p+n+s}, ++ all the squarefree factors of \spad{denom(n)} are normal w.r.t. D, ++ \spad{denom(s)} is special w.r.t. D, ++ and n and s are proper fractions (no pole at infinity). ++ D is the derivation to use. Implementation ==> add normalDenom(f, derivation) == split(denom f, derivation).normal split(p, derivation) == pbar := (gcd(p, derivation p) exquo gcd(p, differentiate p))::UP zero? degree pbar => [p, 1] rec := split((p exquo pbar)::UP, derivation) [rec.normal, pbar * rec.special] splitSquarefree(p, derivation) == s:Factored(UP) := 1 n := s q := squareFree p for rec in factors q repeat r := rec.factor g := gcd(r, derivation r) if not ground? g then s := s * sqfrFactor(g, rec.exponent) h := (r exquo g)::UP if not ground? h then n := n * sqfrFactor(h, rec.exponent) [n, unit(q) * s] decompose(f, derivation) == qr := divide(numer f, denom f) -- rec.normal * rec.special = denom f rec := split(denom f, derivation) -- eeu.coef1 * rec.normal + eeu.coef2 * rec.special = qr.remainder -- and degree(eeu.coef1) < degree(rec.special) -- and degree(eeu.coef2) < degree(rec.normal) -- qr.remainder/denom(f) = eeu.coef1 / rec.special + eeu.coef2 / rec.normal eeu := extendedEuclidean(rec.normal, rec.special, qr.remainder)::Record(coef1:UP, coef2:UP) [qr.quotient, eeu.coef2 / rec.normal, eeu.coef1 / rec.special] @ \section{package INTHERTR TranscendentalHermiteIntegration} <>= )abbrev package INTHERTR TranscendentalHermiteIntegration ++ Hermite integration, transcendental case ++ Author: Manuel Bronstein ++ Date Created: 1987 ++ Date Last Updated: 12 August 1992 ++ Description: Hermite integration, transcendental case. TranscendentalHermiteIntegration(F, UP): Exports == Implementation where F : Field UP : UnivariatePolynomialCategory F N ==> NonNegativeInteger RF ==> Fraction UP REC ==> Record(answer:RF, lognum:UP, logden:UP) HER ==> Record(answer:RF, logpart:RF, specpart:RF, polypart:UP) Exports ==> with HermiteIntegrate: (RF, UP -> UP) -> HER ++ HermiteIntegrate(f, D) returns \spad{[g, h, s, p]} ++ such that \spad{f = Dg + h + s + p}, ++ h has a squarefree denominator normal w.r.t. D, ++ and all the squarefree factors of the denominator of s are ++ special w.r.t. D. Furthermore, h and s have no polynomial parts. ++ D is the derivation to use on \spadtype{UP}. Implementation ==> add import MonomialExtensionTools(F, UP) normalHermiteIntegrate: (RF,UP->UP) -> Record(answer:RF,lognum:UP,logden:UP) HermiteIntegrate(f, derivation) == rec := decompose(f, derivation) hi := normalHermiteIntegrate(rec.normal, derivation) qr := divide(hi.lognum, hi.logden) [hi.answer, qr.remainder / hi.logden, rec.special, qr.quotient + rec.poly] -- Hermite Reduction on f, every squarefree factor of denom(f) is normal wrt D -- this is really a "parallel" Hermite reduction, in the sense that -- every multiple factor of the denominator gets reduced at each pass -- so if the denominator is P1 P2**2 ... Pn**n, this requires O(n) -- reduction steps instead of O(n**2), like Mack's algorithm -- (D.Mack, On Rational Integration, Univ. of Utah C.S. Tech.Rep. UCP-38,1975) -- returns [g, b, d] s.t. f = g' + b/d and d is squarefree and normal wrt D normalHermiteIntegrate(f, derivation) == a := numer f q := denom f p:UP := 0 mult:UP := 1 qhat := (q exquo (g0 := g := gcd(q, differentiate q)))::UP while(degree(qbar := g) > 0) repeat qbarhat := (qbar exquo (g := gcd(qbar, differentiate qbar)))::UP qtil:= - ((qhat * (derivation qbar)) exquo qbar)::UP bc := extendedEuclidean(qtil, qbarhat, a)::Record(coef1:UP, coef2:UP) qr := divide(bc.coef1, qbarhat) a := bc.coef2 + qtil * qr.quotient - derivation(qr.remainder) * (qhat exquo qbarhat)::UP p := p + mult * qr.remainder mult:= mult * qbarhat [p / g0, a, qhat] @ \section{package INTTR TranscendentalIntegration} <>= )abbrev package INTTR TranscendentalIntegration ++ Risch algorithm, transcendental case ++ Author: Manuel Bronstein ++ Date Created: 1987 ++ Date Last Updated: 24 October 1995 ++ Description: ++ This package provides functions for the transcendental ++ case of the Risch algorithm. -- Internally used by the integrator TranscendentalIntegration(F, UP): Exports == Implementation where F : Field UP : UnivariatePolynomialCategory F N ==> NonNegativeInteger Z ==> Integer Q ==> Fraction Z GP ==> LaurentPolynomial(F, UP) UP2 ==> SparseUnivariatePolynomial UP RF ==> Fraction UP UPR ==> SparseUnivariatePolynomial RF IR ==> IntegrationResult RF LOG ==> Record(scalar:Q, coeff:UPR, logand:UPR) LLG ==> List Record(coeff:RF, logand:RF) NE ==> Record(integrand:RF, intvar:RF) NL ==> Record(mainpart:RF, limitedlogs:LLG) UPF ==> Record(answer:UP, a0:F) RFF ==> Record(answer:RF, a0:F) IRF ==> Record(answer:IR, a0:F) NLF ==> Record(answer:NL, a0:F) GPF ==> Record(answer:GP, a0:F) UPUP==> Record(elem:UP, notelem:UP) GPGP==> Record(elem:GP, notelem:GP) RFRF==> Record(elem:RF, notelem:RF) FF ==> Record(ratpart:F, coeff:F) FFR ==> Record(ratpart:RF, coeff:RF) UF ==> Union(FF, "failed") UF2 ==> Union(List F, "failed") REC ==> Record(ir:IR, specpart:RF, polypart:UP) PSOL==> Record(ans:F, right:F, sol?:Boolean) FAIL==> error "Sorry - cannot handle that integrand yet" Exports ==> with primintegrate : (RF, UP -> UP, F -> UF) -> IRF ++ primintegrate(f, ', foo) returns \spad{[g, a]} such that ++ \spad{f = g' + a}, and \spad{a = 0} or \spad{a} has no integral in UP. ++ Argument foo is an extended integration function on F. expintegrate : (RF, UP -> UP, (Z, F) -> PSOL) -> IRF ++ expintegrate(f, ', foo) returns \spad{[g, a]} such that ++ \spad{f = g' + a}, and \spad{a = 0} or \spad{a} has no integral in F; ++ Argument foo is a Risch differential equation solver on F; tanintegrate : (RF, UP -> UP, (Z, F, F) -> UF2) -> IRF ++ tanintegrate(f, ', foo) returns \spad{[g, a]} such that ++ \spad{f = g' + a}, and \spad{a = 0} or \spad{a} has no integral in F; ++ Argument foo is a Risch differential system solver on F; primextendedint:(RF, UP -> UP, F->UF, RF) -> Union(RFF,FFR,"failed") ++ primextendedint(f, ', foo, g) returns either \spad{[v, c]} such that ++ \spad{f = v' + c g} and \spad{c' = 0}, or \spad{[v, a]} such that ++ \spad{f = g' + a}, and \spad{a = 0} or \spad{a} has no integral in UP. ++ Returns "failed" if neither case can hold. ++ Argument foo is an extended integration function on F. expextendedint:(RF,UP->UP,(Z,F)->PSOL, RF) -> Union(RFF,FFR,"failed") ++ expextendedint(f, ', foo, g) returns either \spad{[v, c]} such that ++ \spad{f = v' + c g} and \spad{c' = 0}, or \spad{[v, a]} such that ++ \spad{f = g' + a}, and \spad{a = 0} or \spad{a} has no integral in F. ++ Returns "failed" if neither case can hold. ++ Argument foo is a Risch differential equation function on F. primlimitedint:(RF, UP -> UP, F->UF, List RF) -> Union(NLF,"failed") ++ primlimitedint(f, ', foo, [u1,...,un]) returns ++ \spad{[v, [c1,...,cn], a]} such that \spad{ci' = 0}, ++ \spad{f = v' + a + reduce(+,[ci * ui'/ui])}, ++ and \spad{a = 0} or \spad{a} has no integral in UP. ++ Returns "failed" if no such v, ci, a exist. ++ Argument foo is an extended integration function on F. explimitedint:(RF, UP->UP,(Z,F)->PSOL,List RF) -> Union(NLF,"failed") ++ explimitedint(f, ', foo, [u1,...,un]) returns ++ \spad{[v, [c1,...,cn], a]} such that \spad{ci' = 0}, ++ \spad{f = v' + a + reduce(+,[ci * ui'/ui])}, ++ and \spad{a = 0} or \spad{a} has no integral in F. ++ Returns "failed" if no such v, ci, a exist. ++ Argument foo is a Risch differential equation function on F. primextintfrac : (RF, UP -> UP, RF) -> Union(FFR, "failed") ++ primextintfrac(f, ', g) returns \spad{[v, c]} such that ++ \spad{f = v' + c g} and \spad{c' = 0}. ++ Error: if \spad{degree numer f >= degree denom f} or ++ if \spad{degree numer g >= degree denom g} or ++ if \spad{denom g} is not squarefree. primlimintfrac : (RF, UP -> UP, List RF) -> Union(NL, "failed") ++ primlimintfrac(f, ', [u1,...,un]) returns \spad{[v, [c1,...,cn]]} ++ such that \spad{ci' = 0} and \spad{f = v' + +/[ci * ui'/ui]}. ++ Error: if \spad{degree numer f >= degree denom f}. primintfldpoly : (UP, F -> UF, F) -> Union(UP, "failed") ++ primintfldpoly(p, ', t') returns q such that \spad{p' = q} or ++ "failed" if no such q exists. Argument \spad{t'} is the derivative of ++ the primitive generating the extension. expintfldpoly : (GP, (Z, F) -> PSOL) -> Union(GP, "failed") ++ expintfldpoly(p, foo) returns q such that \spad{p' = q} or ++ "failed" if no such q exists. ++ Argument foo is a Risch differential equation function on F. monomialIntegrate : (RF, UP -> UP) -> REC ++ monomialIntegrate(f, ') returns \spad{[ir, s, p]} such that ++ \spad{f = ir' + s + p} and all the squarefree factors of the ++ denominator of s are special w.r.t the derivation '. monomialIntPoly : (UP, UP -> UP) -> Record(answer:UP, polypart:UP) ++ monomialIntPoly(p, ') returns [q, r] such that ++ \spad{p = q' + r} and \spad{degree(r) < degree(t')}. ++ Error if \spad{degree(t') < 2}. Implementation ==> add import SubResultantPackage(UP, UP2) import MonomialExtensionTools(F, UP) import TranscendentalHermiteIntegration(F, UP) import CommuteUnivariatePolynomialCategory(F, UP, UP2) primintegratepoly : (UP, F -> UF, F) -> Union(UPF, UPUP) expintegratepoly : (GP, (Z, F) -> PSOL) -> Union(GPF, GPGP) expextintfrac : (RF, UP -> UP, RF) -> Union(FFR, "failed") explimintfrac : (RF, UP -> UP, List RF) -> Union(NL, "failed") limitedLogs : (RF, RF -> RF, List RF) -> Union(LLG, "failed") logprmderiv : (RF, UP -> UP) -> RF logexpderiv : (RF, UP -> UP, F) -> RF tanintegratespecial: (RF, RF -> RF, (Z, F, F) -> UF2) -> Union(RFF, RFRF) UP2UP2 : UP -> UP2 UP2UPR : UP -> UPR UP22UPR : UP2 -> UPR notelementary : REC -> IR kappa : (UP, UP -> UP) -> UP dummy:RF := 0 logprmderiv(f, derivation) == differentiate(f, derivation) / f UP2UP2 p == map(#1::UP, p)$UnivariatePolynomialCategoryFunctions2(F, UP, UP, UP2) UP2UPR p == map(#1::UP::RF, p)$UnivariatePolynomialCategoryFunctions2(F, UP, RF, UPR) UP22UPR p == map(#1::RF, p)$SparseUnivariatePolynomialFunctions2(UP, RF) -- given p in k[z] and a derivation on k[t] returns the coefficient lifting -- in k[z] of the restriction of D to k. kappa(p, derivation) == ans:UP := 0 while p ~= 0 repeat ans := ans + derivation(leadingCoefficient(p)::UP)*monomial(1,degree p) p := reductum p ans -- works in any monomial extension monomialIntegrate(f, derivation) == zero? f => [0, 0, 0] r := HermiteIntegrate(f, derivation) zero?(inum := numer(r.logpart)) => [r.answer::IR, r.specpart, r.polypart] iden := denom(r.logpart) x := monomial(1, 1)$UP resultvec := subresultantVector(UP2UP2 inum - (x::UP2) * UP2UP2 derivation iden, UP2UP2 iden) respoly := primitivePart leadingCoefficient resultvec 0 rec := splitSquarefree(respoly, kappa(#1, derivation)) logs:List(LOG) := [ [1, UP2UPR(term.factor), UP22UPR swap primitivePart(resultvec(term.exponent),term.factor)] for term in factors(rec.special)] dlog := one? derivation x => r.logpart differentiate(mkAnswer(0, logs, empty()), differentiate(#1, derivation)) (u := retractIfCan(p := r.logpart - dlog)@Union(UP, "failed")) case UP => [mkAnswer(r.answer, logs, empty), r.specpart, r.polypart + u::UP] [mkAnswer(r.answer, logs, [[p, dummy]]), r.specpart, r.polypart] -- returns [q, r] such that p = q' + r and degree(r) < degree(dt) -- must have degree(derivation t) >= 2 monomialIntPoly(p, derivation) == (d := degree(dt := derivation monomial(1,1))::Z) < 2 => error "monomIntPoly: monomial must have degree 2 or more" l := leadingCoefficient dt ans:UP := 0 while (n := 1 + degree(p)::Z - d) > 0 repeat ans := ans + (term := monomial(leadingCoefficient(p) / (n * l), n::N)) p := p - derivation term -- degree(p) must drop here [ans, p] -- returns either -- (q in GP, a in F) st p = q' + a, and a=0 or a has no integral in F -- or (q in GP, r in GP) st p = q' + r, and r has no integral elem/UP expintegratepoly(p, FRDE) == coef0:F := 0 notelm := answr := 0$GP while p ~= 0 repeat ans1 := FRDE(n := degree p, a := leadingCoefficient p) answr := answr + monomial(ans1.ans, n) if ~ans1.sol? then -- Risch d.e. has no complete solution missing := a - ans1.right if zero? n then coef0 := missing else notelm := notelm + monomial(missing, n) p := reductum p zero? notelm => [answr, coef0] [answr, notelm] -- f is either 0 or of the form p(t)/(1 + t**2)**n -- returns either -- (q in RF, a in F) st f = q' + a, and a=0 or a has no integral in F -- or (q in RF, r in RF) st f = q' + r, and r has no integral elem/UP tanintegratespecial(f, derivation, FRDE) == ans:RF := 0 p := monomial(1, 2)$UP + 1 while (n := degree(denom f) quo 2) ~= 0 repeat r := numer(f) rem p a := coefficient(r, 1) b := coefficient(r, 0) (u := FRDE(n, a, b)) case "failed" => return [ans, f] l := u::List(F) term:RF := (monomial(first l, 1)$UP + second(l)::UP) / denom f ans := ans + term f := f - derivation term -- the order of the pole at 1+t^2 drops zero?(c0 := retract(retract(f)@UP)@F) or (u := FRDE(0, c0, 0)) case "failed" => [ans, c0] [ans + first(u::List(F))::UP::RF, 0::F] -- returns (v in RF, c in RF) s.t. f = v' + cg, and c' = 0, or "failed" -- g must have a squarefree denominator (always possible) -- g must have no polynomial part and no pole above t = 0 -- f must have no polynomial part and no pole above t = 0 expextintfrac(f, derivation, g) == zero? f => [0, 0] degree numer f >= degree denom f => error "Not a proper fraction" order(denom f,monomial(1,1)) ~= 0 => error "Not integral at t = 0" r := HermiteIntegrate(f, derivation) zero? g => r.logpart ~= 0 => "failed" [r.answer, 0] degree numer g >= degree denom g => error "Not a proper fraction" order(denom g,monomial(1,1)) ~= 0 => error "Not integral at t = 0" differentiate(c := r.logpart / g, derivation) ~= 0 => "failed" [r.answer, c] limitedLogs(f, logderiv, lu) == zero? f => empty() empty? lu => "failed" empty? rest lu => logderiv(c0 := f / logderiv(u0 := first lu)) ~= 0 => "failed" [[c0, u0]] num := numer f den := denom f l1:List Record(logand2:RF, contrib:UP) := [[u, numer v] for u in lu | one? denom(v := den * logderiv u)] rows := max(degree den, 1 + reduce(max, [degree(u.contrib) for u in l1], 0)$List(N)) m:Matrix(F) := zero(rows, cols := 1 + #l1) for i in 0..rows-1 repeat for pp in l1 for j in minColIndex m .. maxColIndex m - 1 repeat qsetelt!(m, i + minRowIndex m, j, coefficient(pp.contrib, i)) qsetelt!(m,i+minRowIndex m, maxColIndex m, coefficient(num, i)) m := rowEchelon m ans := empty()$LLG for i in minRowIndex m .. maxRowIndex m | qelt(m, i, maxColIndex m) ~= 0 repeat OK := false for pp in l1 for j in minColIndex m .. maxColIndex m - 1 while not OK repeat if qelt(m, i, j) ~= 0 then OK := true c := qelt(m, i, maxColIndex m) / qelt(m, i, j) logderiv(c0 := c::UP::RF) ~= 0 => return "failed" ans := concat([c0, pp.logand2], ans) not OK => return "failed" ans -- returns q in UP s.t. p = q', or "failed" primintfldpoly(p, extendedint, t') == (u := primintegratepoly(p, extendedint, t')) case UPUP => "failed" u.a0 ~= 0 => "failed" u.answer -- returns q in GP st p = q', or "failed" expintfldpoly(p, FRDE) == (u := expintegratepoly(p, FRDE)) case GPGP => "failed" u.a0 ~= 0 => "failed" u.answer -- returns (v in RF, c1...cn in RF, a in F) s.t. ci' = 0, -- and f = v' + a + +/[ci * ui'/ui] -- and a = 0 or a has no integral in UP primlimitedint(f, derivation, extendedint, lu) == qr := divide(numer f, denom f) (u1 := primlimintfrac(qr.remainder / (denom f), derivation, lu)) case "failed" => "failed" (u2 := primintegratepoly(qr.quotient, extendedint, retract derivation monomial(1, 1))) case UPUP => "failed" [[u1.mainpart + u2.answer::RF, u1.limitedlogs], u2.a0] -- returns (v in RF, c1...cn in RF, a in F) s.t. ci' = 0, -- and f = v' + a + +/[ci * ui'/ui] -- and a = 0 or a has no integral in F explimitedint(f, derivation, FRDE, lu) == qr := separate(f)$GP (u1 := explimintfrac(qr.fracPart,derivation, lu)) case "failed" => "failed" (u2 := expintegratepoly(qr.polyPart, FRDE)) case GPGP => "failed" [[u1.mainpart + convert(u2.answer)@RF, u1.limitedlogs], u2.a0] -- returns [v, c1...cn] s.t. f = v' + +/[ci * ui'/ui] -- f must have no polynomial part (degree numer f < degree denom f) primlimintfrac(f, derivation, lu) == zero? f => [0, empty()] degree numer f >= degree denom f => error "Not a proper fraction" r := HermiteIntegrate(f, derivation) zero?(r.logpart) => [r.answer, empty()] (u := limitedLogs(r.logpart, logprmderiv(#1, derivation), lu)) case "failed" => "failed" [r.answer, u::LLG] -- returns [v, c1...cn] s.t. f = v' + +/[ci * ui'/ui] -- f must have no polynomial part (degree numer f < degree denom f) -- f must be integral above t = 0 explimintfrac(f, derivation, lu) == zero? f => [0, empty()] degree numer f >= degree denom f => error "Not a proper fraction" order(denom f, monomial(1,1)) > 0 => error "Not integral at t = 0" r := HermiteIntegrate(f, derivation) zero?(r.logpart) => [r.answer, empty()] eta' := coefficient(derivation monomial(1, 1), 1) (u := limitedLogs(r.logpart, logexpderiv(#1,derivation,eta'), lu)) case "failed" => "failed" [r.answer - eta'::UP * +/[((degree numer(v.logand))::Z - (degree denom(v.logand))::Z) * v.coeff for v in u], u::LLG] logexpderiv(f, derivation, eta') == (differentiate(f, derivation) / f) - (((degree numer f)::Z - (degree denom f)::Z) * eta')::UP::RF notelementary rec == rec.ir + integral(rec.polypart::RF + rec.specpart, monomial(1,1)$UP :: RF) -- returns -- (g in IR, a in F) st f = g'+ a, and a=0 or a has no integral in UP primintegrate(f, derivation, extendedint) == rec := monomialIntegrate(f, derivation) not elem?(i1 := rec.ir) => [notelementary rec, 0] (u2 := primintegratepoly(rec.polypart, extendedint, retract derivation monomial(1, 1))) case UPUP => [i1 + u2.elem::RF::IR + integral(u2.notelem::RF, monomial(1,1)$UP :: RF), 0] [i1 + u2.answer::RF::IR, u2.a0] -- returns -- (g in IR, a in F) st f = g' + a, and a = 0 or a has no integral in F expintegrate(f, derivation, FRDE) == rec := monomialIntegrate(f, derivation) not elem?(i1 := rec.ir) => [notelementary rec, 0] -- rec.specpart is either 0 or of the form p(t)/t**n special := rec.polypart::GP + (numer(rec.specpart)::GP exquo denom(rec.specpart)::GP)::GP (u2 := expintegratepoly(special, FRDE)) case GPGP => [i1 + convert(u2.elem)@RF::IR + integral(convert(u2.notelem)@RF, monomial(1,1)$UP :: RF), 0] [i1 + convert(u2.answer)@RF::IR, u2.a0] -- returns -- (g in IR, a in F) st f = g' + a, and a = 0 or a has no integral in F tanintegrate(f, derivation, FRDE) == rec := monomialIntegrate(f, derivation) not elem?(i1 := rec.ir) => [notelementary rec, 0] r := monomialIntPoly(rec.polypart, derivation) t := monomial(1, 1)$UP c := coefficient(r.polypart, 1) / leadingCoefficient(derivation t) derivation(c::UP) ~= 0 => [i1 + mkAnswer(r.answer::RF, empty(), [[r.polypart::RF + rec.specpart, dummy]$NE]), 0] logs:List(LOG) := zero? c => empty() [[1, monomial(1,1)$UPR - (c/(2::F))::UP::RF::UPR, (1 + t**2)::RF::UPR]] c0 := coefficient(r.polypart, 0) (u := tanintegratespecial(rec.specpart, differentiate(#1, derivation), FRDE)) case RFRF => [i1 + mkAnswer(r.answer::RF + u.elem, logs, [[u.notelem,dummy]$NE]), c0] [i1 + mkAnswer(r.answer::RF + u.answer, logs, empty()), u.a0 + c0] -- returns either (v in RF, c in RF) s.t. f = v' + cg, and c' = 0 -- or (v in RF, a in F) s.t. f = v' + a -- and a = 0 or a has no integral in UP primextendedint(f, derivation, extendedint, g) == fqr := divide(numer f, denom f) gqr := divide(numer g, denom g) (u1 := primextintfrac(fqr.remainder / (denom f), derivation, gqr.remainder / (denom g))) case "failed" => "failed" zero?(gqr.remainder) => -- the following FAIL cannot occur if the primitives are all logs degree(gqr.quotient) > 0 => FAIL (u3 := primintegratepoly(fqr.quotient, extendedint, retract derivation monomial(1, 1))) case UPUP => "failed" [u1.ratpart + u3.answer::RF, u3.a0] (u2 := primintfldpoly(fqr.quotient - retract(u1.coeff)@UP * gqr.quotient, extendedint, retract derivation monomial(1, 1))) case "failed" => "failed" [u2::UP::RF + u1.ratpart, u1.coeff] -- returns either (v in RF, c in RF) s.t. f = v' + cg, and c' = 0 -- or (v in RF, a in F) s.t. f = v' + a -- and a = 0 or a has no integral in F expextendedint(f, derivation, FRDE, g) == qf := separate(f)$GP qg := separate g (u1 := expextintfrac(qf.fracPart, derivation, qg.fracPart)) case "failed" => "failed" zero?(qg.fracPart) => --the following FAIL's cannot occur if the primitives are all logs retractIfCan(qg.polyPart)@Union(F,"failed") case "failed"=> FAIL (u3 := expintegratepoly(qf.polyPart,FRDE)) case GPGP => "failed" [u1.ratpart + convert(u3.answer)@RF, u3.a0] (u2 := expintfldpoly(qf.polyPart - retract(u1.coeff)@UP :: GP * qg.polyPart, FRDE)) case "failed" => "failed" [convert(u2::GP)@RF + u1.ratpart, u1.coeff] -- returns either -- (q in UP, a in F) st p = q'+ a, and a=0 or a has no integral in UP -- or (q in UP, r in UP) st p = q'+ r, and r has no integral elem/UP primintegratepoly(p, extendedint, t') == zero? p => [0, 0$F] ans:UP := 0 while (d := degree p) > 0 repeat (ans1 := extendedint leadingCoefficient p) case "failed" => return([ans, p]) p := reductum p - monomial(d * t' * ans1.ratpart, (d - 1)::N) ans := ans + monomial(ans1.ratpart, d) + monomial(ans1.coeff / (d + 1)::F, d + 1) (ans1:= extendedint(rp := retract(p)@F)) case "failed" => [ans,rp] [monomial(ans1.coeff, 1) + ans1.ratpart::UP + ans, 0$F] -- returns (v in RF, c in RF) s.t. f = v' + cg, and c' = 0 -- g must have a squarefree denominator (always possible) -- g must have no polynomial part (degree numer g < degree denom g) -- f must have no polynomial part (degree numer f < degree denom f) primextintfrac(f, derivation, g) == zero? f => [0, 0] degree numer f >= degree denom f => error "Not a proper fraction" r := HermiteIntegrate(f, derivation) zero? g => r.logpart ~= 0 => "failed" [r.answer, 0] degree numer g >= degree denom g => error "Not a proper fraction" differentiate(c := r.logpart / g, derivation) ~= 0 => "failed" [r.answer, c] @ \section{package INTRAT RationalIntegration} <>= )abbrev package INTRAT RationalIntegration ++ Rational function integration ++ Author: Manuel Bronstein ++ Date Created: 1987 ++ Date Last Updated: 24 October 1995 ++ Description: ++ This package provides functions for the base ++ case of the Risch algorithm. -- Used internally bt the integration packages RationalIntegration(F, UP): Exports == Implementation where F : Join(Field, CharacteristicZero, RetractableTo Integer) UP: UnivariatePolynomialCategory F RF ==> Fraction UP IR ==> IntegrationResult RF LLG ==> List Record(coeff:RF, logand:RF) URF ==> Union(Record(ratpart:RF, coeff:RF), "failed") U ==> Union(Record(mainpart:RF, limitedlogs:LLG), "failed") Exports ==> with integrate : RF -> IR ++ integrate(f) returns g such that \spad{g' = f}. infieldint : RF -> Union(RF, "failed") ++ infieldint(f) returns g such that \spad{g' = f} or "failed" ++ if the integral of f is not a rational function. extendedint: (RF, RF) -> URF ++ extendedint(f, g) returns fractions \spad{[h, c]} such that ++ \spad{c' = 0} and \spad{h' = f - cg}, ++ if \spad{(h, c)} exist, "failed" otherwise. limitedint : (RF, List RF) -> U ++ \spad{limitedint(f, [g1,...,gn])} returns ++ fractions \spad{[h,[[ci, gi]]]} ++ such that the gi's are among \spad{[g1,...,gn]}, \spad{ci' = 0}, and ++ \spad{(h+sum(ci log(gi)))' = f}, if possible, "failed" otherwise. Implementation ==> add import TranscendentalIntegration(F, UP) infieldint f == rec := baseRDE(0, f)$TranscendentalRischDE(F, UP) rec.nosol => "failed" rec.ans integrate f == rec := monomialIntegrate(f, differentiate) integrate(rec.polypart)::RF::IR + rec.ir limitedint(f, lu) == quorem := divide(numer f, denom f) (u := primlimintfrac(quorem.remainder / (denom f), differentiate, lu)) case "failed" => "failed" [u.mainpart + integrate(quorem.quotient)::RF, u.limitedlogs] extendedint(f, g) == fqr := divide(numer f, denom f) gqr := divide(numer g, denom g) (i1 := primextintfrac(fqr.remainder / (denom f), differentiate, gqr.remainder / (denom g))) case "failed" => "failed" i2:=integrate(fqr.quotient-retract(i1.coeff)@UP *gqr.quotient)::RF [i2 + i1.ratpart, i1.coeff] @ \section{package INTRF RationalFunctionIntegration} <>= )abbrev package INTRF RationalFunctionIntegration ++ Integration of rational functions ++ Author: Manuel Bronstein ++ Date Created: 1987 ++ Date Last Updated: 29 Mar 1990 ++ Keywords: polynomial, fraction, integration. ++ Description: ++ This package provides functions for the integration ++ of rational functions. ++ Examples: )r INTRF INPUT RationalFunctionIntegration(F): Exports == Implementation where F: Join(IntegralDomain, RetractableTo Integer, CharacteristicZero) SE ==> Symbol P ==> Polynomial F Q ==> Fraction P UP ==> SparseUnivariatePolynomial Q QF ==> Fraction UP LGQ ==> List Record(coeff:Q, logand:Q) UQ ==> Union(Record(ratpart:Q, coeff:Q), "failed") ULQ ==> Union(Record(mainpart:Q, limitedlogs:LGQ), "failed") Exports ==> with internalIntegrate: (Q, SE) -> IntegrationResult Q ++ internalIntegrate(f, x) returns g such that \spad{dg/dx = f}. infieldIntegrate : (Q, SE) -> Union(Q, "failed") ++ infieldIntegrate(f, x) returns a fraction ++ g such that \spad{dg/dx = f} ++ if g exists, "failed" otherwise. limitedIntegrate : (Q, SE, List Q) -> ULQ ++ \spad{limitedIntegrate(f, x, [g1,...,gn])} returns fractions ++ \spad{[h, [[ci,gi]]]} such that the gi's are among ++ \spad{[g1,...,gn]}, ++ \spad{dci/dx = 0}, and \spad{d(h + sum(ci log(gi)))/dx = f} ++ if possible, "failed" otherwise. extendedIntegrate: (Q, SE, Q) -> UQ ++ extendedIntegrate(f, x, g) returns fractions \spad{[h, c]} such that ++ \spad{dc/dx = 0} and \spad{dh/dx = f - cg}, if \spad{(h, c)} exist, ++ "failed" otherwise. Implementation ==> add import RationalIntegration(Q, UP) import IntegrationResultFunctions2(QF, Q) import PolynomialCategoryQuotientFunctions(IndexedExponents SE, SE, F, P, Q) infieldIntegrate(f, x) == map(multivariate(#1, x), infieldint univariate(f, x)) internalIntegrate(f, x) == map(multivariate(#1, x), integrate univariate(f, x)) extendedIntegrate(f, x, g) == map(multivariate(#1, x), extendedint(univariate(f, x), univariate(g, x))) limitedIntegrate(f, x, lu) == map(multivariate(#1, x), limitedint(univariate(f, x), [univariate(u, x) for u in lu])) @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> -- SPAD files for the integration world should be compiled in the -- following order: -- -- intaux rderf INTRF curve curvepkg divisor pfo -- intalg intaf efstruc rdeef intpm intef irexpand integrat <> <> <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}