\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra irexpand.spad} \author{Manuel Bronstein} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package IR2F IntegrationResultToFunction} <>= )abbrev package IR2F IntegrationResultToFunction ++ Conversion of integration results to top-level expressions ++ Author: Manuel Bronstein ++ Date Created: 4 February 1988 ++ Date Last Updated: 9 October 1991 ++ 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, LinearlyExplicitRingOver 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 -> List F ++ expand(i) returns the list of possible real functions ++ corresponding to i. complexExpand: IR -> F ++ complexExpand(i) returns the expanded complex function ++ corresponding to i. Implementation ==> add import AlgebraicManipulations(R, F) import 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, List K) -> F ilog0 : (F, F, UP, UP, F) -> F nlogs : LOG -> List LOG lg2func : LOG -> List F quadratic : (UP, UP) -> List F mkRealFunc : List LOG -> 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 == j := split i pairsum([IR2F j], mkRealFunc logpart j) 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 + c -- Expands sum_{p(t) = 0} t log(lg(t)) quadratic(p, lg) == 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, setUnion(setUnion(kernels a, kernels b), kernels p0)) positive? sqr.sgn => [pp] negative? sqr.sgn => [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) 2 * atan(a/b) -- transforms i log((a + i b) / (a - i b)) into a sum of real -- arc-tangents using Rioboo's algorithm -- lk is a list of kernels which are parameters for the integral ilog(a, b, lk) == l := setDifference(setUnion(variables numer a, variables numer b), setUnion(lk, setUnion(variables denom a, variables denom b))) empty? l => tantrick(a, b) k := "max"/l 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 == zero?(d := degree(p := lg.coeff)) => error "poly has degree 0" one? d => [linear(p, lg.logand)] d = 2 => quadratic(p, lg.logand) 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]) [lg2cfunc lg] lg2cfunc lg == +/[cmplex(alpha, lg.logand) for alpha in zerosOf(lg.coeff)] mkRealFunc l == ans := empty()$List(F) for lg in l repeat ans := pairsum(ans, pairprod(lg.scalar::F, lg2func lg)) 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 negative?(s::Z) => 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)]) 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) one?(rec.exponent) => [[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 factors ffactor(primitivePart(lg.coeff) )$FunctionSpaceUnivariatePolynomialFactor(R, F, UP)] @ \section{package IRRF2F IntegrationResultRFToFunction} <>= )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 ++ Date Last Updated: 4 October 1993 IntegrationResultRFToFunction(R): Exports == Implementation where R: Join(GcdDomain, RetractableTo Integer, LinearlyExplicitRingOver 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 -> List F ++ expand(i) returns the list of possible real functions ++ 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 IntegrationTools(R, F) import TrigonometricManipulations(R, F) import IntegrationResultToFunction(R, F) toEF: IR -> IntegrationResult F toEF i == map(#1::F, i)$IntegrationResultFunctions2(RF, F) expand i == expand toEF i complexExpand i == complexExpand toEF i split i == map(retract, split toEF i)$IntegrationResultFunctions2(F, RF) if R has CharacteristicZero then import 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)] empty? rest l => first l l @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. --Copyright (C) 2007-2009, Gabriel Dos Reis. --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 intef IREXPAND integrat <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}