1)abbrev package IR2F IntegrationResultToFunction 2++ Conversion of integration results to top-level expressions 3++ Author: Manuel Bronstein 4++ Date Created: 4 February 1988 5++ Description: 6++ This package allows a sum of logs over the roots of a polynomial 7++ to be expressed as explicit logarithms and arc tangents, provided 8++ that the indexing polynomial can be factored into quadratics. 9++ Keywords: integration, expansion, function. 10IntegrationResultToFunction(R, F) : Exports == Implementation where 11 R : Join(GcdDomain, RetractableTo Integer, Comparable, 12 LinearlyExplicitOver Integer) 13 F : Join(AlgebraicallyClosedFunctionSpace R, 14 TranscendentalFunctionCategory) 15 16 N ==> NonNegativeInteger 17 Z ==> Integer 18 Q ==> Fraction Z 19 K ==> Kernel F 20 P ==> SparseMultivariatePolynomial(R, K) 21 UP ==> SparseUnivariatePolynomial F 22 IR ==> IntegrationResult F 23 REC ==> Record(ans1 : F, ans2 : F) 24 LOG ==> Record(scalar : Q, coeff : UP, logand : UP) 25 26 Exports ==> with 27 split : IR -> IR 28 ++ split(u(x) + sum_{P(a)=0} Q(a, x)) returns 29 ++ \spad{u(x) + sum_{P1(a)=0} Q(a, x) + ... + sum_{Pn(a)=0} Q(a, x)} 30 ++ where P1, ..., Pn are the factors of P. 31 expand : (IR, Symbol) -> List F 32 ++ expand(i, x) returns the list of possible real functions 33 ++ of x corresponding to i. 34 complexExpand : IR -> F 35 ++ complexExpand(i) returns the expanded complex function 36 ++ corresponding to i. 37 38 Implementation ==> add 39 import from AlgebraicManipulations(R, F) 40 import from ElementaryFunctionSign(R, F) 41 42 IR2F : IR -> F 43 insqrt : F -> Record(sqrt : REC, sgn : Z) 44 pairsum : (List F, List F) -> List F 45 pairprod : (F, List F) -> List F 46 quadeval : (UP, F, F, F) -> REC 47 linear : (UP, UP) -> F 48 tantrick : (F, F) -> F 49 ilog : (F, F, Symbol) -> F 50 ilog0 : (F, F, UP, UP, F) -> F 51 nlogs : LOG -> List LOG 52 lg2func : (LOG, Symbol) -> List F 53 quadratic : (UP, UP, Symbol) -> List F 54 mkRealFunc : (List LOG, Symbol) -> List F 55 lg2cfunc : LOG -> F 56 loglist : (Q, UP, UP) -> List LOG 57 cmplex : (F, UP) -> F 58 evenRoots : F -> List F 59 compatible? : (List F, List F) -> Boolean 60 61 cmplex(alpha, p) == alpha * log p alpha 62 IR2F i == retract mkAnswer(ratpart i, empty(), notelem i) 63 pairprod(x, l) == [x * y for y in l] 64 65 evenRoots x == 66 [first argument k for k in tower x | 67 is?(k, 'nthRoot) and even?(retract(second argument k)@Z) 68 and (not empty? variables first argument k)] 69 70 expand(i, x) == 71 j := split i 72 pairsum([IR2F j], mkRealFunc(logpart(j), x)) 73 74 split i == 75 mkAnswer(ratpart i, concat [nlogs l for l in logpart i], notelem i) 76 77 complexExpand i == 78 j := split i 79 IR2F j + +/[lg.scalar::F * lg2cfunc lg for lg in logpart j] 80 81-- p = a t^2 + b t + p0 82-- Expands sum_{p(t) = 0} t log(lg(t)) 83 quadratic(p, lg, x) == 84 zero?(delta := (b := coefficient(p, 1))^2 - 4 * 85 (a := coefficient(p, 2)) * (p0 := coefficient(p, 0))) => 86 [linear(monomial(1, 1) + (b / a)::UP, lg)] 87 e := (q := quadeval(lg, c := - b * (d := inv(2*a)), d, delta)).ans1 88 lgp := c * log(nrm := (e^2 - delta * (f := q.ans2)^2)) 89 s := (sqr := insqrt delta).sqrt 90 pp := nn := 0$F 91 if sqr.sgn >= 0 then 92 sqrp := s.ans1 * rootSimp sqrt(s.ans2) 93 pp := lgp + d * sqrp * log(((2 * e * f) / nrm) * sqrp 94 + (e^2 + delta * f^2) / nrm) 95 if sqr.sgn <= 0 then 96 sqrn := s.ans1 * rootSimp sqrt(-s.ans2) 97 nn := lgp + d * sqrn * ilog(e, f * sqrn, x) 98 sqr.sgn > 0 => [pp] 99 sqr.sgn < 0 => [nn] 100 [pp, nn] 101 102-- returns 2 atan(a/b) or 2 atan(-b/a) whichever looks better 103-- they differ by a constant so it's ok to do it from an IR 104 tantrick(a, b) == 105 retractIfCan(a)@Union(Q, "failed") case Q => 2 * atan(-b/a) 106 sb := sign b 107 if sb case Z then return 2 * atan(a/b) 108 sa := sign a 109 if sa case Z then return 2 * atan(-b/a) 110 -- print("potentially noncontinuous"::OutputForm) 111 2 * atan(a/b) 112 113 var_kers(p : P, x : Symbol) : List(K) == 114 [k for k in variables(p) | D(k::F, x) ~= 0] 115 116 -- transforms i log((a + i b) / (a - i b)) into a sum of real 117 -- arc-tangents using Rioboo's algorithm 118 -- We only apply transform if we hope for improvement 119 ilog(a, b, x) == 120 b = 0 => 0 121 not(empty?(var_kers(denom a, x))) or 122 not(empty?(var_kers(denom b, x))) => tantrick(a, b) 123 l := setUnion(var_kers(numer a, x), var_kers(numer b, x)) 124 #l ~= 1 => tantrick(a, b) 125 k := first(l) 126 opk := operator(k) 127 is?(opk, 'tan) or has?(opk, '%alg) => tantrick(a, b) 128 ilog0(a, b, numer univariate(a, k), numer univariate(b, k), k::F) 129 130-- transforms i log((a + i b) / (a - i b)) into a sum of real 131-- arc-tangents using Rioboo's algorithm 132-- the arc-tangents will not have k in the denominator 133-- we always keep upa(k) = a and upb(k) = b 134 ilog0(a, b, upa, upb, k) == 135 if degree(upa) < degree(upb) then 136 (upa, upb) := (-upb, upa) 137 (a, b) := (-b, a) 138 zero? degree upb => tantrick(a, b) 139 r := extendedEuclidean(upa, upb) 140 (g:= retractIfCan(r.generator)@Union(F,"failed")) case "failed" => 141 tantrick(a, b) 142 if degree(r.coef1) >= degree upb then 143 qr := divide(r.coef1, upb) 144 r.coef1 := qr.remainder 145 r.coef2 := r.coef2 + qr.quotient * upa 146 aa := (r.coef2) k 147 bb := -(r.coef1) k 148 tantrick(aa * a + bb * b, g::F) + ilog0(aa, bb, r.coef2, -r.coef1, k) 149 150 lg2func(lg, x) == 151 zero?(d := degree(p := lg.coeff)) => error "poly has degree 0" 152 (d = 1) => [linear(p, lg.logand)] 153 d = 2 => quadratic(p, lg.logand, x) 154 odd? d and 155 ((r := retractIfCan(reductum p)@Union(F, "failed")) case F) => 156 pairsum([cmplex(alpha := rootSimp zeroOf p, lg.logand)], 157 lg2func([lg.scalar, 158 (p exquo (monomial(1, 1)$UP - alpha::UP))::UP, 159 lg.logand], x)) 160 [lg2cfunc lg] 161 162 lg2cfunc lg == 163 +/[cmplex(alpha, lg.logand) for alpha in zerosOf(lg.coeff)] 164 165 mkRealFunc(l, x) == 166 ans := empty()$List(F) 167 for lg in l repeat 168 ans := pairsum(ans, pairprod(lg.scalar::F, lg2func(lg, x))) 169 ans 170 171-- returns a log(b) 172 linear(p, lg) == 173 alpha := - coefficient(p, 0) / coefficient(p, 1) 174 alpha * log lg alpha 175 176-- returns (c, d) s.t. p(a + b t) = c + d t, where t^2 = delta 177 quadeval(p, a, b, delta) == 178 zero? p => [0, 0] 179 bi := c := d := 0$F 180 ai := 1$F 181 v := vectorise(p, 1 + degree p) 182 for i in minIndex v .. maxIndex v repeat 183 c := c + qelt(v, i) * ai 184 d := d + qelt(v, i) * bi 185 temp := a * ai + b * bi * delta 186 bi := a * bi + b * ai 187 ai := temp 188 [c, d] 189 190 compatible?(lx, ly) == 191 empty? ly => true 192 for x in lx repeat 193 for y in ly repeat 194 ((s := sign(x*y)) case Z) and (s::Z < 0) => return false 195 true 196 197 pairsum(lx, ly) == 198 empty? lx => ly 199 empty? ly => lx 200 l := empty()$List(F) 201 for x in lx repeat 202 ls := evenRoots x 203 if not empty?(ln := 204 [x + y for y in ly | compatible?(ls, evenRoots y)])$List(F) then 205 l := removeDuplicates concat(l, ln) 206 l 207 208-- returns [[a, b], s] where sqrt(y) = a sqrt(b) and 209-- s = 1 if b > 0, -1 if b < 0, 0 if the sign of b cannot be determined 210 insqrt y == 211 rec := froot(y, 2)$PolynomialRoots(IndexedExponents K, K, R, P, F) 212 ((rec.exponent) = 1) => [[rec.coef * rec.radicand, 1], 1] 213 rec.exponent ~= 2 => error "Should not happen" 214 [[rec.coef, rec.radicand], 215 ((s := sign(rec.radicand)) case "failed" => 0; s::Z)] 216 217 nlogs lg == 218 [[f.exponent * lg.scalar, f.factor, lg.logand] for f in factorList 219 factorPolynomial(primitivePart(lg.coeff) 220 )$ExpressionFactorPolynomial(R, F)] 221 222)abbrev package IRRF2F IntegrationResultRFToFunction 223++ Conversion of integration results to top-level expressions 224++ Author: Manuel Bronstein 225++ Description: 226++ This package allows a sum of logs over the roots of a polynomial 227++ to be expressed as explicit logarithms and arc tangents, provided 228++ that the indexing polynomial can be factored into quadratics. 229++ Date Created: 21 August 1988 230IntegrationResultRFToFunction(R) : Exports == Implementation where 231 R : Join(GcdDomain, RetractableTo Integer, Comparable, 232 LinearlyExplicitOver Integer) 233 234 RF ==> Fraction Polynomial R 235 F ==> Expression R 236 IR ==> IntegrationResult RF 237 238 Exports ==> with 239 split : IR -> IR 240 ++ split(u(x) + sum_{P(a)=0} Q(a, x)) returns 241 ++ \spad{u(x) + sum_{P1(a)=0} Q(a, x) + ... + sum_{Pn(a)=0} Q(a, x)} 242 ++ where P1, ..., Pn are the factors of P. 243 expand : (IR, Symbol) -> List F 244 ++ expand(i, x) returns the list of possible real functions 245 ++ of x corresponding to i. 246 complexExpand : IR -> F 247 ++ complexExpand(i) returns the expanded complex function 248 ++ corresponding to i. 249 if R has CharacteristicZero then 250 integrate : (RF, Symbol) -> Union(F, List F) 251 ++ integrate(f, x) returns the integral of \spad{f(x)dx} 252 ++ where x is viewed as a real variable. 253 complexIntegrate : (RF, Symbol) -> F 254 ++ complexIntegrate(f, x) returns the integral of \spad{f(x)dx} 255 ++ where x is viewed as a complex variable. 256 257 Implementation ==> add 258 import from IntegrationTools(R, F) 259 import from TrigonometricManipulations(R, F) 260 import from IntegrationResultToFunction(R, F) 261 262 toEF : IR -> IntegrationResult F 263 264 toEF i == map(z1 +-> z1::F, i)$IntegrationResultFunctions2(RF, F) 265 expand(i, x) == expand(toEF(i), x) 266 complexExpand i == complexExpand toEF i 267 268 split i == 269 map(retract, split toEF i)$IntegrationResultFunctions2(F, RF) 270 271 if R has CharacteristicZero then 272 import from RationalFunctionIntegration(R) 273 274 complexIntegrate(f, x) == complexExpand internalIntegrate(f, x) 275 276-- do not use real integration if R is complex 277 if R has imaginary : () -> R then integrate(f, x) == complexIntegrate(f, x) 278 else 279 integrate(f, x) == 280 l := [mkPrim(real g, x) for g in expand(internalIntegrate(f, x), x)] 281 empty? rest l => first l 282 l 283 284--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. 285--All rights reserved. 286-- 287--Redistribution and use in source and binary forms, with or without 288--modification, are permitted provided that the following conditions are 289--met: 290-- 291-- - Redistributions of source code must retain the above copyright 292-- notice, this list of conditions and the following disclaimer. 293-- 294-- - Redistributions in binary form must reproduce the above copyright 295-- notice, this list of conditions and the following disclaimer in 296-- the documentation and/or other materials provided with the 297-- distribution. 298-- 299-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the 300-- names of its contributors may be used to endorse or promote products 301-- derived from this software without specific prior written permission. 302-- 303--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 304--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 305--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 306--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 307--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 308--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 309--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 310--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 311--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 312--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 313--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 314