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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra lodof.spad}
\author{Manuel Bronstein, Fritz Schwarz}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain SETMN SetOfMIntegersInOneToN}
<<domain SETMN SetOfMIntegersInOneToN>>=
)abbrev domain SETMN SetOfMIntegersInOneToN
++ Author: Manuel Bronstein
++ Date Created: 10 January 1994
++ Date Last Updated: 10 January 1994
++ Description:
++ \spadtype{SetOfMIntegersInOneToN} implements the subsets of M integers
++ in the interval \spad{[1..n]}
SetOfMIntegersInOneToN(m, n): Exports == Implementation where
  PI ==> PositiveInteger
  N  ==> NonNegativeInteger
  U  ==> Union(%, "failed")
  n,m: PI
 
  Exports ==> Finite with
    incrementKthElement: (%, PI) -> U
      ++ incrementKthElement(S,k) increments the k^{th} element of S,
      ++ and returns "failed" if the result is not a set of M integers
      ++ in \spad{1..n} any more.
    replaceKthElement:   (%, PI, PI) -> U
      ++ replaceKthElement(S,k,p) replaces the k^{th} element of S by p,
      ++ and returns "failed" if the result is not a set of M integers
      ++ in \spad{1..n} any more.
    elements: % -> List PI
      ++ elements(S) returns the list of the elements of S in increasing order.
    setOfMinN: List PI -> %
      ++ setOfMinN([a_1,...,a_m]) returns the set {a_1,...,a_m}.
      ++ Error if {a_1,...,a_m} is not a set of M integers in \spad{1..n}.
    enumerate: () -> Vector %
      ++ enumerate() returns a vector of all the sets of M integers in
      ++ \spad{1..n}.
    member?:   (PI, %) -> Boolean
      ++ member?(p, s) returns true is p is in s, false otherwise.
    delta: (%, PI, PI) -> N
      ++ delta(S,k,p) returns the number of elements of S which are strictly
      ++ between p and the k^{th} element of S.
 
  Implementation ==> add
    Rep := Record(bits:Bits, pos:N)
 
    reallyEnumerate: () -> Vector %
    enum: (N, N, PI) -> List Bits
 
    all:Reference Vector % := ref empty()
    sz:Reference N := ref 0
 
    s1 = s2                == s1.bits =$Bits s2.bits
    coerce(s:%):OutputForm == brace [i::OutputForm for i in elements s]
    random()               == index((1 + (random()$Integer rem size()))::PI)
    reallyEnumerate()      == [[b, i] for b in enum(m, n, n) for i in 1..]
    member?(p, s)          == s.bits.p
 
    enumerate() ==
      if empty? deref all then setref(all,reallyEnumerate())
      deref all
 
-- enumerates the sets of p integers in 1..q, returns them as sets in 1..n
-- must have p <= q
    enum(p, q, n) ==
      zero? p or zero? q => empty()
      p = q =>
        b := new(n, false)$Bits
        for i in 1..p repeat b.i := true
        [b]
      q1 := (q - 1)::N
      l := enum((p - 1)::N, q1, n)
      if empty? l then l := [new(n, false)$Bits]
      for s in l repeat s.q := true
      concat!(enum(p, q1, n), l)
 
    size() ==
      if zero? deref sz then
        setref(sz,binomial(n, m)$IntegerCombinatoricFunctions(Integer) :: N)
      deref sz
 
    lookup s ==
      if empty? deref all then setref(all,reallyEnumerate())
      if zero?(s.pos) then s.pos := position(s, deref all) :: N
      s.pos :: PI
      
    index p ==
      p > size() => error "index: argument too large"
      if empty? deref all then setref(all,reallyEnumerate())
      deref(all).p
 
    setOfMinN l ==
      s := new(n, false)$Bits
      count:N := 0
      for i in l repeat
        count := count + 1
        count > m or zero? i or i > n or s.i =>
          error "setOfMinN: improper set of integers"
        s.i := true
      count < m => error "setOfMinN: improper set of integers"
      [s, 0]
 
    elements s ==
      b := s.bits
      l:List PI := empty()
      found:N := 0
      i:PI := 1
      while found < m repeat
          if b.i then
              l := concat(i, l)
              found := found + 1
          i := i + 1
      reverse! l
 
    incrementKthElement(s, k) ==
      b := s.bits
      found:N := 0
      i:N := 1
      while found < k repeat
          if b.i then found := found + 1
          i := i + 1
      i > n or b.i => "failed"
      newb := copy b
      newb.i := true
      newb.((i-1)::N) := false
      [newb, 0]
 
    delta(s, k, p) ==
      b := s.bits
      count:N := found:N := 0
      i:PI := 1
      while found < k repeat
          if b.i then
             found := found + 1
             if i > p and found < k then count := count + 1
          i := i + 1
      count
 
    replaceKthElement(s, k, p) ==
      b := s.bits
      found:N := 0
      i:PI := 1
      while found < k repeat
          if b.i then found := found + 1
          if found < k then i := i + 1
      b.p and i ~= p => "failed"
      newb := copy b
      newb.p := true
      newb.i := false
      [newb, (i = p => s.pos; 0)]

@
\section{package PREASSOC PrecomputedAssociatedEquations}
<<package PREASSOC PrecomputedAssociatedEquations>>=
)abbrev package PREASSOC PrecomputedAssociatedEquations
++ Author: Manuel Bronstein
++ Date Created: 13 January 1994
++ Date Last Updated: 3 February 1994
++ Description:
++ \spadtype{PrecomputedAssociatedEquations} stores some generic
++ precomputations which speed up the computations of the
++ associated equations needed for factoring operators.
PrecomputedAssociatedEquations(R, L): Exports == Implementation where
  R: IntegralDomain
  L: LinearOrdinaryDifferentialOperatorCategory R
 
  PI  ==> PositiveInteger
  N   ==> NonNegativeInteger
  A   ==> PrimitiveArray R
  U   ==> Union(Matrix R, "failed")
 
  Exports ==> with
    firstUncouplingMatrix: (L, PI) -> U
      ++ firstUncouplingMatrix(op, m) returns the matrix A such that
      ++ \spad{A w = (W',W'',...,W^N)} in the corresponding associated
      ++ equations for right-factors of order m of op.
      ++ Returns "failed" if the matrix A has not been precomputed for
      ++ the particular combination \spad{degree(L), m}.
 
  Implementation ==> add
    A32:  L -> U
    A42:  L -> U
    A425: (A, A, A) -> List R
    A426: (A, A, A) -> List R
    makeMonic: L -> Union(A, "failed")
 
    diff:L := D()
 
    firstUncouplingMatrix(op, m) ==
      n := degree op
      n = 3 and m = 2 => A32 op
      n = 4 and m = 2 => A42 op
      "failed"
        
    makeMonic op ==
      lc := leadingCoefficient op
      a:A := new(n := degree op, 0)
      for i in 0..(n-1)::N repeat
        (u := coefficient(op, i) exquo lc) case "failed" => return "failed"
        a.i := - (u::R)
      a
    
    A32 op ==
      (u := makeMonic op) case "failed" => "failed"
      a := u::A
      matrix [[0, 1, 0], [a.1, a.2, 1],
              [diff(a.1) + a.1 * a.2 - a.0, diff(a.2) + a.2**2 + a.1, 2 * a.2]]
 
    A42 op ==
      (u := makeMonic op) case "failed" => "failed"
      a := u::A
      a':A := new(4, 0)
      a'':A := new(4, 0)
      for i in 0..3 repeat
        a'.i := diff(a.i)
        a''.i := diff(a'.i)
      matrix [[0, 1, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0], [a.1,a.2,0,a.3,2::R,0],
              [a'.1 + a.1 * a.3 - 2 * a.0, a'.2 + a.2 * a.3 + a.1, 3 * a.2,
               a'.3 + a.3 ** 2 + a.2, 3 * a.3, 2::R],
                A425(a, a', a''), A426(a, a', a'')]
 
    A425(a, a', a'') ==
      [a''.1 + 2 * a.1 * a'.3 + a.3 * a'.1 - 2 * a'.0 + a.1 * a.3 ** 2
       - 3 * a.0 * a.3 + a.1 * a.2,
        a''.2 + 2 * a.2 * a'.3 + a.3 * a'.2 + 2 * a'.1 + a.2 * a.3 ** 2
         + a.1 * a.3 + a.2 ** 2 - 4 * a.0,
          4 * a'.2 + 4 * a.2 * a.3 - a.1,
           a''.3 + 3 * a.3 * a'.3 + 2 * a'.2 + a.3 ** 3 + 2 * a.2 * a.3 + a.1,
            4 * a'.3 + 4 * a.3 ** 2 + 4 * a.2, 5 * a.3]
              
    A426(a, a', a'') ==
      [diff(a''.1) + 3 * a.1 * a''.3 + a.3 * a''.1 - 2 * a''.0
       + (3 * a'.1 + 5 * a.1 * a.3 - 7 * a.0) * a'.3 + 3 * a.1 * a'.2
        + (a.3 ** 2 + a.2) * a'.1 - 3 * a.3 * a'.0 + a.1 * a.3 ** 3
         - 4 * a.0 * a.3 ** 2 + 2 * a.1 * a.2 * a.3 - 4 * a.0 * a.2 + a.1 ** 2,
          diff(a''.2) + 3 * a.2 * a''.3 + a.3 * a''.2 + 3 * a''.1
           + (3*a'.2 + 5*a.2 * a.3 + 3 * a.1) * a'.3 + (a.3**2 + 4*a.2)*a'.2
            + 2 * a.3 * a'.1 - 6 * a'.0 + a.2 * a.3 ** 3 + a.1 * a.3 ** 2
             + (2 * a.2**2 - 8 * a.0) * a.3 + 2 * a.1 * a.2,
              5 * a''.2 + 10 * a.2 * a'.3 + 5 * a.3 * a'.2 + a'.1
               + 5 * a.2 * a.3 ** 2 - 4 * a.1 * a.3 + 5 * a.2**2 - 4 * a.0,
                diff(a''.3) + 4 * a.3 * a''.3 + 3*a''.2 + 3 * a'.3**2
                 + (6 * a.3**2 + 4 * a.2) * a'.3 + 5 * a.3 * a'.2 + 3 * a'.1
                  + a.3**4 + 3 * a.2 * a.3**2 + 2 * a.1 * a.3 + a.2**2 - 4*a.0,
                   5 * a''.3 + 15 * a.3 * a'.3 + 10 * a'.2 + 5 * a.3**3
                    + 10 * a.2 * a.3, 9 * a'.3 + 9 * a.3**2 + 4 * a.2]

@
\section{package ASSOCEQ AssociatedEquations}
<<package ASSOCEQ AssociatedEquations>>=
)abbrev package ASSOCEQ AssociatedEquations
++ Author: Manuel Bronstein
++ Date Created: 10 January 1994
++ Date Last Updated: 3 February 1994
++ Description:
++ \spadtype{AssociatedEquations} provides functions to compute the
++ associated equations needed for factoring operators
AssociatedEquations(R, L):Exports == Implementation where
  R: IntegralDomain
  L: LinearOrdinaryDifferentialOperatorCategory R
 
  PI  ==> PositiveInteger
  N   ==> NonNegativeInteger
  MAT ==> Matrix R
  REC ==> Record(minor: List PI, eq: L, minors: List List PI, ops: List L)
 
  Exports ==> with
    associatedSystem: (L, PI) -> Record(mat: MAT, vec:Vector List PI)
      ++ associatedSystem(op, m) returns \spad{[M,w]} such that the
      ++ m-th associated equation system to L is \spad{w' = M w}.
    uncouplingMatrices: MAT -> Vector MAT
      ++ uncouplingMatrices(M) returns \spad{[A_1,...,A_n]} such that if
      ++ \spad{y = [y_1,...,y_n]} is a solution of \spad{y' = M y}, then
      ++ \spad{[$y_j',y_j'',...,y_j^{(n)}$] = $A_j y$} for all j's.
    if R has Field then
        associatedEquations: (L, PI) -> REC
          ++ associatedEquations(op, m) returns \spad{[w, eq, lw, lop]}
          ++ such that \spad{eq(w) = 0} where w is the given minor, and
          ++ \spad{lw_i = lop_i(w)} for all the other minors.
 
  Implementation ==> add
    makeMatrix: (Vector MAT, N) -> MAT
 
    diff:L := D()
 
    makeMatrix(v, n) == matrix [parts row(v.i, n) for i in 1..#v]
 
    associatedSystem(op, m) ==
      eq: Vector R
      S := SetOfMIntegersInOneToN(m, n := degree(op)::PI)
      w := enumerate()$S
      s := size()$S
      ww:Vector List PI := new(s, empty())
      M:MAT := new(s, s, 0)
      m1 := (m::Integer - 1)::PI
      an := leadingCoefficient op
      a:Vector(R) := [- (coefficient(op, j) exquo an)::R for j in 0..n - 1]
      for i in 1..s repeat
          eq := new(s, 0)
          wi := w.i
          ww.i := elements wi
          for k in 1..m1 repeat
              u := incrementKthElement(wi, k::PI)$S
              if u case S then eq(lookup(u::S)) := 1
          if member?(n, wi) then
              for j in 1..n | a.j ~= 0 repeat
                  u := replaceKthElement(wi, m, j::PI)
                  if u case S then
                    eq(lookup(u::S)) := (odd? delta(wi, m, j::PI) => -a.j; a.j)
          else
              u := incrementKthElement(wi, m)$S
              if u case S then eq(lookup(u::S)) := 1
          setRow!(M, i, eq)
      [M, ww]
 
    uncouplingMatrices m ==
      n := nrows m
      v:Vector MAT := new(n, zero(1, 0)$MAT)
      v.1 := mi := m
      for i in 2..n repeat v.i := mi := map(diff #1, mi) + mi * m
      [makeMatrix(v, i) for i in 1..n]
 
    if R has Field then
        import PrecomputedAssociatedEquations(R, L)
 
        makeop:    Vector R -> L
        makeeq:    (Vector List PI, MAT, N, N) -> REC
        computeIt: (L, PI, N) -> REC
 
        makeeq(v, m, i, n) ==
          [v.i, makeop row(m, i) - 1, [v.j for j in 1..n | j ~= i],
                                    [makeop row(m, j) for j in 1..n | j ~= i]]
 
        associatedEquations(op, m) ==
          (u := firstUncouplingMatrix(op, m)) case "failed" => computeIt(op,m,1)
          (v := inverse(u::MAT)) case "failed" => computeIt(op, m, 2)
          S := SetOfMIntegersInOneToN(m, degree(op)::PI)
          w := enumerate()$S
          s := size()$S
          ww:Vector List PI := new(s, empty())
          for i in 1..s repeat ww.i := elements(w.i)
          makeeq(ww, v::MAT, 1, s)
 
        computeIt(op, m, k) ==
          rec := associatedSystem(op, m)
          a := uncouplingMatrices(rec.mat)
          n := #a
          for i in k..n repeat
            (u := inverse(a.i)) case MAT => return makeeq(rec.vec,u::MAT,i,n)
          error "associatedEquations: full degenerate case"
 
        makeop v ==
          op:L := 0
          for i in 1..#v repeat op := op + monomial(v i, i)
          op

@
\section{package LODOF LinearOrdinaryDifferentialOperatorFactorizer}
<<package LODOF LinearOrdinaryDifferentialOperatorFactorizer>>=
)abbrev package LODOF LinearOrdinaryDifferentialOperatorFactorizer
++ Author: Fritz Schwarz, Manuel Bronstein
++ Date Created: 1988
++ Date Last Updated: 3 February 1994
++ Description:
++ \spadtype{LinearOrdinaryDifferentialOperatorFactorizer} provides a
++ factorizer for linear ordinary differential operators whose coefficients
++ are rational functions.
++ Keywords: differential equation, ODE, LODO, factoring
LinearOrdinaryDifferentialOperatorFactorizer(F, UP): Exports == Impl where
  F : Join(Field, CharacteristicZero,
           RetractableTo Integer, RetractableTo Fraction Integer)
  UP: UnivariatePolynomialCategory F
 
  RF ==> Fraction UP
  L  ==> LinearOrdinaryDifferentialOperator1 RF
 
  Exports ==> with
    factor: (L, UP -> List F) -> List L
      ++ factor(a, zeros) returns the factorisation of a.
      ++ \spad{zeros} is a zero finder in \spad{UP}.
    if F has AlgebraicallyClosedField then
      factor: L -> List L
        ++ factor(a) returns the factorisation of a.
      factor1: L -> List L
        ++ factor1(a) returns the factorisation of a,
        ++ assuming that a has no first-order right factor.
 
  Impl ==> add
    import RationalLODE(F, UP)
    import RationalRicDE(F, UP)
--  import AssociatedEquations RF
 
    dd := D()$L
 
    expsol     : (L, UP -> List F, UP -> Factored UP) -> Union(RF, "failed")
    expsols    : (L, UP -> List F, UP -> Factored UP, Boolean) -> List RF
    opeval     : (L, L) -> L
    recurfactor: (L, L, UP -> List F, UP -> Factored UP, Boolean) -> List L
    rfactor    : (L, L, UP -> List F, UP -> Factored UP, Boolean) -> List L
    rightFactor: (L, NonNegativeInteger, UP -> List F, UP -> Factored UP)
                                                          -> Union(L, "failed")
    innerFactor: (L, UP -> List F, UP -> Factored UP, Boolean) -> List L
 
    factor(l, zeros) == innerFactor(l, zeros, squareFree, true)
 
    expsol(l, zeros, ezfactor) ==
      empty?(sol := expsols(l, zeros, ezfactor, false)) => "failed"
      first sol
 
    expsols(l, zeros, ezfactor, all?) ==
      sol := [differentiate(f)/f for f in ratDsolve(l, 0).basis | f ~= 0]
      not(all? or empty? sol) => sol
      concat(sol, ricDsolve(l, zeros, ezfactor))
 
-- opeval(l1, l2) returns l1(l2)
    opeval(l1, l2) ==
      ans:L := 0
      l2n:L := 1
      for i in 0..degree l1 repeat
        ans := ans + coefficient(l1, i) * l2n
        l2n := l2 * l2n
      ans
      
    recurfactor(l, r, zeros, ezfactor, adj?) ==
      q := rightExactQuotient(l, r)::L
      if adj? then q := adjoint q
      innerFactor(q, zeros, ezfactor, true)
 
    rfactor(op, r, zeros, ezfactor, adj?) ==
      degree r > 1 or not one? leadingCoefficient r =>
        recurfactor(op, r, zeros, ezfactor, adj?)
      op1 := opeval(op, dd - coefficient(r, 0)::L)
      map!(opeval(#1, r), recurfactor(op1, dd, zeros, ezfactor, adj?))
 
-- r1? is true means look for 1st-order right-factor also
    innerFactor(l, zeros, ezfactor, r1?) ==
      (n := degree l) <= 1 => [l]
      ll := adjoint l
      for i in 1..(n quo 2) repeat
        (r1? or (i > 1)) and ((u := rightFactor(l,i,zeros,ezfactor)) case L) =>
           return concat!(rfactor(l, u::L, zeros, ezfactor, false), u::L)
        (2 * i < n) and ((u := rightFactor(ll, i, zeros, ezfactor)) case L) =>
           return concat(adjoint(u::L), rfactor(ll, u::L, zeros,ezfactor,true))
      [l]
 
    rightFactor(l, n, zeros, ezfactor) ==
      one? n =>
        (u := expsol(l, zeros, ezfactor)) case "failed" => "failed"
        D() - u::RF::L
--    rec := associatedEquations(l, n::PositiveInteger)
--    empty?(sol := expsols(rec.eq, zeros, ezfactor, true)) => "failed"
      "failed"
 
    if F has AlgebraicallyClosedField then
      zro1: UP -> List F
      zro : (UP, UP -> Factored UP) -> List F
 
      zro(p, ezfactor) ==
        concat [zro1(r.factor) for r in factors ezfactor p]
 
      zro1 p ==
        [zeroOf(map(#1, p)$UnivariatePolynomialCategoryFunctions2(F, UP,
                                             F, SparseUnivariatePolynomial F))]
 
      if F is AlgebraicNumber then
        import AlgFactor UP
 
        factor l  == innerFactor(l, zro(#1, factor), factor, true)
        factor1 l == innerFactor(l, zro(#1, factor), factor, false)
 
      else
        factor l  == innerFactor(l, zro(#1, squareFree), squareFree, true)
        factor1 l == innerFactor(l, zro(#1, squareFree), squareFree, false)

@
\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
-- kovacic.spad  lodof.spad  odeef.spad
 
<<domain SETMN SetOfMIntegersInOneToN>>
<<package PREASSOC PrecomputedAssociatedEquations>>
<<package ASSOCEQ AssociatedEquations>>
<<package LODOF LinearOrdinaryDifferentialOperatorFactorizer>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}