\documentclass{article}
\usepackage{open-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, 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, 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}
<<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 : 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}
<<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}