skpol.spad

)abbrev category MORECAT MultivariateSkewPolynomialCategory
MultivariateSkewPolynomialCategory(R : Ring, E : OrderedAbelianMonoidSup,
    Var : OrderedSet) : Category == MaybeSkewPolynomialCategory(R, E, Var)


)abbrev category SSKPOLC SolvableSkewPolynomialCategory
++ Description:
++   This is the category of polynomials in noncommutative
++   variables over noncommutative rings.  We do not
++   assume that variables and elements of the base
++   ring commute.  We assume that the polynomial ring
++   is of solvable type, so noncommutative version
++   of Buchberger algorithm works.

SolvableSkewPolynomialCategory(R : Ring,
  Expon : OrderedAbelianMonoidSup) : Category == Join(Ring, LeftModule(R)) with
    leadingCoefficient : % -> R
      ++ leadingCoefficient(p) returns the coefficient of the highest
      ++ degree term of p.

    leadingMonomial : % -> %
      ++ leadingMonomial(p) returns the monomial of p with the highest degree.

    degree : % -> Expon
      ++ degree(p) returns the maximum of the exponents of the terms of p.

    reductum : % -> %
      ++ reductum(u) returns u minus its leading monomial
      ++ returns zero if handed the zero element.

    monomial : (R, Expon) -> %
      ++ monomial(r, e) makes a term from a coefficient r and an exponent e.



)abbrev package NGROEB NGroebnerPackage

++ Description:
++   This is package computes noncommutative Groebner basis.
++   Based on commutative version.  Note that this package
++   accepts rings as base domain, however computed basis
++   is over left fraction field.  Computations are done
++   in fraction free way (coefficients stay in base ring).

NGroebnerPackage(Dom, Expon, Dpol) : T == C where

  Dom : LeftOreRing
  Expon : OrderedAbelianMonoidSup
  Dpol :  SolvableSkewPolynomialCategory(Dom, Expon)
  NNI    ==> NonNegativeInteger
  critPair ==> Record( lcmfij : Expon, totdeg : NonNegativeInteger,
                      poli : Dpol, polj : Dpol )
  sugarPol ==> Record( totdeg : NonNegativeInteger, pol : Dpol)

  T == with
     groebner : List(Dpol) -> List(Dpol)
       ++ groebner(lp) computes a groebner basis for a polynomial ideal
       ++ generated by the list of polynomials lp.

     redPol : (Dpol, List(Dpol))  -> Dpol
       ++ normalForm(poly, gb) reduces the polynomial poly modulo the
       ++ precomputed groebner basis gb giving up to a constant factor
       ++ a canonical representative of the residue class.

     hMonic : Dpol  -> Dpol
       ++ hMonic(p) tries to remove content from p
     virtualDegree : Dpol -> NonNegativeInteger
       ++ virtualDegree
     sPol : critPair  -> Dpol
       ++ sPol

  C == add

   import from OutputForm

   credPol : (Dpol, List(Dpol))  -> Dpol
   critM : (Expon, Expon) -> Boolean
   redPo : (Dpol, List(Dpol) )  ->  Record(poly : Dpol, mult : Dom)
   -- hMonic:  Dpol  -> Dpol
   updatF : (Dpol, NNI, List(sugarPol) ) -> List(sugarPol)
   -- sPol:  critPair  -> Dpol
   updatD : (List(critPair), List(critPair)) -> List(critPair)
   critpOrder : (critPair, critPair) -> Boolean
   makeCrit : (sugarPol, Dpol, NonNegativeInteger) -> critPair
   -- virtualDegree : Dpol -> NonNegativeInteger
   lcmCo : (Dom, Dom) -> Record(co1 : Dom, co2 : Dom)

   lcmCo(c1, c2) ==
       lcc := lcmCoef(c1, c2)
       [lcc.coeff1, lcc.coeff2]

   ------  Definition of intermediate functions
   if Dpol has totalDegree : Dpol -> NonNegativeInteger then
     virtualDegree p == totalDegree p
   else
     virtualDegree p == 0

   ------  ordering of critpairs

   critpOrder(cp1, cp2) ==
     cp1.totdeg < cp2.totdeg => true
     cp2.totdeg < cp1.totdeg => false
     cp1.lcmfij < cp2.lcmfij

   ------    creating a critical pair

   makeCrit(sp1, p2, totdeg2) ==
     p1 := sp1.pol
     deg := sup(degree(p1), degree(p2))
     e1 := subtractIfCan(deg, degree(p1))::Expon
     e2 := subtractIfCan(deg, degree(p2))::Expon
     tdeg := max(sp1.totdeg + virtualDegree(monomial(1, e1)),
                 totdeg2 + virtualDegree(monomial(1, e2)))
     [deg, tdeg, p1, p2]$critPair

   ------    calculate basis

   gbasis(Pol : List(Dpol)) : List(Dpol) ==
     D, D1 : List(critPair)
     ---------   create D and Pol

     Pol1 := sort((z1, z2) +-> degree z1 > degree z2, Pol)
     basPols := updatF(hMonic(first Pol1), virtualDegree(first Pol1), [])
     Pol1 := rest(Pol1)
     D := []
     while not(empty?(Pol1)) repeat
        h := hMonic(first(Pol1))
        Pol1 := rest(Pol1)
        toth := virtualDegree h
        D1 := [makeCrit(x, h, toth) for x in basPols]
        D := updatD(sort(critpOrder, D1), D)
        basPols := updatF(h, toth, basPols)
     D := sort(critpOrder, D)

     redPols := [x.pol for x in basPols]
     while not(empty?(D)) repeat
         D0 := first D
         s := hMonic(sPol(D0))
         D := rest(D)
         h := hMonic(redPol(s, redPols))
         h = 0  => iterate
         degree(h) = 0 =>
           D := []
           basPols := updatF(h, 0, [])
           break

         D1 := [makeCrit(x, h, D0.totdeg) for x in basPols]
         D := updatD(sort(critpOrder, D1), D)

         basPols := updatF(h, D0.totdeg, basPols)
         redPols := concat(redPols, h)
     Pol := [x.pol for x in basPols]
     Pol

            ----  calculate minimal basis for ordered F

   minGbasis(F : List(Dpol)) : List(Dpol) ==
       empty?(F) => []
       newbas := minGbasis rest F
       cons(hMonic credPol( first(F), newbas), newbas)

   groebner(F) ==
       minGbasis(sort((x, y) +-> degree x > degree y, gbasis(F)))

            --- calculate S-polynomial of a critical pair

   sPol(p : critPair)==
      Tij := p.lcmfij
      fi := p.poli
      fj := p.polj
      fi := monomial(1, subtractIfCan(Tij, degree fi)::Expon)*fi
      fj := monomial(1, subtractIfCan(Tij, degree fj)::Expon)*fj
      cc := lcmCo(leadingCoefficient fi, leadingCoefficient fj)
      cc.co1*fi - cc.co2*fj

   redPo(s : Dpol, F : List(Dpol)) ==
      m : Dom := 1
      Fh := F
      while not(s = 0 or empty?(F)) repeat
        f1 := first(F)
        s1 := degree(s)
        e : Union(Expon, "failed")
        d1 := degree(f1)
        (e := subtractIfCan(s1, d1)) case Expon  =>
           f2 := monomial(1, e)*f1
           cc := lcmCo(leadingCoefficient s, leadingCoefficient f2)
           s := cc.co1*s - cc.co2*f2
           if degree(s) = s1 then
               print(message("no progress in reduction"))
               print(s1::OutputForm)
               print(d1::OutputForm)
               print(e::OutputForm)
               print(degree(f2)::OutputForm)
               print(cc::OutputForm)
               print(s::OutputForm)
               error "no progress in reduction"
           m := m*cc.co1
           F := Fh
        F := rest F
      [s, m]

   redPol(s : Dpol, F : List(Dpol)) ==
      credPol(redPo(s, F).poly, F)

   if Dpol has primitivePart : Dpol -> Dpol then
       hMonic(p) == primitivePart p
   else
       hMonic(p) == p

            ---  reduce all terms of h mod F  (iterative version )

   credPol(h : Dpol, F : List(Dpol) ) ==
        empty?(F) => h
        h0 : Dpol := monomial(leadingCoefficient h, degree h)
        while (h := reductum h) ~= 0 repeat
           hred := redPo(h, F)
           h := hred.poly
           h0 := (hred.mult)*h0 + monomial(leadingCoefficient(h), degree h)
        h0

             --- concat ordered critical pair lists D1 and D2

   updatD(D1 : List(critPair), D2 : List(critPair)) ==
      empty?(D1) => D2
      empty?(D2) => D1
      dl1 := first(D1)
      dl2 := first(D2)
      critpOrder(dl1, dl2) => cons(dl1, updatD(D1.rest, D2))
      cons(dl2, updatD(D1, D2.rest))

            --- crit M - true, if lexp#2 multiple of lexp#1

   critM(e1 : Expon, e2 : Expon) ==
         en : Union(Expon, "failed")
         (en := subtractIfCan(e2, e1)) case Expon

   --- concat F and h and erase polys in F with lexp
   --- beeing multiple of lexp h
   --- Can do this because such poly can be recoverd
   --- from critical pair

   updatF(h : Dpol, deg : NNI, F : List(sugarPol)) ==
       empty?(F) => [[deg, h]]
       f1:= first(F)
       critM(degree(h), degree(f1.pol))  => updatF(h, deg, rest(F))
       cons(f1, updatF(h, deg, rest(F)))



)abbrev domain SKSMP SparseMultivariateSkewPolynomial
++ Description:
++   SparseMultivariateSkewPolynomial(R, Var, sigma, delta) defines
++   a mutivariate Ore ring over R in variables from V. \spad{sigma(v)}
++   gives automorphism of R corresponding to variable \spad{v} and
++   \spad{delta(v)} gives corresponding derivative.
++
SparseMultivariateSkewPolynomial(R : Ring, Var : OrderedSet,
      sigma : Var -> Automorphism R , delta : Var -> (R -> R)
    ) : Exports == Implementation where
  Exports ==> MultivariateSkewPolynomialCategory(R, IndexedExponents(Var),
                                                 Var) with
    D : Var -> %
      ++ D(v) returns operator corresponding to derivative with
      ++ respect to v in R.
    if Var has variable : Symbol -> Var then
        Delta : Symbol -> %
          ++ Delta(s) returns operator corresponding to derivative with
          ++ respect to s in R.

  Smp ==> SparseMultivariatePolynomial(R, Var)
  Implementation ==> Smp add

    D(v : Var) : % == monomial(1, v, 1)

    if Var has with (variable : Symbol -> Union(Var, "failed")) then

        Delta(s : Symbol) : % ==
            vu := variable(s)
            vu case "failed" => error "missing variable"
            D(vu@Var)

    Sup := SparseUnivariatePolynomial(%)
    Supp := SparseUnivariatePolynomial(Smp)

    -- Fake domain so that we can use UnivariateSkewPolynomialCategoryOps
    id : Automorphism(%) := morphism((x +-> x)@(% -> %))
    Upol := SparseUnivariateSkewPolynomial(%, id, (x : %) : % +-> 0)

    lift_map(f : R -> R) : % -> % ==
        x +-> map(f, x)

    lift_morphism(m : Automorphism R) : Automorphism % ==
        morphism(lift_map((x : R) : R +-> m(x)),
                 lift_map((x : R) : R +-> inv(m)(x)))

    import from UnivariateSkewPolynomialCategoryOps(%, Upol)

    x : % * y : % ==
        ground?(x) => retract(x)@R*y
        v := mainVariable(x)::Var
        xu := univariate(x pretend Smp, v)$Smp pretend Upol
        yu := univariate(y pretend Smp, v)$Smp pretend Upol
        ru := times(xu, yu, lift_morphism(sigma(v)), lift_map(delta(v)))
        multivariate(ru pretend Supp, v)$Smp pretend %

    x : % * r : R == x*(r::%)

    if R has IntegralDomain then
        x : % exquo y : % ==
            ground?(y) => x exquo retract(y)@R
            ground?(x) => "failed"
            maxd := totalDegree(x)
            ly := leadingMonomial(y) pretend Smp
            res : % := 0
            while x ~= 0 repeat
                lx1 := leadingMonomial(x)
                totalDegree(lx1) > maxd => return "failed"
                lx := lx1 pretend Smp
                cu := lx exquo ly
                cu case "failed" => return "failed"
                cc := (cu@Smp) pretend %
                res := res + cc
                x := x - cc*y
            res

    x : % ^ n : PositiveInteger ==
        res := x
        for i in 2..n repeat res := res*x
        res

    x : % ^ n : NonNegativeInteger ==
        n = 0 => 1
        x^(qcoerce(n)@PositiveInteger)

    coerce(x : %) : OutputForm ==
        ground?(x) => (retract(x)@R)::OutputForm
        v := mainVariable(x)::Var
        xu := univariate(x pretend Smp, v) pretend Upol
        outputForm(xu, sub('D::OutputForm, v::OutputForm))



)abbrev package PDOHLP PartialDifferentialOperatorHelper
PartialDifferentialOperatorHelper(R : PartialDifferentialRing(Var),
                                  Var : OrderedSet) : with
    id_map : Var -> Automorphism(R)
    diff_map : Var -> (R -> R)

 == add

    id_map(v) == 1$Automorphism(R)

    diff_map(v) == (x : R) : R +-> differentiate(x, v)



)abbrev domain PDO PartialDifferentialOperator
++ Description:
++   PartialDifferentialOperator(R, V) defines a ring of partial
++   differential operators in variables from V and with coefficients
++   in a partial differential ring R.  Multiplication of operators
++   corresponds to composition of operators.
++
PartialDifferentialOperator(R : PartialDifferentialRing(Var),
    Var : OrderedSet) : Exports == Implementation where
  HLP ==> PartialDifferentialOperatorHelper(R, Var)
  Exports ==> MultivariateSkewPolynomialCategory(R,
                  IndexedExponents(Var), Var) with
    D : Var -> %
      ++ D(v) returns the operator corresponding to derivative with
      ++ respect to v in R.
    adjoint : % -> %
      ++ adjoint(p) returns the adjoint of operator p.
    if R has LeftOreRing then
        SolvableSkewPolynomialCategory(R, IndexedExponents(Var))
  Implementation ==> SparseMultivariateSkewPolynomial(R, Var, id_map$HLP,
                        diff_map$HLP) add

    Smp ==> SparseMultivariatePolynomial(R, Var)
    Sup := SparseUnivariatePolynomial(%)

    adjoint(x : %) : % ==
        ground?(x) => x
        v := mainVariable(x)::Var
        xu := univariate(x pretend Smp, v)$Smp pretend Sup
        res : % := 0
        xu := map(adjoint, xu)$Sup
        while xu ~= 0 repeat
            d := degree(xu)
            sign := (odd?(d) => -1 ; 1)
            res := res + sign*monomial(1, v, d)*leadingCoefficient(xu)
            xu := reductum(xu)
        res

)abbrev package OREMAT OrePolynomialMatrixOperations
OrePolynomialMatrixOperations(F, LO) : Exports == Implementation where
  F : Field
  LO :  UnivariateSkewPolynomialCategory(F)
  M   ==> Matrix LO
  V   ==> Vector F
  Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F)
  Param_Rec_V ==> Record(ratpart : V, coeffs : Vector F)
  FPL ==> Record(particular : List Param_Rec_F, basis : List F)
  VPL ==> Record(particular : List Param_Rec_V, basis : List V)

  Exports ==> with

    rowEchelon : M -> M
      ++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
    rowEchelon : (M, F, List V) -> Record(mat : M, vecs : List V)
      ++ \spad{rowEchelon(m, c, lv)} returns [m2, lv2] such that
      ++ m2 is the row echelon form of the matrix m. lv2
      ++ is transformed lv using c as parameter to \spad{apply}.
    solve : (M, F, List V, (LO, List F) -> FPL) -> Union(VPL, "failed")
      ++ solve(m, c, lv, solf) returns "failed" of
      ++ [[vp1, ..., vp_m], [b1, ..., bl]]
      ++ such that bi-s are basis of solutions of homogeneous system
      ++ m bi = 0.  Each vpi = [r, [c1, ..., cn]] is a particular solution
      ++ of a parametric matrix equation m r = \sum ci vi where
      ++ lv = [v1, ..., vn].  solf is scalar solver, c is a
      ++ parameter to \spad{apply} (needed for action of LO on F).
      ++ "failed" means that system is underdetermined.

  Implementation ==> add

    import from LinearCombinationUtilities(F, SparseUnivariatePolynomial F)

    rowEchelon(m : M) == rowEchelon(m, 0, []).mat

    rowEchelon(m, cc, lv) ==
        mm := copy m
        lw := [copy(vv) for vv in lv]
        nc := ncols(mm)
        nr := nrows(mm)
        for w in lw repeat
            if #w ~= nc then
                error "rowEchelon: elements of lv must have size = ncols(m)"
        i : Integer := 1
        for j in 1..nc repeat
            if i > nr then break
            pivk : Integer := -1
            dpiv : Integer := -1
            mij : LO
            for k in i..nr repeat
                if not zero?(mkj := mm(k, j)) and (dp1 := degree(mkj);
                     (dpiv = -1 or dp1 < dpiv)) then
                    mij := mkj
                    pivk := k
                    dpiv := dp1
            dpiv = - 1 => iterate
            mm := swapRows!(mm, i, pivk)
            for w in lw repeat
                swap!(w, i, pivk)
            for k in (i + 1)..nr | not zero?(mkj := mm(k, j)) repeat
                (mij, c, d, u, v) := right_ext_ext_GCD(mij, mkj)
                for k1 in (j + 1)..nc repeat
                    el1 := qelt(mm, i, k1)
                    el2 := qelt(mm, k, k1)
                    qsetelt!(mm, i, k1, c*el1 + d*el2)
                    qsetelt!(mm, k, k1, u*el1 + v*el2)
                for w in lw repeat
                    elf1 := qelt(w, i)
                    elf2 := qelt(w, k)
                    qsetelt!(w, i, apply(c, cc, elf1) + apply(d, cc, elf2))
                    qsetelt!(w, k, apply(u, cc, elf1) + apply(v, cc, elf2))
                qsetelt!(mm, i, j, mij)
                qsetelt!(mm, k, j, 0)
            for k in 1..(i - 1) repeat
                zero?(mkj := mm(k, j)) => iterate
                (q, r) := rightDivide(mkj, mij)
                qsetelt!(mm, k, j, r)
                for k1 in (j + 1)..nc repeat
                    qsetelt!(mm, k, k1, qelt(mm, k, k1) - q*qelt(mm, i, k1))
                for w in lw repeat
                    qsetelt!(w, k, qelt(w, k) - apply(q, cc, qelt(w, i)))
            i := i + 1
        [mm, lw]

    triangular_solve(m : M, cc : F, lv : List V,
                     solf : (LO, List F) -> FPL) : Union(VPL, "failed") ==
        nc := ncols(m)
        nr := nrows(m)
        nc > nr => "failed"
        n := #lv
        bvl : List V := []
        cb := [new(n, 0)$Vector(F) for i in 1..n]
        cba := [new(nc, 0)$Vector(F) for bv in cb]
        for i in 1..n for bv in cb repeat
            bv(i) := 1
        for i in nc..1 by -1 repeat
            m(i, i) = 0 => return "failed"
            lvi0 := [v(i) for v in lv]
            lvi1 := [lin_comb(bv, lvi0) for bv in cb]
            nlvi := #lvi1
            nbas := #bvl
            lvi1 := concat(lvi1, [0 for bv in bvl])
            lvi2 : List(F) := []
            for bv in concat(cba, bvl) for vvi in lvi1 repeat
                for j in i+1..nc repeat
                    vvi := vvi - apply(m(i, j), cc, bv(j))
                lvi2 := cons(vvi, lvi2)
            lvi := reverse!(lvi2)
            resi := solf(m(i, i), lvi)
            nbvl : List V := []
            for be in resi.basis repeat
                bv := new(nc, 0)
                bv(i) := be
                nbvl := cons(bv, nbvl)
            prl := resi.particular
            nrl := #prl
            ncvl := [pr.coeffs for pr in prl]
            nsl1 := [pr.ratpart for pr in prl]
            cm1 := matrix([parts(ncv) for ncv in ncvl])$Matrix(F)
            cm2 := rowEchelon(horizConcat(cm1, scalarMatrix(nrl, 1)))
            nsl2 : List(F) := []
            noff := nlvi + nbas + 1
            for j in 1..nrl repeat
                ss : F := 0
                for k in noff.. for slk in nsl1 repeat
                    ss := ss + cm2(j, k)*slk
                nsl2 := cons(ss, nsl2)
            nrl1 : NonNegativeInteger := 0
            for j in nrl..1 by -1 for slj in nsl2 repeat
                for k in 1..nlvi repeat
                    if cm2(j, k) ~= 0 then
                        nrl1 := j
                        break
                nrl1 ~= 0 => break
                nbcv := new(nbas, 0)$Vector(F)
                for k in 1..nbas repeat
                    nbcv(k) := cm2(j, k + nlvi)
                bv := lin_comb(nbcv, bvl)
                bv(i) := slj
                nbvl := cons(bv, nbvl)
            nsl2 := reverse!(nsl2)
            ncb : List(Vector(F)) := []
            ncba : List(Vector(F)) := []
            for j in 1..nrl1 for slj in nsl2 repeat
                nbcv1 := new(nbas, 0)$Vector(F)
                for k in 1..nbas repeat
                    nbcv1(k) := cm2(j, k + nlvi)
                nbcv2 := new(nlvi, 0)$Vector(F)
                for k in 1..nlvi repeat
                    nbcv2(k) := cm2(j, k)
                bv := lin_comb(nbcv2, cb)
                ncb := cons(bv, ncb)
                na := lin_comb(nbcv2, cba)
                na := lin_comb!(nbcv1, na, bvl)
                na(i) := slj
                ncba := cons(na, ncba)
            cb := reverse!(ncb)
            cba := reverse!(ncba)
            bvl := nbvl
        [[[ba, bv] for ba in cba for bv in cb], bvl]

    solve(m, cc, lv, solf) ==
        rec := rowEchelon(m, cc, lv)
        triangular_solve(rec.mat, cc, rec.vecs, solf)