path: root/src/algebra/oderf.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra oderf.spad}
\author{Manuel Bronstein}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package BALFACT BalancedFactorisation}
<<package BALFACT BalancedFactorisation>>=
)abbrev package BALFACT BalancedFactorisation
++ Author: Manuel Bronstein
++ Date Created: 1 March 1991
++ Date Last Updated: 11 October 1991
++ Description: This package provides balanced factorisations of polynomials.
BalancedFactorisation(R, UP): Exports == Implementation where
  R  : Join(GcdDomain, CharacteristicZero)
  UP : UnivariatePolynomialCategory R

  Exports ==> with
    balancedFactorisation: (UP, UP) -> Factored UP
      ++ balancedFactorisation(a, b) returns
      ++ a factorisation \spad{a = p1^e1 ... pm^em} such that each
      ++ \spad{pi} is balanced with respect to b.
    balancedFactorisation: (UP, List UP) -> Factored UP
      ++ balancedFactorisation(a, [b1,...,bn]) returns
      ++ a factorisation \spad{a = p1^e1 ... pm^em} such that each
      ++ pi is balanced with respect to \spad{[b1,...,bm]}.

  Implementation ==> add
    balSqfr : (UP, Integer, List UP) -> Factored UP
    balSqfr1: (UP, Integer,      UP) -> Factored UP

    balancedFactorisation(a:UP, b:UP) == balancedFactorisation(a, [b])

    balSqfr1(a, n, b) ==
      g := gcd(a, b)
      fa := sqfrFactor((a exquo g)::UP, n)
      ground? g => fa
      fa * balSqfr1(g, n, (b exquo (g ** order(b, g)))::UP)

    balSqfr(a, n, l) ==
      b := first l
      empty? rest l => balSqfr1(a, n, b)
      */[balSqfr1(f.factor, n, b) for f in factors balSqfr(a,n,rest l)]

    balancedFactorisation(a:UP, l:List UP) ==
      empty?(ll := select(#1 ~= 0, l)) =>
        error "balancedFactorisation: 2nd argument is empty or all 0"
      sa := squareFree a
      unit(sa) * */[balSqfr(f.factor,f.exponent,ll) for f in factors sa])

@
\section{package BOUNDZRO BoundIntegerRoots}
<<package BOUNDZRO BoundIntegerRoots>>=
)abbrev package BOUNDZRO BoundIntegerRoots
++ Author: Manuel Bronstein
++ Date Created: 11 March 1991
++ Date Last Updated: 18 November 1991
++ Description:
++   \spadtype{BoundIntegerRoots} provides functions to
++   find lower bounds on the integer roots of a polynomial.
BoundIntegerRoots(F, UP): Exports == Implementation where
  F  : Join(Field, RetractableTo Fraction Integer)
  UP : UnivariatePolynomialCategory F

  Z   ==> Integer
  Q   ==> Fraction Z
  K   ==> Kernel F
  UPQ ==> SparseUnivariatePolynomial Q
  ALGOP ==> "%alg"

  Exports ==> with
    integerBound: UP -> Z
      ++ integerBound(p) returns a lower bound on the negative integer
      ++ roots of p, and 0 if p has no negative integer roots.

  Implementation ==> add
    import RationalFactorize(UPQ)
    import UnivariatePolynomialCategoryFunctions2(F, UP, Q, UPQ)

    qbound : (UP, UPQ) -> Z
    zroot1 : UP -> Z
    qzroot1: UPQ -> Z
    negint : Q -> Z

-- returns 0 if p has no integer root < 0, its negative integer root otherwise
    qzroot1 p == negint(- leadingCoefficient(reductum p) / leadingCoefficient p)

-- returns 0 if p has no integer root < 0, its negative integer root otherwise
    zroot1 p ==
      z := - leadingCoefficient(reductum p) / leadingCoefficient p
      (r := retractIfCan(z)@Union(Q, "failed")) case Q => negint(r::Q)
      0

-- returns 0 if r is not a negative integer, r otherwise
    negint r ==
      ((u := retractIfCan(r)@Union(Z, "failed")) case Z) and (u::Z < 0) => u::Z
      0

    if F has ExpressionSpace then
      bringDown: F -> Q

-- the random substitution used by bringDown is NOT always a ring-homorphism
-- (because of potential algebraic kernels), but is ALWAYS a Z-linear map.
-- this guarantees that bringing down the coefficients of (x + n) q(x) for an
-- integer n yields a polynomial h(x) which is divisible by x + n
-- the only problem is that evaluating with random numbers can cause a
-- division by 0. We should really be able to trap this error later and
-- reevaluate with a new set of random numbers    MB 11/91
      bringDown f ==
        t := tower f
        retract eval(f, t, [random()$Q :: F for k in t])

      integerBound p ==
        one? degree p => zroot1 p
        q1 := map(bringDown, p)
        q2 := map(bringDown, p)
        qbound(p, gcd(q1, q2))

    else
      integerBound p ==
        one? degree p => zroot1 p
        qbound(p, map(retract(#1)@Q, p))

-- we can probably do better here (i.e. without factoring)
    qbound(p, q) ==
      bound:Z := 0
      for rec in factors factor q repeat
        if one?(degree(rec.factor)) and ((r := qzroot1(rec.factor)) < bound)
           and zero? p(r::Q::F) then bound := r
      bound

@
\section{package ODEPRIM PrimitiveRatDE}
<<package ODEPRIM PrimitiveRatDE>>=
)abbrev package ODEPRIM PrimitiveRatDE
++ Author: Manuel Bronstein
++ Date Created: 1 March 1991
++ Date Last Updated: 1 February 1994
++ Description:
++  \spad{PrimitiveRatDE} provides functions for in-field solutions of linear
++   ordinary differential equations, in the transcendental case.
++   The derivation to use is given by the parameter \spad{L}.
PrimitiveRatDE(F, UP, L, LQ): Exports == Implementation where
  F  : Join(Field, CharacteristicZero, RetractableTo Fraction Integer)
  UP : UnivariatePolynomialCategory F
  L  : LinearOrdinaryDifferentialOperatorCategory UP
  LQ : LinearOrdinaryDifferentialOperatorCategory Fraction UP

  N   ==> NonNegativeInteger
  Z   ==> Integer
  RF  ==> Fraction UP
  UP2 ==> SparseUnivariatePolynomial UP
  REC ==> Record(center:UP, equation:UP)

  Exports ==> with
    denomLODE: (L, RF) -> Union(UP, "failed")
      ++ denomLODE(op, g) returns a polynomial d such that
      ++ any rational solution of \spad{op y = g}
      ++ is of the form \spad{p/d} for some polynomial p, and
      ++ "failed", if the equation has no rational solution.
    denomLODE: (L, List RF) -> UP
      ++ denomLODE(op, [g1,...,gm]) returns a polynomial
      ++ d such that any rational solution of \spad{op y = c1 g1 + ... + cm gm}
      ++ is of the form \spad{p/d} for some polynomial p.
    indicialEquations: L -> List REC
      ++ indicialEquations op returns \spad{[[d1,e1],...,[dq,eq]]} where
      ++ the \spad{d_i}'s are the affine singularities of \spad{op},
      ++ and the \spad{e_i}'s are the indicial equations at each \spad{d_i}.
    indicialEquations: (L, UP) -> List REC
      ++ indicialEquations(op, p) returns \spad{[[d1,e1],...,[dq,eq]]} where
      ++ the \spad{d_i}'s are the affine singularities of \spad{op}
      ++ above the roots of \spad{p},
      ++ and the \spad{e_i}'s are the indicial equations at each \spad{d_i}.
    indicialEquation: (L, F) -> UP
      ++ indicialEquation(op, a) returns the indicial equation of \spad{op}
      ++ at \spad{a}.
    indicialEquations: LQ -> List REC
      ++ indicialEquations op returns \spad{[[d1,e1],...,[dq,eq]]} where
      ++ the \spad{d_i}'s are the affine singularities of \spad{op},
      ++ and the \spad{e_i}'s are the indicial equations at each \spad{d_i}.
    indicialEquations: (LQ, UP) -> List REC
      ++ indicialEquations(op, p) returns \spad{[[d1,e1],...,[dq,eq]]} where
      ++ the \spad{d_i}'s are the affine singularities of \spad{op}
      ++ above the roots of \spad{p},
      ++ and the \spad{e_i}'s are the indicial equations at each \spad{d_i}.
    indicialEquation: (LQ, F) -> UP
      ++ indicialEquation(op, a) returns the indicial equation of \spad{op}
      ++ at \spad{a}.
    splitDenominator: (LQ, List RF) -> Record(eq:L, rh:List RF)
      ++ splitDenominator(op, [g1,...,gm]) returns \spad{op0, [h1,...,hm]}
      ++ such that the equations \spad{op y = c1 g1 + ... + cm gm} and
      ++ \spad{op0 y = c1 h1 + ... + cm hm} have the same solutions.

  Implementation ==> add
    import BoundIntegerRoots(F, UP)
    import BalancedFactorisation(F, UP)
    import InnerCommonDenominator(UP, RF, List UP, List RF)
    import UnivariatePolynomialCategoryFunctions2(F, UP, UP, UP2)

    tau          : (UP, UP, UP, N) -> UP
    NPbound      : (UP, L, UP) -> N
    hdenom       : (L, UP, UP) -> UP
    denom0       : (Z, L, UP, UP, UP) -> UP
    indicialEq   : (UP, List N, List UP) -> UP
    separateZeros: (UP, UP) -> UP
    UPfact       : N -> UP
    UP2UP2       : UP -> UP2
    indeq        : (UP, L) -> UP
    NPmulambda   : (UP, L) -> Record(mu:Z, lambda:List N, func:List UP)

    diff := D()$L

    UP2UP2 p                    == map(#1::UP, p)
    indicialEquations(op:L)     == indicialEquations(op, leadingCoefficient op)
    indicialEquation(op:L, a:F) == indeq(monomial(1, 1) - a::UP, op)

    splitDenominator(op, lg) ==
      cd := splitDenominator coefficients op
      f  := cd.den / gcd(cd.num)
      l:L := 0
      while op ~= 0 repeat
          l  := l + monomial(retract(f * leadingCoefficient op), degree op)
          op := reductum op
      [l, [f * g for g in lg]]

    tau(p, pp, q, n) ==
      ((pp ** n) * ((q exquo (p ** order(q, p)))::UP)) rem p

    indicialEquations(op:LQ) ==
      indicialEquations(splitDenominator(op, empty()).eq)

    indicialEquations(op:LQ, p:UP) ==
      indicialEquations(splitDenominator(op, empty()).eq, p)

    indicialEquation(op:LQ, a:F) ==
      indeq(monomial(1, 1) - a::UP, splitDenominator(op, empty()).eq)

-- returns z(z-1)...(z-(n-1))
    UPfact n ==
      zero? n => 1
      z := monomial(1, 1)$UP
      */[z - i::F::UP for i in 0..(n-1)::N]

    indicialEq(c, lamb, lf) ==
      cp := diff c
      cc := UP2UP2 c
      s:UP2 := 0
      for i in lamb for f in lf repeat
        s := s + (UPfact i) * UP2UP2 tau(c, cp, f, i)
      primitivePart resultant(cc, s)

    NPmulambda(c, l) ==
      lamb:List(N) := [d := degree l]
      lf:List(UP) := [a := leadingCoefficient l]
      mup := d::Z - order(a, c)
      while (l := reductum l) ~= 0 repeat
          a := leadingCoefficient l
          if (m := (d := degree l)::Z - order(a, c)) > mup then
              mup := m
              lamb := [d]
              lf := [a]
          else if (m = mup) then
              lamb := concat(d, lamb)
              lf := concat(a, lf)
      [mup, lamb, lf]

-- e = 0 means homogeneous equation
    NPbound(c, l, e) ==
      rec := NPmulambda(c, l)
      n := max(0, - integerBound indicialEq(c, rec.lambda, rec.func))
      zero? e => n::N
      max(n, order(e, c)::Z - rec.mu)::N

    hdenom(l, d, e) ==
      */[dd.factor ** NPbound(dd.factor, l, e)
                    for dd in factors balancedFactorisation(d, coefficients l)]

    denom0(n, l, d, e, h) ==
      hdenom(l, d, e) * */[hh.factor ** max(0, order(e, hh.factor) - n)::N
                                 for hh in factors balancedFactorisation(h, e)]

-- returns a polynomials whose zeros are the zeros of e which are not
-- zeros of d
    separateZeros(d, e) ==
      ((g := squareFreePart e) exquo gcd(g, squareFreePart d))::UP

    indeq(c, l) ==
      rec := NPmulambda(c, l)
      indicialEq(c, rec.lambda, rec.func)

    indicialEquations(op:L, p:UP) ==
      [[dd.factor, indeq(dd.factor, op)]
                   for dd in factors balancedFactorisation(p, coefficients op)]

-- cannot return "failed" in the homogeneous case
    denomLODE(l:L, g:RF) ==
      d := leadingCoefficient l
      zero? g => hdenom(l, d, 0)
      h := separateZeros(d, e := denom g)
      n := degree l
      (e exquo (h**(n + 1))) case "failed" => "failed"
      denom0(n, l, d, e, h)

    denomLODE(l:L, lg:List RF) ==
      empty? lg => denomLODE(l, 0)::UP
      d := leadingCoefficient l
      h := separateZeros(d, e := "lcm"/[denom g for g in lg])
      denom0(degree l, l, d, e, h)

@
\section{package UTSODETL UTSodetools}
<<package UTSODETL UTSodetools>>=
)abbrev package UTSODETL UTSodetools
++ Author: Manuel Bronstein
++ Date Created: 31 January 1994
++ Date Last Updated: 3 February 1994
++ Description:
++  \spad{RUTSodetools} provides tools to interface with the series
++   ODE solver when presented with linear ODEs.
UTSodetools(F, UP, L, UTS): Exports == Implementation where
  F  : Ring
  UP : UnivariatePolynomialCategory F
  L  : LinearOrdinaryDifferentialOperatorCategory UP
  UTS: UnivariateTaylorSeriesCategory F

  Exports ==> with
      UP2UTS:   UP -> UTS
          ++ UP2UTS(p) converts \spad{p} to a Taylor series.
      UTS2UP:   (UTS, NonNegativeInteger) -> UP
          ++ UTS2UP(s, n) converts the first \spad{n} terms of \spad{s}
          ++ to a univariate polynomial.
      LODO2FUN: L -> (List UTS -> UTS)
          ++ LODO2FUN(op) returns the function to pass to the series ODE
          ++ solver in order to solve \spad{op y = 0}.
      if F has IntegralDomain then
          RF2UTS: Fraction UP -> UTS
              ++ RF2UTS(f) converts \spad{f} to a Taylor series.

  Implementation ==> add
      fun: (Vector UTS, List UTS) -> UTS

      UP2UTS p ==
        q := p(monomial(1, 1) + center(0)::UP)
        +/[monomial(coefficient(q, i), i)$UTS for i in 0..degree q]

      UTS2UP(s, n) ==
        xmc     := monomial(1, 1)$UP - center(0)::UP
        xmcn:UP := 1
        ans:UP  := 0
        for i in 0..n repeat
            ans  := ans + coefficient(s, i) * xmcn
            xmcn := xmc * xmcn
        ans

      LODO2FUN op ==
          a := recip(UP2UTS(- leadingCoefficient op))::UTS
          n := (degree(op) - 1)::NonNegativeInteger
          v := [a * UP2UTS coefficient(op, i) for i in 0..n]$Vector(UTS)
          fun(v, #1)

      fun(v, l) ==
          ans:UTS := 0
          for b in l for i in 1.. repeat ans := ans + v.i * b
          ans

      if F has IntegralDomain then
          RF2UTS f == UP2UTS(numer f) * recip(UP2UTS denom f)::UTS

@
\section{package ODERAT RationalLODE}
<<package ODERAT RationalLODE>>=
)abbrev package ODERAT RationalLODE
++ Author: Manuel Bronstein
++ Date Created: 13 March 1991
++ Date Last Updated: 13 April 1994
++ Description:
++  \spad{RationalLODE} provides functions for in-field solutions of linear
++   ordinary differential equations, in the rational case.
RationalLODE(F, UP): Exports == Implementation where
  F  : Join(Field, CharacteristicZero, RetractableTo Integer,
                                       RetractableTo Fraction Integer)
  UP : UnivariatePolynomialCategory F

  N   ==> NonNegativeInteger
  Z   ==> Integer
  RF  ==> Fraction UP
  U   ==> Union(RF, "failed")
  V   ==> Vector F
  M   ==> Matrix F
  LODO ==> LinearOrdinaryDifferentialOperator1 RF
  LODO2==> LinearOrdinaryDifferentialOperator2(UP, RF)

  Exports ==> with
    ratDsolve: (LODO, RF) -> Record(particular: U, basis: List RF)
      ++ ratDsolve(op, g) returns \spad{["failed", []]} if the equation
      ++ \spad{op y = g} has no rational solution. Otherwise, it returns
      ++ \spad{[f, [y1,...,ym]]} where f is a particular rational solution
      ++ and the yi's form a basis for the rational solutions of the
      ++ homogeneous equation.
    ratDsolve: (LODO, List RF) -> Record(basis:List RF, mat:Matrix F)
      ++ ratDsolve(op, [g1,...,gm]) returns \spad{[[h1,...,hq], M]} such
      ++ that any rational solution of \spad{op y = c1 g1 + ... + cm gm}
      ++ is of the form \spad{d1 h1 + ... + dq hq} where
      ++ \spad{M [d1,...,dq,c1,...,cm] = 0}.
    ratDsolve: (LODO2, RF) -> Record(particular: U, basis: List RF)
      ++ ratDsolve(op, g) returns \spad{["failed", []]} if the equation
      ++ \spad{op y = g} has no rational solution. Otherwise, it returns
      ++ \spad{[f, [y1,...,ym]]} where f is a particular rational solution
      ++ and the yi's form a basis for the rational solutions of the
      ++ homogeneous equation.
    ratDsolve: (LODO2, List RF) -> Record(basis:List RF, mat:Matrix F)
      ++ ratDsolve(op, [g1,...,gm]) returns \spad{[[h1,...,hq], M]} such
      ++ that any rational solution of \spad{op y = c1 g1 + ... + cm gm}
      ++ is of the form \spad{d1 h1 + ... + dq hq} where
      ++ \spad{M [d1,...,dq,c1,...,cm] = 0}.
    indicialEquationAtInfinity: LODO -> UP
      ++ indicialEquationAtInfinity op returns the indicial equation of
      ++ \spad{op} at infinity.
    indicialEquationAtInfinity: LODO2 -> UP
      ++ indicialEquationAtInfinity op returns the indicial equation of
      ++ \spad{op} at infinity.

  Implementation ==> add
    import BoundIntegerRoots(F, UP)
    import RationalIntegration(F, UP)
    import PrimitiveRatDE(F, UP, LODO2, LODO)
    import LinearSystemMatrixPackage(F, V, V, M)
    import InnerCommonDenominator(UP, RF, List UP, List RF)

    nzero?             : V -> Boolean
    evenodd            : N -> F
    UPfact             : N -> UP
    infOrder           : RF -> Z
    infTau             : (UP, N) -> F
    infBound           : (LODO2, List RF) -> N
    regularPoint       : (LODO2, List RF) -> Z
    infIndicialEquation: (List N, List UP) -> UP
    makeDot            : (Vector F, List RF) -> RF
    unitlist           : (N, N) -> List F
    infMuLambda: LODO2 -> Record(mu:Z, lambda:List N, func:List UP)
    ratDsolve0: (LODO2, RF) -> Record(particular: U, basis: List RF)
    ratDsolve1: (LODO2, List RF) -> Record(basis:List RF, mat:Matrix F)
    candidates: (LODO2,List RF,UP) -> Record(basis:List RF,particular:List RF)

    dummy := new()$Symbol

    infOrder f == (degree denom f) - (degree numer f)
    evenodd n  == (even? n => 1; -1)

    ratDsolve1(op, lg) ==
      d := denomLODE(op, lg)
      rec := candidates(op, lg, d)
      l := concat([op q for q in rec.basis],
                  [op(rec.particular.i) - lg.i for i in 1..#(rec.particular)])
      sys1 := reducedSystem(matrix [l])@Matrix(UP)
      [rec.basis, reducedSystem sys1]

    ratDsolve0(op, g) ==
      zero? degree op => [inv(leadingCoefficient(op)::RF) * g, empty()]
      minimumDegree op > 0 =>
        sol := ratDsolve0(monicRightDivide(op, monomial(1, 1)).quotient, g)
        b:List(RF) := [1]
        for f in sol.basis repeat
          if (uu := infieldint f) case RF then b := concat(uu::RF, b)
        sol.particular case "failed" => ["failed", b]
        [infieldint(sol.particular::RF), b]
      (u := denomLODE(op, g)) case "failed" => ["failed", empty()]
      rec := candidates(op, [g], u::UP)
      l := lb := lsol := empty()$List(RF)
      for q in rec.basis repeat
          if zero?(opq := op q) then lsol := concat(q, lsol)
          else (l := concat(opq, l); lb := concat(q, lb))
      h:RF := (zero? g => 0; first(rec.particular))
      empty? l =>
          zero? g => [0, lsol]
          [(g = op h => h; "failed"), lsol]
      m:M
      v:V
      if zero? g then
          m := reducedSystem(reducedSystem(matrix [l])@Matrix(UP))@M
          v := new(ncols m, 0)$V
      else
          sys1 := reducedSystem(matrix [l], vector [g - op h]
                               )@Record(mat: Matrix UP, vec: Vector UP)
          sys2 := reducedSystem(sys1.mat, sys1.vec)@Record(mat:M, vec:V)
          m := sys2.mat
          v := sys2.vec
      sol := solve(m, v)
      part:U :=
        zero? g => 0
        sol.particular case "failed" => "failed"
        makeDot(sol.particular::V, lb) + first(rec.particular)
      [part,
       concat_!(lsol, [makeDot(v, lb) for v in sol.basis | nzero? v])]

    indicialEquationAtInfinity(op:LODO2) ==
      rec := infMuLambda op
      infIndicialEquation(rec.lambda, rec.func)

    indicialEquationAtInfinity(op:LODO) ==
      rec := splitDenominator(op, empty())
      indicialEquationAtInfinity(rec.eq)

    regularPoint(l, lg) ==
      a := leadingCoefficient(l) * commonDenominator lg
      coefficient(a, 0) ~= 0 => 0
      for i in 1.. repeat
        a(j := i::F) ~= 0 => return i
        a(-j) ~= 0 => return(-i)

    unitlist(i, q) ==
      v := new(q, 0)$Vector(F)
      v.i := 1
      parts v

    candidates(op, lg, d) ==
      n := degree d + infBound(op, lg)
      m := regularPoint(op, lg)
      uts := UnivariateTaylorSeries(F, dummy, m::F)
      tools := UTSodetools(F, UP, LODO2, uts)
      solver := UnivariateTaylorSeriesODESolver(F, uts)
      dd := UP2UTS(d)$tools
      f := LODO2FUN(op)$tools
      q := degree op
      e := unitlist(1, q)
      hom := [UTS2UP(dd * ode(f, unitlist(i, q))$solver, n)$tools /$RF d
                   for i in 1..q]$List(RF)
      a1 := inv(leadingCoefficient(op)::RF)
      part := [UTS2UP(dd * ode(RF2UTS(a1 * g)$tools + f #1, e)$solver, n)$tools
                /$RF d for g in lg | g ~= 0]$List(RF)
      [hom, part]

    nzero? v ==
      for i in minIndex v .. maxIndex v repeat
        not zero? qelt(v, i) => return true
      false

-- returns z(z+1)...(z+(n-1))
    UPfact n ==
      zero? n => 1
      z := monomial(1, 1)$UP
      */[z + i::F::UP for i in 0..(n-1)::N]

    infMuLambda l ==
      lamb:List(N) := [d := degree l]
      lf:List(UP) := [a := leadingCoefficient l]
      mup := degree(a)::Z - d
      while (l := reductum l) ~= 0 repeat
          a := leadingCoefficient l
          if (m := degree(a)::Z - (d := degree l)) > mup then
            mup := m
            lamb := [d]
            lf := [a]
          else if (m = mup) then
            lamb := concat(d, lamb)
            lf := concat(a, lf)
      [mup, lamb, lf]

    infIndicialEquation(lambda, lf) ==
      ans:UP := 0
      for i in lambda for f in lf repeat
        ans := ans + evenodd i * leadingCoefficient f * UPfact i
      ans

    infBound(l, lg) ==
      rec := infMuLambda l
      n := min(- degree(l)::Z - 1,
               integerBound infIndicialEquation(rec.lambda, rec.func))
      while not(empty? lg) and zero? first lg repeat lg := rest lg
      empty? lg => (-n)::N
      m := infOrder first lg
      for g in rest lg repeat
        if not(zero? g) and (mm := infOrder g) < m then m := mm
      (-min(n, rec.mu - degree(leadingCoefficient l)::Z + m))::N

    makeDot(v, bas) ==
      ans:RF := 0
      for i in 1.. for b in bas repeat ans := ans + v.i::UP * b
      ans

    ratDsolve(op:LODO, g:RF) ==
      rec := splitDenominator(op, [g])
      ratDsolve0(rec.eq, first(rec.rh))

    ratDsolve(op:LODO, lg:List RF) ==
      rec := splitDenominator(op, lg)
      ratDsolve1(rec.eq, rec.rh)

    ratDsolve(op:LODO2, g:RF) ==
      unit?(c := content op) => ratDsolve0(op, g)
      ratDsolve0((op exquo c)::LODO2, inv(c::RF) * g)

    ratDsolve(op:LODO2, lg:List RF) ==
      unit?(c := content op) => ratDsolve1(op, lg)
      ratDsolve1((op exquo c)::LODO2, [inv(c::RF) * g for g in lg])

@
\section{package ODETOOLS ODETools}
<<package ODETOOLS ODETools>>=
)abbrev package ODETOOLS ODETools
++ Author: Manuel Bronstein
++ Date Created: 20 March 1991
++ Date Last Updated: 2 February 1994
++ Description:
++   \spad{ODETools} provides tools for the linear ODE solver.
ODETools(F, LODO): Exports == Implementation where
  N ==> NonNegativeInteger
  L ==> List F
  V ==> Vector F
  M ==> Matrix F

  F:    Field
  LODO: LinearOrdinaryDifferentialOperatorCategory F

  Exports ==> with
    wronskianMatrix: L -> M
      ++ wronskianMatrix([f1,...,fn]) returns the \spad{n x n} matrix m
      ++ whose i^th row is \spad{[f1^(i-1),...,fn^(i-1)]}.
    wronskianMatrix: (L, N) -> M
      ++ wronskianMatrix([f1,...,fn], q, D) returns the \spad{q x n} matrix m
      ++ whose i^th row is \spad{[f1^(i-1),...,fn^(i-1)]}.
    variationOfParameters: (LODO, F, L) -> Union(V, "failed")
      ++ variationOfParameters(op, g, [f1,...,fm])
      ++ returns \spad{[u1,...,um]} such that a particular solution of the
      ++ equation \spad{op y = g} is \spad{f1 int(u1) + ... + fm int(um)}
      ++ where \spad{[f1,...,fm]} are linearly independent and \spad{op(fi)=0}.
      ++ The value "failed" is returned if \spad{m < n} and no particular
      ++ solution is found.
    particularSolution: (LODO, F, L, F -> F) -> Union(F, "failed")
      ++ particularSolution(op, g, [f1,...,fm], I) returns a particular
      ++ solution h of the equation \spad{op y = g} where \spad{[f1,...,fm]}
      ++ are linearly independent and \spad{op(fi)=0}.
      ++ The value "failed" is returned if no particular solution is found.
      ++ Note: the method of variations of parameters is used.

  Implementation ==> add
    import LinearSystemMatrixPackage(F, V, V, M)

    diff := D()$LODO

    wronskianMatrix l == wronskianMatrix(l, #l)

    wronskianMatrix(l, q) ==
      v:V := vector l
      m:M := zero(q, #v)
      for i in minRowIndex m .. maxRowIndex m repeat
        setRow_!(m, i, v)
        v := map_!(diff #1, v)
      m

    variationOfParameters(op, g, b) ==
      empty? b => "failed"
      v:V := new(n := degree op, 0)
      qsetelt_!(v, maxIndex v, g / leadingCoefficient op)
      particularSolution(wronskianMatrix(b, n), v)

    particularSolution(op, g, b, integration) ==
      zero? g => 0
      (sol := variationOfParameters(op, g, b)) case "failed" => "failed"
      ans:F := 0
      for f in b for i in minIndex(s := sol::V) .. repeat
        ans := ans + integration(qelt(s, i)) * f
      ans

@
\section{package ODEINT ODEIntegration}
<<package ODEINT ODEIntegration>>=
)abbrev package ODEINT ODEIntegration
++ Author: Manuel Bronstein
++ Date Created: 4 November 1991
++ Date Last Updated: 2 February 1994
++ Description:
++ \spadtype{ODEIntegration} provides an interface to the integrator.
++ This package is intended for use
++ by the differential equations solver but not at top-level.
ODEIntegration(R, F): Exports == Implementation where
  R: Join(EuclideanDomain, RetractableTo Integer,
          LinearlyExplicitRingOver Integer, CharacteristicZero)
  F: Join(AlgebraicallyClosedFunctionSpace R, TranscendentalFunctionCategory,
                                              PrimitiveFunctionCategory)

  Q   ==> Fraction Integer
  UQ  ==> Union(Q, "failed")
  SY  ==> Symbol
  K   ==> Kernel F
  P   ==> SparseMultivariatePolynomial(R, K)
  REC ==> Record(coef:Q, logand:F)

  Exports ==> with
    int   : (F, SY) -> F
      ++ int(f, x) returns the integral of f with respect to x.
    expint: (F, SY) -> F
      ++ expint(f, x) returns e^{the integral of f with respect to x}.
    diff  : SY -> (F -> F)
      ++ diff(x) returns the derivation with respect to x.

  Implementation ==> add
    import FunctionSpaceIntegration(R, F)
    import ElementaryFunctionStructurePackage(R, F)

    isQ   : List F -> UQ
    isQlog: F -> Union(REC, "failed")
    mkprod: List REC -> F

    diff x == differentiate(#1, x)

-- This is the integration function to be used for quadratures
    int(f, x) ==
      (u := integrate(f, x)) case F => u::F
      first(u::List(F))

-- mkprod([q1, f1],...,[qn,fn]) returns */(fi^qi) but groups the
-- qi having the same denominator together
    mkprod l ==
      empty? l => 1
      rec := first l
      d := denom(rec.coef)
      ll := select(denom(#1.coef) = d, l)
      nthRoot(*/[r.logand ** numer(r.coef) for r in ll], d) *
        mkprod setDifference(l, ll)

-- computes exp(int(f,x)) in a non-naive way
    expint(f, x) ==
      a := int(f, x)
      (u := validExponential(tower a, a, x)) case F => u::F
      da := denom a
      l :=
        (v := isPlus(na := numer a)) case List(P) => v::List(P)
        [na]
      exponent:P := 0
      lrec:List(REC) := empty()
      for term in l repeat
        if (w := isQlog(term / da)) case REC then
          lrec := concat(w::REC, lrec)
        else
          exponent := exponent + term
      mkprod(lrec) * exp(exponent / da)

-- checks if all the elements of l are rational numbers, returns their product
    isQ l ==
      prod:Q := 1
      for x in l repeat
        (u := retractIfCan(x)@UQ) case "failed" => return "failed"
        prod := prod * u::Q
      prod

-- checks if a non-sum expr is of the form c * log(g) for a rational number c
    isQlog f ==
      is?(f, 'log) => [1, first argument(retract(f)@K)]
      (v := isTimes f) case List(F) and (#(l := v::List(F)) <= 3) =>
          l := reverse! sort!(before?,l)
          is?(first l, 'log) and ((u := isQ rest l) case Q) =>
              [u::Q, first argument(retract(first(l))@K)]
          "failed"
      "failed"

@
\section{package ODECONST ConstantLODE}
<<package ODECONST ConstantLODE>>=
)abbrev package ODECONST ConstantLODE
++ Author: Manuel Bronstein
++ Date Created: 18 March 1991
++ Date Last Updated: 3 February 1994
++ Description: Solution of linear ordinary differential equations, constant coefficient case.
ConstantLODE(R, F, L): Exports == Implementation where
  R: Join(EuclideanDomain, RetractableTo Integer,
          LinearlyExplicitRingOver Integer, CharacteristicZero)
  F: Join(AlgebraicallyClosedFunctionSpace R,
          TranscendentalFunctionCategory, PrimitiveFunctionCategory)
  L: LinearOrdinaryDifferentialOperatorCategory F

  Z   ==> Integer
  SY  ==> Symbol
  K   ==> Kernel F
  V   ==> Vector F
  M   ==> Matrix F
  SUP ==> SparseUnivariatePolynomial F

  Exports ==> with
    constDsolve: (L, F, SY) -> Record(particular:F, basis:List F)
      ++ constDsolve(op, g, x) returns \spad{[f, [y1,...,ym]]}
      ++ where f is a particular solution of the equation \spad{op y = g},
      ++ and the \spad{yi}'s form a basis for the solutions of \spad{op y = 0}.

  Implementation ==> add
    import ODETools(F, L)
    import ODEIntegration(R, F)
    import ElementaryFunctionSign(R, F)
    import AlgebraicManipulations(R, F)
    import FunctionSpaceIntegration(R, F)
    import FunctionSpaceUnivariatePolynomialFactor(R, F, SUP)

    homoBasis: (L, F) -> List F
    quadSol  : (SUP, F) -> List F
    basisSqfr: (SUP, F) -> List F
    basisSol : (SUP, Z, F) -> List F

    constDsolve(op, g, x) ==
      b := homoBasis(op, x::F)
      [particularSolution(op, g, b, int(#1, x))::F, b]

    homoBasis(op, x) ==
      p:SUP := 0
      while op ~= 0 repeat
          p  := p + monomial(leadingCoefficient op, degree op)
          op := reductum op
      b:List(F) := empty()
      for ff in factors ffactor p repeat
        b := concat_!(b, basisSol(ff.factor, dec(ff.exponent), x))
      b

    basisSol(p, n, x) ==
      l := basisSqfr(p, x)
      zero? n => l
      ll := copy l
      xn := x::F
      for i in 1..n repeat
        l := concat_!(l, [xn * f for f in ll])
        xn := x * xn
      l

    basisSqfr(p, x) ==
      one?(d := degree p) =>
        [exp(- coefficient(p, 0) * x / leadingCoefficient p)]
      d = 2 => quadSol(p, x)
      [exp(a * x) for a in rootsOf p]

    quadSol(p, x) ==
      (u := sign(delta := (b := coefficient(p, 1))**2 - 4 *
        (a := leadingCoefficient p) * (c := coefficient(p, 0))))
          case Z and negative?(u::Z) =>
            y := x / (2 * a)
            r := - b * y
            i := rootSimp(sqrt(-delta)) * y
            [exp(r) * cos(i), exp(r) * sin(i)]
      [exp(a * x) for a in zerosOf p]

@
\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>>

-- Compile order for the differential equation solver:
-- oderf.spad  odealg.spad  nlode.spad  nlinsol.spad  riccati.spad  odeef.spad

<<package BALFACT BalancedFactorisation>>
<<package BOUNDZRO BoundIntegerRoots>>
<<package ODEPRIM PrimitiveRatDE>>
<<package UTSODETL UTSodetools>>
<<package ODERAT RationalLODE>>
<<package ODETOOLS ODETools>>
<<package ODEINT ODEIntegration>>
<<package ODECONST ConstantLODE>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}