)abbrev package IR2F IntegrationResultToFunction ++ Conversion of integration results to top-level expressions ++ Author: Manuel Bronstein ++ Date Created: 4 February 1988 ++ Description: ++ This package allows a sum of logs over the roots of a polynomial ++ to be expressed as explicit logarithms and arc tangents, provided ++ that the indexing polynomial can be factored into quadratics. ++ Keywords: integration, expansion, function. IntegrationResultToFunction(R, F) : Exports == Implementation where R : Join(GcdDomain, RetractableTo Integer, Comparable, LinearlyExplicitOver Integer) F : Join(AlgebraicallyClosedFunctionSpace R, TranscendentalFunctionCategory) N ==> NonNegativeInteger Z ==> Integer Q ==> Fraction Z K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial F IR ==> IntegrationResult F REC ==> Record(ans1 : F, ans2 : F) LOG ==> Record(scalar : Q, coeff : UP, logand : UP) Exports ==> with split : IR -> IR ++ split(u(x) + sum_{P(a)=0} Q(a, x)) returns ++ \spad{u(x) + sum_{P1(a)=0} Q(a, x) + ... + sum_{Pn(a)=0} Q(a, x)} ++ where P1, ..., Pn are the factors of P. expand : (IR, Symbol) -> List F ++ expand(i, x) returns the list of possible real functions ++ of x corresponding to i. complexExpand : IR -> F ++ complexExpand(i) returns the expanded complex function ++ corresponding to i. Implementation ==> add import from AlgebraicManipulations(R, F) import from ElementaryFunctionSign(R, F) IR2F : IR -> F insqrt : F -> Record(sqrt : REC, sgn : Z) pairsum : (List F, List F) -> List F pairprod : (F, List F) -> List F quadeval : (UP, F, F, F) -> REC linear : (UP, UP) -> F tantrick : (F, F) -> F ilog : (F, F, Symbol) -> F ilog0 : (F, F, UP, UP, F) -> F nlogs : LOG -> List LOG lg2func : (LOG, Symbol) -> List F quadratic : (UP, UP, Symbol) -> List F mkRealFunc : (List LOG, Symbol) -> List F lg2cfunc : LOG -> F loglist : (Q, UP, UP) -> List LOG cmplex : (F, UP) -> F evenRoots : F -> List F compatible? : (List F, List F) -> Boolean cmplex(alpha, p) == alpha * log p alpha IR2F i == retract mkAnswer(ratpart i, empty(), notelem i) pairprod(x, l) == [x * y for y in l] evenRoots x == [first argument k for k in tower x | is?(k, 'nthRoot) and even?(retract(second argument k)@Z) and (not empty? variables first argument k)] expand(i, x) == j := split i pairsum([IR2F j], mkRealFunc(logpart(j), x)) split i == mkAnswer(ratpart i, concat [nlogs l for l in logpart i], notelem i) complexExpand i == j := split i IR2F j + +/[lg.scalar::F * lg2cfunc lg for lg in logpart j] -- p = a t^2 + b t + p0 -- Expands sum_{p(t) = 0} t log(lg(t)) quadratic(p, lg, x) == zero?(delta := (b := coefficient(p, 1))^2 - 4 * (a := coefficient(p, 2)) * (p0 := coefficient(p, 0))) => [linear(monomial(1, 1) + (b / a)::UP, lg)] e := (q := quadeval(lg, c := - b * (d := inv(2*a)), d, delta)).ans1 lgp := c * log(nrm := (e^2 - delta * (f := q.ans2)^2)) s := (sqr := insqrt delta).sqrt pp := nn := 0$F if sqr.sgn >= 0 then sqrp := s.ans1 * rootSimp sqrt(s.ans2) pp := lgp + d * sqrp * log(((2 * e * f) / nrm) * sqrp + (e^2 + delta * f^2) / nrm) if sqr.sgn <= 0 then sqrn := s.ans1 * rootSimp sqrt(-s.ans2) nn := lgp + d * sqrn * ilog(e, f * sqrn, x) sqr.sgn > 0 => [pp] sqr.sgn < 0 => [nn] [pp, nn] -- returns 2 atan(a/b) or 2 atan(-b/a) whichever looks better -- they differ by a constant so it's ok to do it from an IR tantrick(a, b) == retractIfCan(a)@Union(Q, "failed") case Q => 2 * atan(-b/a) sb := sign b if sb case Z then return 2 * atan(a/b) sa := sign a if sa case Z then return 2 * atan(-b/a) -- print("potentially noncontinuous"::OutputForm) 2 * atan(a/b) var_kers(p : P, x : Symbol) : List(K) == [k for k in variables(p) | D(k::F, x) ~= 0] -- transforms i log((a + i b) / (a - i b)) into a sum of real -- arc-tangents using Rioboo's algorithm -- We only apply transform if we hope for improvement ilog(a, b, x) == b = 0 => 0 not(empty?(var_kers(denom a, x))) or not(empty?(var_kers(denom b, x))) => tantrick(a, b) l := setUnion(var_kers(numer a, x), var_kers(numer b, x)) #l ~= 1 => tantrick(a, b) k := first(l) opk := operator(k) is?(opk, 'tan) or has?(opk, '%alg) => tantrick(a, b) ilog0(a, b, numer univariate(a, k), numer univariate(b, k), k::F) -- transforms i log((a + i b) / (a - i b)) into a sum of real -- arc-tangents using Rioboo's algorithm -- the arc-tangents will not have k in the denominator -- we always keep upa(k) = a and upb(k) = b ilog0(a, b, upa, upb, k) == if degree(upa) < degree(upb) then (upa, upb) := (-upb, upa) (a, b) := (-b, a) zero? degree upb => tantrick(a, b) r := extendedEuclidean(upa, upb) (g:= retractIfCan(r.generator)@Union(F,"failed")) case "failed" => tantrick(a, b) if degree(r.coef1) >= degree upb then qr := divide(r.coef1, upb) r.coef1 := qr.remainder r.coef2 := r.coef2 + qr.quotient * upa aa := (r.coef2) k bb := -(r.coef1) k tantrick(aa * a + bb * b, g::F) + ilog0(aa, bb, r.coef2, -r.coef1, k) lg2func(lg, x) == zero?(d := degree(p := lg.coeff)) => error "poly has degree 0" (d = 1) => [linear(p, lg.logand)] d = 2 => quadratic(p, lg.logand, x) odd? d and ((r := retractIfCan(reductum p)@Union(F, "failed")) case F) => pairsum([cmplex(alpha := rootSimp zeroOf p, lg.logand)], lg2func([lg.scalar, (p exquo (monomial(1, 1)$UP - alpha::UP))::UP, lg.logand], x)) [lg2cfunc lg] lg2cfunc lg == +/[cmplex(alpha, lg.logand) for alpha in zerosOf(lg.coeff)] mkRealFunc(l, x) == ans := empty()$List(F) for lg in l repeat ans := pairsum(ans, pairprod(lg.scalar::F, lg2func(lg, x))) ans -- returns a log(b) linear(p, lg) == alpha := - coefficient(p, 0) / coefficient(p, 1) alpha * log lg alpha -- returns (c, d) s.t. p(a + b t) = c + d t, where t^2 = delta quadeval(p, a, b, delta) == zero? p => [0, 0] bi := c := d := 0$F ai := 1$F v := vectorise(p, 1 + degree p) for i in minIndex v .. maxIndex v repeat c := c + qelt(v, i) * ai d := d + qelt(v, i) * bi temp := a * ai + b * bi * delta bi := a * bi + b * ai ai := temp [c, d] compatible?(lx, ly) == empty? ly => true for x in lx repeat for y in ly repeat ((s := sign(x*y)) case Z) and (s::Z < 0) => return false true pairsum(lx, ly) == empty? lx => ly empty? ly => lx l := empty()$List(F) for x in lx repeat ls := evenRoots x if not empty?(ln := [x + y for y in ly | compatible?(ls, evenRoots y)])$List(F) then l := removeDuplicates concat(l, ln) l -- returns [[a, b], s] where sqrt(y) = a sqrt(b) and -- s = 1 if b > 0, -1 if b < 0, 0 if the sign of b cannot be determined insqrt y == rec := froot(y, 2)$PolynomialRoots(IndexedExponents K, K, R, P, F) ((rec.exponent) = 1) => [[rec.coef * rec.radicand, 1], 1] rec.exponent ~= 2 => error "Should not happen" [[rec.coef, rec.radicand], ((s := sign(rec.radicand)) case "failed" => 0; s::Z)] nlogs lg == [[f.exponent * lg.scalar, f.factor, lg.logand] for f in factorList factorPolynomial(primitivePart(lg.coeff) )$ExpressionFactorPolynomial(R, F)] )abbrev package IRRF2F IntegrationResultRFToFunction ++ Conversion of integration results to top-level expressions ++ Author: Manuel Bronstein ++ Description: ++ This package allows a sum of logs over the roots of a polynomial ++ to be expressed as explicit logarithms and arc tangents, provided ++ that the indexing polynomial can be factored into quadratics. ++ Date Created: 21 August 1988 IntegrationResultRFToFunction(R) : Exports == Implementation where R : Join(GcdDomain, RetractableTo Integer, Comparable, LinearlyExplicitOver Integer) RF ==> Fraction Polynomial R F ==> Expression R IR ==> IntegrationResult RF Exports ==> with split : IR -> IR ++ split(u(x) + sum_{P(a)=0} Q(a, x)) returns ++ \spad{u(x) + sum_{P1(a)=0} Q(a, x) + ... + sum_{Pn(a)=0} Q(a, x)} ++ where P1, ..., Pn are the factors of P. expand : (IR, Symbol) -> List F ++ expand(i, x) returns the list of possible real functions ++ of x corresponding to i. complexExpand : IR -> F ++ complexExpand(i) returns the expanded complex function ++ corresponding to i. if R has CharacteristicZero then integrate : (RF, Symbol) -> Union(F, List F) ++ integrate(f, x) returns the integral of \spad{f(x)dx} ++ where x is viewed as a real variable. complexIntegrate : (RF, Symbol) -> F ++ complexIntegrate(f, x) returns the integral of \spad{f(x)dx} ++ where x is viewed as a complex variable. Implementation ==> add import from IntegrationTools(R, F) import from TrigonometricManipulations(R, F) import from IntegrationResultToFunction(R, F) toEF : IR -> IntegrationResult F toEF i == map(z1 +-> z1::F, i)$IntegrationResultFunctions2(RF, F) expand(i, x) == expand(toEF(i), x) complexExpand i == complexExpand toEF i split i == map(retract, split toEF i)$IntegrationResultFunctions2(F, RF) if R has CharacteristicZero then import from RationalFunctionIntegration(R) complexIntegrate(f, x) == complexExpand internalIntegrate(f, x) -- do not use real integration if R is complex if R has imaginary : () -> R then integrate(f, x) == complexIntegrate(f, x) else integrate(f, x) == l := [mkPrim(real g, x) for g in expand(internalIntegrate(f, x), x)] empty? rest l => first l l --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.