\documentclass{article} \usepackage{axiom} \begin{document} \title{\$SPAD/src/algebra intaf.spad} \author{Manuel Bronstein} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package INTG0 GenusZeroIntegration} <>= )abbrev package INTG0 GenusZeroIntegration ++ Rationalization of several types of genus 0 integrands; ++ Author: Manuel Bronstein ++ Date Created: 11 October 1988 ++ Date Last Updated: 24 June 1994 ++ Description: ++ This internal package rationalises integrands on curves of the form: ++ \spad{y\^2 = a x\^2 + b x + c} ++ \spad{y\^2 = (a x + b) / (c x + d)} ++ \spad{f(x, y) = 0} where f has degree 1 in x ++ The rationalization is done for integration, limited integration, ++ extended integration and the risch differential equation; GenusZeroIntegration(R, F, L): Exports == Implementation where R: Join(GcdDomain, RetractableTo Integer, OrderedSet, CharacteristicZero, LinearlyExplicitRingOver Integer) F: Join(FunctionSpace R, AlgebraicallyClosedField, TranscendentalFunctionCategory) L: SetCategory SY ==> Symbol Q ==> Fraction Integer K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial F RF ==> Fraction UP UPUP ==> SparseUnivariatePolynomial RF IR ==> IntegrationResult F LOG ==> Record(coeff:F, logand:F) U1 ==> Union(F, "failed") U2 ==> Union(Record(ratpart:F, coeff:F),"failed") U3 ==> Union(Record(mainpart:F, limitedlogs:List LOG), "failed") REC ==> Record(coeff:F, var:List K, val:List F) ODE ==> Record(particular: Union(F, "failed"), basis: List F) LODO==> LinearOrdinaryDifferentialOperator1 RF Exports ==> with palgint0 : (F, K, K, F, UP) -> IR ++ palgint0(f, x, y, d, p) returns the integral of \spad{f(x,y)dx} ++ where y is an algebraic function of x satisfying ++ \spad{d(x)\^2 y(x)\^2 = P(x)}. palgint0 : (F, K, K, K, F, RF) -> IR ++ palgint0(f, x, y, z, t, c) returns the integral of \spad{f(x,y)dx} ++ where y is an algebraic function of x satisfying ++ \spad{f(x,y)dx = c f(t,y) dy}; c and t are rational functions of y. ++ Argument z is a dummy variable not appearing in \spad{f(x,y)}. palgextint0: (F, K, K, F, F, UP) -> U2 ++ palgextint0(f, x, y, g, d, p) returns functions \spad{[h, c]} such ++ that \spad{dh/dx = f(x,y) - c g}, where y is an algebraic function ++ of x satisfying \spad{d(x)\^2 y(x)\^2 = P(x)}, ++ or "failed" if no such functions exist. palgextint0: (F, K, K, F, K, F, RF) -> U2 ++ palgextint0(f, x, y, g, z, t, c) returns functions \spad{[h, d]} such ++ that \spad{dh/dx = f(x,y) - d g}, where y is an algebraic function ++ of x satisfying \spad{f(x,y)dx = c f(t,y) dy}, and c and t are rational ++ functions of y. ++ Argument z is a dummy variable not appearing in \spad{f(x,y)}. ++ The operation returns "failed" if no such functions exist. palglimint0: (F, K, K, List F, F, UP) -> U3 ++ palglimint0(f, x, y, [u1,...,un], d, p) returns functions ++ \spad{[h,[[ci, ui]]]} such that the ui's are among \spad{[u1,...,un]} ++ and \spad{d(h + sum(ci log(ui)))/dx = f(x,y)} if such functions exist, ++ and "failed" otherwise. ++ Argument y is an algebraic function of x satisfying ++ \spad{d(x)\^2y(x)\^2 = P(x)}. palglimint0: (F, K, K, List F, K, F, RF) -> U3 ++ palglimint0(f, x, y, [u1,...,un], z, t, c) returns functions ++ \spad{[h,[[ci, ui]]]} such that the ui's are among \spad{[u1,...,un]} ++ and \spad{d(h + sum(ci log(ui)))/dx = f(x,y)} if such functions exist, ++ and "failed" otherwise. ++ Argument y is an algebraic function of x satisfying ++ \spad{f(x,y)dx = c f(t,y) dy}; c and t are rational functions of y. palgRDE0 : (F, F, K, K, (F, F, SY) -> U1, F, UP) -> U1 ++ palgRDE0(f, g, x, y, foo, d, p) returns a function \spad{z(x,y)} ++ such that \spad{dz/dx + n * df/dx z(x,y) = g(x,y)} if such a z exists, ++ and "failed" otherwise. ++ Argument y is an algebraic function of x satisfying ++ \spad{d(x)\^2y(x)\^2 = P(x)}. ++ Argument foo, called by \spad{foo(a, b, x)}, is a function that solves ++ \spad{du/dx + n * da/dx u(x) = u(x)} ++ for an unknown \spad{u(x)} not involving y. palgRDE0 : (F, F, K, K, (F, F, SY) -> U1, K, F, RF) -> U1 ++ palgRDE0(f, g, x, y, foo, t, c) returns a function \spad{z(x,y)} ++ such that \spad{dz/dx + n * df/dx z(x,y) = g(x,y)} if such a z exists, ++ and "failed" otherwise. ++ Argument y is an algebraic function of x satisfying ++ \spad{f(x,y)dx = c f(t,y) dy}; c and t are rational functions of y. ++ Argument \spad{foo}, called by \spad{foo(a, b, x)}, is a function that ++ solves \spad{du/dx + n * da/dx u(x) = u(x)} ++ for an unknown \spad{u(x)} not involving y. univariate: (F, K, K, UP) -> UPUP ++ univariate(f,k,k,p) \undocumented multivariate: (UPUP, K, F) -> F ++ multivariate(u,k,f) \undocumented lift: (UP, K) -> UPUP ++ lift(u,k) \undocumented if L has LinearOrdinaryDifferentialOperatorCategory F then palgLODE0 : (L, F, K, K, F, UP) -> ODE ++ palgLODE0(op, g, x, y, d, p) returns the solution of \spad{op f = g}. ++ Argument y is an algebraic function of x satisfying ++ \spad{d(x)\^2y(x)\^2 = P(x)}. palgLODE0 : (L, F, K, K, K, F, RF) -> ODE ++ palgLODE0(op,g,x,y,z,t,c) returns the solution of \spad{op f = g} ++ Argument y is an algebraic function of x satisfying ++ \spad{f(x,y)dx = c f(t,y) dy}; c and t are rational functions of y. Implementation ==> add import RationalIntegration(F, UP) import AlgebraicManipulations(R, F) import IntegrationResultFunctions2(RF, F) import ElementaryFunctionStructurePackage(R, F) import SparseUnivariatePolynomialFunctions2(F, RF) import PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) mkRat : (F, REC, List K) -> RF mkRatlx : (F, K, K, F, K, RF) -> RF quadsubst: (K, K, F, UP) -> Record(diff:F, subs:REC, newk:List K) kerdiff : (F, F) -> List K checkroot: (F, List K) -> F univ : (F, List K, K) -> RF dummy := kernel(new()$SY)@K kerdiff(sa, a) == setDifference(kernels sa, kernels a) checkroot(f, l) == (empty? l => f; rootNormalize(f, first l)) univ(c, l, x) == univariate(checkroot(c, l), x) univariate(f, x, y, p) == lift(univariate(f, y, p), x) lift(p, k) == map(univariate(#1, k), p) palgint0(f, x, y, den, radi) == -- y is a square root so write f as f1 y + f0 and integrate separately ff := univariate(f, x, y, minPoly y) f0 := reductum ff pr := quadsubst(x, y, den, radi) map(#1(x::F), integrate(retract(f0)@RF)) + map(#1(pr.diff), integrate mkRat(multivariate(leadingMonomial ff,x,y::F), pr.subs, pr.newk)) -- the algebraic relation is (den * y)**2 = p where p is a * x**2 + b * x + c -- if p is squarefree, then parametrize in the following form: -- u = y - x \sqrt{a} -- x = (u^2 - c) / (b - 2 u \sqrt{a}) = h(u) -- dx = h'(u) du -- y = (u + a h(u)) / den = g(u) -- if a is a perfect square, -- u = (y - \sqrt{c}) / x -- x = (b - 2 u \sqrt{c}) / (u^2 - a) = h(u) -- dx = h'(u) du -- y = (u h(u) + \sqrt{c}) / den = g(u) -- otherwise. -- if p is a square p = a t^2, then we choose only one branch for now: -- u = x -- x = u = h(u) -- dx = du -- y = t \sqrt{a} / den = g(u) -- returns [u(x,y), [h'(u), [x,y], [h(u), g(u)], l] in both cases, -- where l is empty if no new square root was needed, -- l := [k] if k is the new square root kernel that was created. quadsubst(x, y, den, p) == u := dummy::F b := coefficient(p, 1) c := coefficient(p, 0) sa := rootSimp sqrt(a := coefficient(p, 2)) zero?(b * b - 4 * a * c) => -- case where p = a (x + b/(2a))^2 [x::F, [1, [x, y], [u, sa * (u + b / (2*a)) / eval(den,x,u)]], empty()] empty? kerdiff(sa, a) => bm2u := b - 2 * u * sa q := eval(den, x, xx := (u**2 - c) / bm2u) yy := (ua := u + xx * sa) / q [y::F - x::F * sa, [2 * ua / bm2u, [x, y], [xx, yy]], empty()] u2ma:= u**2 - a sc := rootSimp sqrt c q := eval(den, x, xx := (b - 2 * u * sc) / u2ma) yy := (ux := xx * u + sc) / q [(y::F - sc) / x::F, [- 2 * ux / u2ma, [x ,y], [xx, yy]], kerdiff(sc, c)] mkRatlx(f,x,y,t,z,dx) == rat := univariate(eval(f, [x, y], [t, z::F]), z) * dx numer(rat) / denom(rat) mkRat(f, rec, l) == rat:=univariate(checkroot(rec.coeff * eval(f,rec.var,rec.val), l), dummy) numer(rat) / denom(rat) palgint0(f, x, y, z, xx, dx) == map(multivariate(#1, y), integrate mkRatlx(f, x, y, xx, z, dx)) palgextint0(f, x, y, g, z, xx, dx) == map(multivariate(#1, y), extendedint(mkRatlx(f,x,y,xx,z,dx), mkRatlx(g,x,y,xx,z,dx))) palglimint0(f, x, y, lu, z, xx, dx) == map(multivariate(#1, y), limitedint(mkRatlx(f, x, y, xx, z, dx), [mkRatlx(u, x, y, xx, z, dx) for u in lu])) palgRDE0(f, g, x, y, rischde, z, xx, dx) == (u := rischde(eval(f, [x, y], [xx, z::F]), multivariate(dx, z) * eval(g, [x, y], [xx, z::F]), symbolIfCan(z)::SY)) case "failed" => "failed" eval(u::F, z, y::F) -- given p = sum_i a_i(X) Y^i, returns sum_i a_i(x) y^i multivariate(p, x, y) == (map(multivariate(#1, x), p)$SparseUnivariatePolynomialFunctions2(RF, F)) (y) palgextint0(f, x, y, g, den, radi) == pr := quadsubst(x, y, den, radi) map(#1(pr.diff), extendedint(mkRat(f, pr.subs, pr.newk), mkRat(g, pr.subs, pr.newk))) palglimint0(f, x, y, lu, den, radi) == pr := quadsubst(x, y, den, radi) map(#1(pr.diff), limitedint(mkRat(f, pr.subs, pr.newk), [mkRat(u, pr.subs, pr.newk) for u in lu])) palgRDE0(f, g, x, y, rischde, den, radi) == pr := quadsubst(x, y, den, radi) (u := rischde(checkroot(eval(f, pr.subs.var, pr.subs.val), pr.newk), checkroot(pr.subs.coeff * eval(g, pr.subs.var, pr.subs.val), pr.newk), symbolIfCan(dummy)::SY)) case "failed" => "failed" eval(u::F, dummy, pr.diff) if L has LinearOrdinaryDifferentialOperatorCategory F then import RationalLODE(F, UP) palgLODE0(eq, g, x, y, den, radi) == pr := quadsubst(x, y, den, radi) d := monomial(univ(inv(pr.subs.coeff), pr.newk, dummy), 1)$LODO di:LODO := 1 -- will accumulate the powers of d op:LODO := 0 -- will accumulate the new LODO for i in 0..degree eq repeat op := op + univ(eval(coefficient(eq, i), pr.subs.var, pr.subs.val), pr.newk, dummy) * di di := d * di rec := ratDsolve(op,univ(eval(g,pr.subs.var,pr.subs.val),pr.newk,dummy)) bas:List(F) := [b(pr.diff) for b in rec.basis] rec.particular case "failed" => ["failed", bas] [((rec.particular)::RF) (pr.diff), bas] palgLODE0(eq, g, x, y, kz, xx, dx) == d := monomial(univariate(inv multivariate(dx, kz), kz), 1)$LODO di:LODO := 1 -- will accumulate the powers of d op:LODO := 0 -- will accumulate the new LODO lk:List(K) := [x, y] lv:List(F) := [xx, kz::F] for i in 0..degree eq repeat op := op + univariate(eval(coefficient(eq, i), lk, lv), kz) * di di := d * di rec := ratDsolve(op, univariate(eval(g, lk, lv), kz)) bas:List(F) := [multivariate(b, y) for b in rec.basis] rec.particular case "failed" => ["failed", bas] [multivariate((rec.particular)::RF, y), bas] @ \section{package INTPAF PureAlgebraicIntegration} <>= )abbrev package INTPAF PureAlgebraicIntegration ++ Integration of pure algebraic functions; ++ Author: Manuel Bronstein ++ Date Created: 27 May 1988 ++ Date Last Updated: 24 June 1994 ++ Description: ++ This package provides functions for integration, limited integration, ++ extended integration and the risch differential equation for ++ pure algebraic integrands; PureAlgebraicIntegration(R, F, L): Exports == Implementation where R: Join(GcdDomain,RetractableTo Integer,OrderedSet, CharacteristicZero, LinearlyExplicitRingOver Integer) F: Join(FunctionSpace R, AlgebraicallyClosedField, TranscendentalFunctionCategory) L: SetCategory SY ==> Symbol N ==> NonNegativeInteger K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial F RF ==> Fraction UP UPUP==> SparseUnivariatePolynomial RF IR ==> IntegrationResult F IR2 ==> IntegrationResultFunctions2(curve, F) ALG ==> AlgebraicIntegrate(R, F, UP, UPUP, curve) LDALG ==> LinearOrdinaryDifferentialOperator1 curve RDALG ==> PureAlgebraicLODE(F, UP, UPUP, curve) LOG ==> Record(coeff:F, logand:F) REC ==> Record(particular:U1, basis:List F) CND ==> Record(left:UP, right:UP) CHV ==> Record(int:UPUP, left:UP, right:UP, den:RF, deg:N) U1 ==> Union(F, "failed") U2 ==> Union(Record(ratpart:F, coeff:F),"failed") U3 ==> Union(Record(mainpart:F, limitedlogs:List LOG), "failed") FAIL==> error "failed - cannot handle that integrand" Exports ==> with palgint : (F, K, K) -> IR ++ palgint(f, x, y) returns the integral of \spad{f(x,y)dx} ++ where y is an algebraic function of x. palgextint: (F, K, K, F) -> U2 ++ palgextint(f, x, y, g) returns functions \spad{[h, c]} such that ++ \spad{dh/dx = f(x,y) - c g}, where y is an algebraic function of x; ++ returns "failed" if no such functions exist. palglimint: (F, K, K, List F) -> U3 ++ palglimint(f, x, y, [u1,...,un]) returns functions ++ \spad{[h,[[ci, ui]]]} such that the ui's are among \spad{[u1,...,un]} ++ and \spad{d(h + sum(ci log(ui)))/dx = f(x,y)} if such functions exist, ++ "failed" otherwise; ++ y is an algebraic function of x. palgRDE : (F, F, F, K, K, (F, F, SY) -> U1) -> U1 ++ palgRDE(nfp, f, g, x, y, foo) returns a function \spad{z(x,y)} ++ such that \spad{dz/dx + n * df/dx z(x,y) = g(x,y)} if such a z exists, ++ "failed" otherwise; ++ y is an algebraic function of x; ++ \spad{foo(a, b, x)} is a function that solves ++ \spad{du/dx + n * da/dx u(x) = u(x)} ++ for an unknown \spad{u(x)} not involving y. ++ \spad{nfp} is \spad{n * df/dx}. if L has LinearOrdinaryDifferentialOperatorCategory F then palgLODE: (L, F, K, K, SY) -> REC ++ palgLODE(op, g, kx, y, x) returns the solution of \spad{op f = g}. ++ y is an algebraic function of x. Implementation ==> add import IntegrationTools(R, F) import RationalIntegration(F, UP) import GenusZeroIntegration(R, F, L) import ChangeOfVariable(F, UP, UPUP) import IntegrationResultFunctions2(F, F) import IntegrationResultFunctions2(RF, F) import SparseUnivariatePolynomialFunctions2(F, RF) import UnivariatePolynomialCommonDenominator(UP, RF, UPUP) import PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) quadIfCan : (K, K) -> Union(Record(coef:F, poly:UP), "failed") linearInXIfCan : (K, K) -> Union(Record(xsub:F, dxsub:RF), "failed") prootintegrate : (F, K, K) -> IR prootintegrate1: (UPUP, K, K, UPUP) -> IR prootextint : (F, K, K, F) -> U2 prootlimint : (F, K, K, List F) -> U3 prootRDE : (F, F, F, K, K, (F, F, SY) -> U1) -> U1 palgRDE1 : (F, F, K, K) -> U1 palgLODE1 : (List F, F, K, K, SY) -> REC palgintegrate : (F, K, K) -> IR palgext : (F, K, K, F) -> U2 palglim : (F, K, K, List F) -> U3 UPUP2F1 : (UPUP, RF, RF, K, K) -> F UPUP2F0 : (UPUP, K, K) -> F RF2UPUP : (RF, UPUP) -> UPUP algaddx : (IR, F) -> IR chvarIfCan : (UPUP, RF, UP, RF) -> Union(UPUP, "failed") changeVarIfCan : (UPUP, RF, N) -> Union(CHV, "failed") rationalInt : (UPUP, N, UP) -> IntegrationResult RF chv : (UPUP, N, F, F) -> RF chv0 : (UPUP, N, F, F) -> F candidates : UP -> List CND dummy := new()$SY dumk := kernel(dummy)@K UPUP2F1(p, t, cf, kx, k) == UPUP2F0(eval(p, t, cf), kx, k) UPUP2F0(p, kx, k) == multivariate(p, kx, k::F) chv(f, n, a, b) == univariate(chv0(f, n, a, b), dumk) RF2UPUP(f, modulus) == bc := extendedEuclidean(map(#1::UP::RF, denom f), modulus, 1)::Record(coef1:UPUP, coef2:UPUP) (map(#1::UP::RF, numer f) * bc.coef1) rem modulus -- returns "failed", or (xx, c) such that f(x, y)dx = f(xx, y) c dy -- if p(x, y) = 0 is linear in x linearInXIfCan(x, y) == a := b := 0$UP p := clearDenominator lift(minPoly y, x) while p ~= 0 repeat degree(q := numer leadingCoefficient p) > 1 => return "failed" a := a + monomial(coefficient(q, 1), d := degree p) b := b - monomial(coefficient(q, 0), d) p := reductum p xx:RF := b / a [xx(dumk::F), differentiate(xx, differentiate)] -- return Int(f(x,y)dx) where y is an n^th root of a rational function in x prootintegrate(f, x, y) == modulus := lift(p := minPoly y, x) rf := reductum(ff := univariate(f, x, y, p)) ((r := retractIfCan(rf)@Union(RF,"failed")) case RF) and rf ~= 0 => -- in this case, ff := lc(ff) y^i + r so we integrate both terms -- separately to gain time map(#1(x::F), integrate(r::RF)) + prootintegrate1(leadingMonomial ff, x, y, modulus) prootintegrate1(ff, x, y, modulus) prootintegrate1(ff, x, y, modulus) == chv:CHV r := radPoly(modulus)::Record(radicand:RF, deg:N) (uu := changeVarIfCan(ff, r.radicand, r.deg)) case CHV => chv := uu::CHV newalg := nthRoot((chv.left)(dumk::F), chv.deg) kz := retract(numer newalg)@K newf := multivariate(chv.int, ku := dumk, newalg) vu := (chv.right)(x::F) vz := (chv.den)(x::F) * (y::F) * denom(newalg)::F map(eval(#1, [ku, kz], [vu, vz]), palgint(newf, ku, kz)) cv := chvar(ff, modulus) r := radPoly(cv.poly)::Record(radicand:RF, deg:N) qprime := differentiate(q := retract(r.radicand)@UP)::RF not zero? qprime and ((u := chvarIfCan(cv.func, 1, q, inv qprime)) case UPUP) => m := monomial(1, r.deg)$UPUP - q::RF::UPUP map(UPUP2F1(RF2UPUP(#1, m), cv.c1, cv.c2, x, y), rationalInt(u::UPUP, r.deg, monomial(1, 1))) curve := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg) algaddx(map(UPUP2F1(lift #1, cv.c1, cv.c2, x, y), palgintegrate(reduce(cv.func), differentiate$UP)$ALG)$IR2, x::F) -- Do the rationalizing change of variable -- Int(f(x, y) dx) --> Int(n u^(n-1) f((u^n - b)/a, u) / a du) where -- u^n = y^n = g(x) = a x + b -- returns the integral as an integral of a rational function in u rationalInt(f, n, g) == not one? degree g => error "rationalInt: radicand must be linear" a := leadingCoefficient g integrate(n * monomial(inv a, (n-1)::N)$UP * chv(f, n, a, leadingCoefficient reductum g)) -- Do the rationalizing change of variable f(x,y) --> f((u^n - b)/a, u) where -- u = y = (a x + b)^(1/n). -- Returns f((u^n - b)/a,u) as an element of F chv0(f, n, a, b) == d := dumk::F (f (d::UP::RF)) ((d ** n - b) / a) -- candidates(p) returns a list of pairs [g, u] such that p(x) = g(u(x)), -- those u's are candidates for change of variables -- currently uses a dumb heuristic where the candidates u's are p itself -- and all the powers x^2, x^3, ..., x^{deg(p)}, -- will use polynomial decomposition in smarter days MB 8/93 candidates p == l:List(CND) := empty() ground? p => l for i in 2..degree p repeat if (u := composite(p, xi := monomial(1, i))) case UP then l := concat([u::UP, xi], l) concat([monomial(1, 1), p], l) -- checks whether Int(p(x, y) dx) can be rewritten as -- Int(r(u, z) du) where u is some polynomial of x, -- z = d y for some polynomial d, and z^m = g(u) -- returns either [r(u, z), g, u, d, m] or "failed" -- we have y^n = radi changeVarIfCan(p, radi, n) == rec := rootPoly(radi, n) for cnd in candidates(rec.radicand) repeat (u := chvarIfCan(p, rec.coef, cnd.right, inv(differentiate(cnd.right)::RF))) case UPUP => return [u::UPUP, cnd.left, cnd.right, rec.coef, rec.exponent] "failed" -- checks whether Int(p(x, y) dx) can be rewritten as -- Int(r(u, z) du) where u is some polynomial of x and z = d y -- we have y^n = a(x)/d(x) -- returns either "failed" or r(u, z) chvarIfCan(p, d, u, u1) == ans:UPUP := 0 while p ~= 0 repeat (v := composite(u1 * leadingCoefficient(p) / d ** degree(p), u)) case "failed" => return "failed" ans := ans + monomial(v::RF, degree p) p := reductum p ans algaddx(i, xx) == elem? i => i mkAnswer(ratpart i, logpart i, [[- ne.integrand / (xx**2), xx] for ne in notelem i]) prootRDE(nfp, f, g, x, k, rde) == modulus := lift(p := minPoly k, x) r := radPoly(modulus)::Record(radicand:RF, deg:N) rec := rootPoly(r.radicand, r.deg) dqdx := inv(differentiate(q := rec.radicand)::RF) ((uf := chvarIfCan(ff := univariate(f,x,k,p),rec.coef,q,1)) case UPUP) and ((ug:=chvarIfCan(gg:=univariate(g,x,k,p),rec.coef,q,dqdx)) case UPUP) => (u := rde(chv0(uf::UPUP, rec.exponent, 1, 0), rec.exponent * (dumk::F) ** (rec.exponent * (rec.exponent - 1)) * chv0(ug::UPUP, rec.exponent, 1, 0), symbolIfCan(dumk)::SY)) case "failed" => "failed" eval(u::F, dumk, k::F) one?(rec.coef) => curve := RadicalFunctionField(F, UP, UPUP, q::RF, rec.exponent) rc := algDsolve(D()$LDALG + reduce(univariate(nfp, x, k, p))::LDALG, reduce univariate(g, x, k, p))$RDALG rc.particular case "failed" => "failed" UPUP2F0(lift((rc.particular)::curve), x, k) palgRDE1(nfp, g, x, k) prootlimint(f, x, k, lu) == modulus := lift(p := minPoly k, x) r := radPoly(modulus)::Record(radicand:RF, deg:N) rec := rootPoly(r.radicand, r.deg) dqdx := inv(differentiate(q := rec.radicand)::RF) (uf := chvarIfCan(ff := univariate(f,x,k,p),rec.coef,q,dqdx)) case UPUP => l := empty()$List(RF) n := rec.exponent * monomial(1, (rec.exponent - 1)::N)$UP for u in lu repeat if ((v:=chvarIfCan(uu:=univariate(u,x,k,p),rec.coef,q,dqdx))case UPUP) then l := concat(n * chv(v::UPUP,rec.exponent, 1, 0), l) else FAIL m := monomial(1, rec.exponent)$UPUP - q::RF::UPUP map(UPUP2F0(RF2UPUP(#1,m), x, k), limitedint(n * chv(uf::UPUP, rec.exponent, 1, 0), reverse_! l)) cv := chvar(ff, modulus) r := radPoly(cv.poly)::Record(radicand:RF, deg:N) dqdx := inv(differentiate(q := retract(r.radicand)@UP)::RF) curve := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg) (ui := palginfieldint(reduce(cv.func), differentiate$UP)$ALG) case "failed" => FAIL [UPUP2F1(lift(ui::curve), cv.c1, cv.c2, x, k), empty()] prootextint(f, x, k, g) == modulus := lift(p := minPoly k, x) r := radPoly(modulus)::Record(radicand:RF, deg:N) rec := rootPoly(r.radicand, r.deg) dqdx := inv(differentiate(q := rec.radicand)::RF) ((uf:=chvarIfCan(ff:=univariate(f,x,k,p),rec.coef,q,dqdx)) case UPUP) and ((ug:=chvarIfCan(gg:=univariate(g,x,k,p),rec.coef,q,dqdx)) case UPUP) => m := monomial(1, rec.exponent)$UPUP - q::RF::UPUP n := rec.exponent * monomial(1, (rec.exponent - 1)::N)$UP map(UPUP2F0(RF2UPUP(#1,m), x, k), extendedint(n * chv(uf::UPUP, rec.exponent, 1, 0), n * chv(ug::UPUP, rec.exponent, 1, 0))) cv := chvar(ff, modulus) r := radPoly(cv.poly)::Record(radicand:RF, deg:N) dqdx := inv(differentiate(q := retract(r.radicand)@UP)::RF) curve := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg) (u := palginfieldint(reduce(cv.func), differentiate$UP)$ALG) case "failed" => FAIL [UPUP2F1(lift(u::curve), cv.c1, cv.c2, x, k), 0] palgRDE1(nfp, g, x, y) == palgLODE1([nfp, 1], g, x, y, symbolIfCan(x)::SY).particular palgLODE1(eq, g, kx, y, x) == modulus:= lift(p := minPoly y, kx) curve := AlgebraicFunctionField(F, UP, UPUP, modulus) neq:LDALG := 0 for f in eq for i in 0.. repeat neq := neq + monomial(reduce univariate(f, kx, y, p), i) empty? remove_!(y, remove_!(kx, varselect(kernels g, x))) => rec := algDsolve(neq, reduce univariate(g, kx, y, p))$RDALG bas:List(F) := [UPUP2F0(lift h, kx, y) for h in rec.basis] rec.particular case "failed" => ["failed", bas] [UPUP2F0(lift((rec.particular)::curve), kx, y), bas] rec := algDsolve(neq, 0) ["failed", [UPUP2F0(lift h, kx, y) for h in rec.basis]] palgintegrate(f, x, k) == modulus:= lift(p := minPoly k, x) cv := chvar(univariate(f, x, k, p), modulus) curve := AlgebraicFunctionField(F, UP, UPUP, cv.poly) knownInfBasis(cv.deg) algaddx(map(UPUP2F1(lift #1, cv.c1, cv.c2, x, k), palgintegrate(reduce(cv.func), differentiate$UP)$ALG)$IR2, x::F) palglim(f, x, k, lu) == modulus:= lift(p := minPoly k, x) cv := chvar(univariate(f, x, k, p), modulus) curve := AlgebraicFunctionField(F, UP, UPUP, cv.poly) knownInfBasis(cv.deg) (u := palginfieldint(reduce(cv.func), differentiate$UP)$ALG) case "failed" => FAIL [UPUP2F1(lift(u::curve), cv.c1, cv.c2, x, k), empty()] palgext(f, x, k, g) == modulus:= lift(p := minPoly k, x) cv := chvar(univariate(f, x, k, p), modulus) curve := AlgebraicFunctionField(F, UP, UPUP, cv.poly) knownInfBasis(cv.deg) (u := palginfieldint(reduce(cv.func), differentiate$UP)$ALG) case "failed" => FAIL [UPUP2F1(lift(u::curve), cv.c1, cv.c2, x, k), 0] palgint(f, x, y) == (v := linearInXIfCan(x, y)) case "failed" => (u := quadIfCan(x, y)) case "failed" => is?(y, 'nthRoot) => prootintegrate(f, x, y) is?(y, 'rootOf) => palgintegrate(f, x, y) FAIL palgint0(f, x, y, u.coef, u.poly) palgint0(f, x, y, dumk, v.xsub, v.dxsub) palgextint(f, x, y, g) == (v := linearInXIfCan(x, y)) case "failed" => (u := quadIfCan(x, y)) case "failed" => is?(y, 'nthRoot) => prootextint(f, x, y, g) is?(y, 'rootOf) => palgext(f, x, y, g) FAIL palgextint0(f, x, y, g, u.coef, u.poly) palgextint0(f, x, y, g, dumk, v.xsub, v.dxsub) palglimint(f, x, y, lu) == (v := linearInXIfCan(x, y)) case "failed" => (u := quadIfCan(x, y)) case "failed" => is?(y, 'nthRoot) => prootlimint(f, x, y, lu) is?(y, 'rootOf) => palglim(f, x, y, lu) FAIL palglimint0(f, x, y, lu, u.coef, u.poly) palglimint0(f, x, y, lu, dumk, v.xsub, v.dxsub) palgRDE(nfp, f, g, x, y, rde) == (v := linearInXIfCan(x, y)) case "failed" => (u := quadIfCan(x, y)) case "failed" => is?(y, 'nthRoot) => prootRDE(nfp, f, g, x, y, rde) palgRDE1(nfp, g, x, y) palgRDE0(f, g, x, y, rde, u.coef, u.poly) palgRDE0(f, g, x, y, rde, dumk, v.xsub, v.dxsub) -- returns "failed", or (d, P) such that (dy)**2 = P(x) -- and degree(P) = 2 quadIfCan(x, y) == (degree(p := minPoly y) = 2) and zero?(coefficient(p, 1)) => d := denom(ff := univariate(- coefficient(p, 0) / coefficient(p, 2), x)) degree(radi := d * numer ff) = 2 => [d(x::F), radi] "failed" "failed" if L has LinearOrdinaryDifferentialOperatorCategory F then palgLODE(eq, g, kx, y, x) == (v := linearInXIfCan(kx, y)) case "failed" => (u := quadIfCan(kx, y)) case "failed" => palgLODE1([coefficient(eq, i) for i in 0..degree eq], g, kx, y, x) palgLODE0(eq, g, kx, y, u.coef, u.poly) palgLODE0(eq, g, kx, y, dumk, v.xsub, v.dxsub) @ \section{package INTAF AlgebraicIntegration} <>= )abbrev package INTAF AlgebraicIntegration ++ Mixed algebraic integration; ++ Author: Manuel Bronstein ++ Date Created: 12 October 1988 ++ Date Last Updated: 4 June 1988 ++ Description: ++ This package provides functions for the integration of ++ algebraic integrands over transcendental functions; AlgebraicIntegration(R, F): Exports == Implementation where R : Join(OrderedSet, IntegralDomain) F : Join(AlgebraicallyClosedField, FunctionSpace R) SY ==> Symbol N ==> NonNegativeInteger K ==> Kernel F P ==> SparseMultivariatePolynomial(R, K) UP ==> SparseUnivariatePolynomial F RF ==> Fraction UP UPUP==> SparseUnivariatePolynomial RF IR ==> IntegrationResult F IR2 ==> IntegrationResultFunctions2(curve, F) ALG ==> AlgebraicIntegrate(R, F, UP, UPUP, curve) FAIL==> error "failed - cannot handle that integrand" Exports ==> with algint: (F, K, K, UP -> UP) -> IR ++ algint(f, x, y, d) returns the integral of \spad{f(x,y)dx} ++ where y is an algebraic function of x; ++ d is the derivation to use on \spad{k[x]}. Implementation ==> add import ChangeOfVariable(F, UP, UPUP) import PolynomialCategoryQuotientFunctions(IndexedExponents K, K, R, P, F) rootintegrate: (F, K, K, UP -> UP) -> IR algintegrate : (F, K, K, UP -> UP) -> IR UPUP2F : (UPUP, RF, K, K) -> F F2UPUP : (F, K, K, UP) -> UPUP UP2UPUP : (UP, K) -> UPUP F2UPUP(f, kx, k, p) == UP2UPUP(univariate(f, k, p), kx) rootintegrate(f, t, k, derivation) == r1 := mkIntegral(modulus := UP2UPUP(p := minPoly k, t)) f1 := F2UPUP(f, t, k, p) monomial(inv(r1.coef), 1) r := radPoly(r1.poly)::Record(radicand:RF, deg:N) q := retract(r.radicand) curve := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg) map(UPUP2F(lift #1, r1.coef, t, k), algintegrate(reduce f1, derivation)$ALG)$IR2 algintegrate(f, t, k, derivation) == r1 := mkIntegral(modulus := UP2UPUP(p := minPoly k, t)) f1 := F2UPUP(f, t, k, p) monomial(inv(r1.coef), 1) modulus:= UP2UPUP(p := minPoly k, t) curve := AlgebraicFunctionField(F, UP, UPUP, r1.poly) map(UPUP2F(lift #1, r1.coef, t, k), algintegrate(reduce f1, derivation)$ALG)$IR2 UP2UPUP(p, k) == map(univariate(#1,k),p)$SparseUnivariatePolynomialFunctions2(F,RF) UPUP2F(p, cf, t, k) == map(multivariate(#1, t), p)$SparseUnivariatePolynomialFunctions2(RF, F) (multivariate(cf, t) * k::F) algint(f, t, y, derivation) == is?(y, 'nthRoot) => rootintegrate(f, t, y, derivation) is?(y, 'rootOf) => algintegrate(f, t, y, derivation) FAIL @ \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 intef irexpand integrat <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}