\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra intpm.spad} \author{Manuel Bronstein} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package INTPM PatternMatchIntegration} <>= )abbrev package INTPM PatternMatchIntegration ++ Author: Manuel Bronstein ++ Date Created: 5 May 1992 ++ Date Last Updated: 27 September 1995 ++ Description: ++ \spadtype{PatternMatchIntegration} provides functions that use ++ the pattern matcher to find some indefinite and definite integrals ++ involving special functions and found in the litterature. PatternMatchIntegration(R, F): Exports == Implementation where R : Join(RetractableTo Integer, GcdDomain, LinearlyExplicitRingOver Integer) F : Join(AlgebraicallyClosedField, TranscendentalFunctionCategory, FunctionSpace R) N ==> NonNegativeInteger Z ==> Integer SY ==> Symbol K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) SUP ==> SparseUnivariatePolynomial F PAT ==> Pattern Z RES ==> PatternMatchResult(Z, F) OFE ==> OrderedCompletion F REC ==> Record(which: Z, exponent: F, coeff: F) ANS ==> Record(special:F, integrand:F) NONE ==> 0 EI ==> 1 ERF ==> 2 SI ==> 3 CI ==> 4 GAM2 ==> 5 CI0 ==> 6 Exports ==> with splitConstant: (F, SY) -> Record(const:F, nconst:F) ++ splitConstant(f, x) returns \spad{[c, g]} such that ++ \spad{f = c * g} and \spad{c} does not involve \spad{t}. if R has ConvertibleTo Pattern Integer and R has PatternMatchable Integer then if F has LiouvillianFunctionCategory then pmComplexintegrate: (F, SY) -> Union(ANS, "failed") ++ pmComplexintegrate(f, x) returns either "failed" or ++ \spad{[g,h]} such that ++ \spad{integrate(f,x) = g + integrate(h,x)}. ++ It only looks for special complex integrals that pmintegrate ++ does not return. pmintegrate: (F, SY) -> Union(ANS, "failed") ++ pmintegrate(f, x) returns either "failed" or \spad{[g,h]} such ++ that \spad{integrate(f,x) = g + integrate(h,x)}. if F has SpecialFunctionCategory then pmintegrate: (F, SY, OFE, OFE) -> Union(F, "failed") ++ pmintegrate(f, x = a..b) returns the integral of ++ \spad{f(x)dx} from a to b ++ if it can be found by the built-in pattern matching rules. Implementation ==> add import PatternMatch(Z, F, F) import ElementaryFunctionSign(R, F) import FunctionSpaceAssertions(R, F) import TrigonometricManipulations(R, F) import FunctionSpaceAttachPredicates(R, F, F) mkalist : RES -> AssociationList(SY, F) pm := new()$SY pmw := new pm pmm := new pm pms := new pm pmc := new pm pma := new pm pmb := new pm c := optional(pmc::F) w := suchThat(optional(pmw::F), empty? variables #1) s := suchThat(optional(pms::F), empty? variables #1 and real? #1) m := suchThat(optional(pmm::F), (retractIfCan(#1)@Union(Z,"failed") case Z) and not before?(#1,0)) spi := sqrt(pi()$F) half := 1::F / 2::F mkalist res == construct destruct res splitConstant(f, x) == not member?(x, variables f) => [f, 1] (retractIfCan(f)@Union(K, "failed")) case K => [1, f] (u := isTimes f) case List(F) => cc := nc := 1$F for g in u::List(F) repeat rec := splitConstant(g, x) cc := cc * rec.const nc := nc * rec.nconst [cc, nc] (u := isPlus f) case List(F) => rec := splitConstant(first(u::List(F)), x) cc := rec.const nc := rec.nconst for g in rest(u::List(F)) repeat rec := splitConstant(g, x) if rec.nconst = nc then cc := cc + rec.const else if rec.nconst = -nc then cc := cc - rec.const else return [1, f] [cc, nc] if (v := isPower f) case Record(val:F, exponent:Z) then vv := v::Record(val:F, exponent:Z) not one?(vv.exponent) => rec := splitConstant(vv.val, x) return [rec.const ** vv.exponent, rec.nconst ** vv.exponent] error "splitConstant: should not happen" if R has ConvertibleTo Pattern Integer and R has PatternMatchable Integer then if F has LiouvillianFunctionCategory then insqrt : F -> F matchei : (F, SY) -> REC matcherfei : (F, SY, Boolean) -> REC matchsici : (F, SY) -> REC matchli : (F, SY) -> List F matchli0 : (F, K, SY) -> List F matchdilog : (F, SY) -> List F matchdilog0: (F, K, SY, P, F) -> List F goodlilog? : (K, P) -> Boolean gooddilog? : (K, P, P) -> Boolean goodlilog?(k, p) == is?(k, 'log) and one? minimumDegree(p, k) gooddilog?(k, p, q) == is?(k, 'log) and one? degree(p, k) and zero? degree(q, k) -- matches the integral to a result of the form d * erf(u) or d * ei(u) -- returns [case, u, d] matcherfei(f, x, comp?) == res0 := new()$RES pat := c * exp(pma::F) failed?(res := patternMatch(f, convert(pat)@PAT, res0)) => comp? => [NONE, 0,0] matchei(f,x) l := mkalist res da := differentiate(a := l.pma, x) d := a * (cc := l.pmc) / da zero? differentiate(d, x) => [EI, a, d] comp? or (((u := sign a) case Z) and negative?(u::Z)) => d := cc * (sa := insqrt(- a)) / da zero? differentiate(d, x) => [ERF, sa, - d * spi] [NONE, 0, 0] [NONE, 0, 0] -- matches the integral to a result of the form d * ei(k * log u) -- returns [case, k * log u, d] matchei(f, x) == res0 := new()$RES a := pma::F pat := c * a**w / log a failed?(res := patternMatch(f, convert(pat)@PAT, res0)) => [NONE, 0, 0] l := mkalist res da := differentiate(a := l.pma, x) d := (cc := l.pmc) / da zero? differentiate(d, x) => [EI, (1 + l.pmw) * log a, d] [NONE, 0, 0] -- matches the integral to a result of the form d * dilog(u) + int(v), -- returns [u,d,v] or [] matchdilog(f, x) == n := numer f df := (d := denom f)::F for k in select!(gooddilog?(#1, n, d), variables n)$List(K) repeat not empty?(l := matchdilog0(f, k, x, n, df)) => return l empty() -- matches the integral to a result of the form d * dilog(a) + int(v) -- where k = log(a) -- returns [a,d,v] or [] matchdilog0(f, k, x, p, q) == zero?(da := differentiate(a := first argument k, x)) => empty() a1 := 1 - a d := coefficient(univariate(p, k), 1)::F * a1 / (q * da) zero? differentiate(d, x) => [a, d, f - d * da * (k::F) / a1] empty() -- matches the integral to a result of the form d * li(u) + int(v), -- returns [u,d,v] or [] matchli(f, x) == d := denom f for k in select!(goodlilog?(#1, d), variables d)$List(K) repeat not empty?(l := matchli0(f, k, x)) => return l empty() -- matches the integral to a result of the form d * li(a) + int(v) -- where k = log(a) -- returns [a,d,v] or [] matchli0(f, k, x) == g := (lg := k::F) * f zero?(da := differentiate(a := first argument k, x)) => empty() zero? differentiate(d := g / da, x) => [a, d, 0] ug := univariate(g, k) (u:=retractIfCan(ug)@Union(SUP,"failed")) case "failed" => empty() degree(p := u::SUP) > 1 => empty() zero? differentiate(d := coefficient(p, 0) / da, x) => [a, d, leadingCoefficient p] empty() -- matches the integral to a result of the form d * Si(u) or d * Ci(u) -- returns [case, u, d] matchsici(f, x) == res0 := new()$RES b := pmb::F t := tan(a := pma::F) patsi := c * t / (patden := b + b * t**2) patci := (c - c * t**2) / patden patci0 := c / patden ci0?:Boolean (ci? := failed?(res := patternMatch(f, convert(patsi)@PAT, res0))) and (ci0?:=failed?(res:=patternMatch(f,convert(patci)@PAT,res0))) and failed?(res := patternMatch(f,convert(patci0)@PAT,res0)) => [NONE, 0, 0] l := mkalist res (b := l.pmb) ~= 2 * (a := l.pma) => [NONE, 0, 0] db := differentiate(b, x) d := (cc := l.pmc) / db zero? differentiate(d, x) => ci? => ci0? => [CI0, b, d / (2::F)] [CI, b, d] [SI, b, d / (2::F)] [NONE, 0, 0] -- returns a simplified sqrt(y) insqrt y == rec := froot(y, 2)$PolynomialRoots(IndexedExponents K, K, R, P, F) one?(rec.exponent) => rec.coef * rec.radicand rec.exponent ~=2 => error "insqrt: hould not happen" rec.coef * sqrt(rec.radicand) pmintegrate(f, x) == not one?((rc := splitConstant(f, x)).const) => (u := pmintegrate(rc.nconst, x)) case "failed" => "failed" rec := u::ANS [rc.const * rec.special, rc.const * rec.integrand] not empty?(l := matchli(f, x)) => [second l * li first l, third l] not empty?(l := matchdilog(f, x)) => [second l * dilog first l, third l] cse := (rec := matcherfei(f, x, false)).which cse = EI => [rec.coeff * Ei(rec.exponent), 0] cse = ERF => [rec.coeff * erf(rec.exponent), 0] cse := (rec := matchsici(f, x)).which cse = SI => [rec.coeff * Si(rec.exponent), 0] cse = CI => [rec.coeff * Ci(rec.exponent), 0] cse = CI0 => [rec.coeff * Ci(rec.exponent) + rec.coeff * log(rec.exponent), 0] "failed" pmComplexintegrate(f, x) == not one?((rc := splitConstant(f, x)).const) => (u := pmintegrate(rc.nconst, x)) case "failed" => "failed" rec := u::ANS [rc.const * rec.special, rc.const * rec.integrand] cse := (rec := matcherfei(f, x, true)).which cse = ERF => [rec.coeff * erf(rec.exponent), 0] "failed" if F has SpecialFunctionCategory then match1 : (F, SY, F, F) -> List F formula1 : (F, SY, F, F) -> Union(F, "failed") -- tries only formula (1) of the Geddes & al, AAECC 1 (1990) paper formula1(f, x, t, cc) == empty?(l := match1(f, x, t, cc)) => "failed" mw := first l zero?(ms := third l) or ((sgs := sign ms) case "failed")=> "failed" ((sgz := sign(z := (mw + 1) / ms)) case "failed") or negative?(sgz::Z) => "failed" mmi := retract(mm := second l)@Z sgs * (last l) * ms**(- mmi - 1) * eval(differentiate(Gamma(x::F), x, mmi::N), [kernel(x)@K], [z]) -- returns [w, m, s, c] or [] -- matches only formula (1) of the Geddes & al, AAECC 1 (1990) paper match1(f, x, t, cc) == res0 := new()$RES pat := cc * log(t)**m * exp(-t**s) not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) => l := mkalist res [0, l.pmm, l.pms, l.pmc] pat := cc * t**w * exp(-t**s) not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) => l := mkalist res [l.pmw, 0, l.pms, l.pmc] pat := cc / t**w * exp(-t**s) not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) => l := mkalist res [- l.pmw, 0, l.pms, l.pmc] pat := cc * t**w * log(t)**m * exp(-t**s) not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) => l := mkalist res [l.pmw, l.pmm, l.pms, l.pmc] pat := cc / t**w * log(t)**m * exp(-t**s) not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) => l := mkalist res [- l.pmw, l.pmm, l.pms, l.pmc] empty() pmintegrate(f, x, a, b) == zero? a and one? whatInfinity b => formula1(f, x, constant(x::F), suchThat(c, freeOf?(#1, x))) "failed" @ \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}