Initial population.

1 files changed, 782 insertions, 0 deletions

diff --git a/src/algebra/intaf.spad.pamphlet b/src/algebra/intaf.spad.pamphlet
new file mode 100644
index 00000000..23f73b17
--- /dev/null
+++ b/src/algebra/intaf.spad.pamphlet
@@ -0,0 +1,782 @@
+\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}
+<<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}
+<<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"
+      not ((degree g) = 1) => 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) =>
+      ((rec.coef) = 1) =>
+        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"::SY) => prootintegrate(f, x, y)
+          is?(y,  "rootOf"::SY) => 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"::SY) => prootextint(f, x, y, g)
+          is?(y,  "rootOf"::SY) => 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"::SY) => prootlimint(f, x, y, lu)
+          is?(y,  "rootOf"::SY) => 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"::SY) => 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}
+<<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"::SY) => rootintegrate(f, t, y, derivation)
+      is?(y, "rootOf"::SY)  => algintegrate(f, t, y, derivation)
+      FAIL
+
+@
+\section{License}
+<<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.
+@
+<<*>>=
+<<license>>
+
+-- 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
+
+<<package INTG0 GenusZeroIntegration>>
+<<package INTPAF PureAlgebraicIntegration>>
+<<package INTAF AlgebraicIntegration>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}