axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / MODOP.spad

)abbrev domain MODOP ModuleOperator
++ Author: Manuel Bronstein
++ Date Created: 15 May 1990
++ Date Last Updated: 17 June 1993
++ Description:
++ Algebra of ADDITIVE operators on a module.

ModuleOperator(R,M) : SIG == CODE where
  R : Ring
  M : LeftModule(R)

  O    ==> OutputForm
  OP   ==> BasicOperator
  FG   ==> FreeGroup OP
  RM   ==> Record(coef:R, monom:FG)
  TERM ==> List RM
  FAB  ==> FreeAbelianGroup TERM
  OPADJ   ==> "%opAdjoint"
  OPEVAL  ==> "%opEval"
  INVEVAL ==> "%invEval"

  SIG ==> Join(Ring, RetractableTo R, RetractableTo OP,
                   Eltable(M, M)) with

    if R has CharacteristicZero then CharacteristicZero

    if R has CharacteristicNonZero then CharacteristicNonZero

    if R has CommutativeRing then

      Algebra(R)

      adjoint : $ -> $
        ++ adjoint(op) returns the adjoint of the operator \spad{op}.

      adjoint : ($, $) -> $
        ++ adjoint(op1, op2) sets the adjoint of op1 to be op2.
        ++ op1 must be a basic operator

      conjug : R -> R
        ++ conjug(x)should be local but conditional

    evaluate : ($, M -> M) -> $
      ++ evaluate(f, u +-> g u) attaches the map g to f.
      ++ f must be a basic operator
      ++ g MUST be additive, \spad{g(a + b) = g(a) + g(b)} for
      ++ any \spad{a}, \spad{b} in M.
      ++ This implies that \spad{g(n a) = n g(a)} for
      ++ any \spad{a} in M and integer \spad{n > 0}.

    evaluateInverse : ($, M -> M) -> $
      ++ evaluateInverse(x,f) is not documented

    "**" : (OP, Integer) -> $
      ++ op**n is not documented

    "**" : ($, Integer) -> $
      ++ op**n is not documented

    opeval : (OP, M) -> M
      ++ opeval should be local but conditional

    makeop : (R, FG) -> $
      ++ makeop should be local but conditional

  CODE ==> FAB add

    import NoneFunctions1($)
    import BasicOperatorFunctions1(M)

    Rep := FAB

    inv      : TERM -> $
    termeval : (TERM, M) -> M
    rmeval   : (RM, M) -> M
    monomeval: (FG, M) -> M
    opInvEval: (OP, M) -> M
    mkop     : (R, FG) -> $
    termprod0: (Integer, TERM, TERM) -> $
    termprod : (Integer, TERM, TERM) -> TERM
    termcopy : TERM -> TERM
    trm2O    : (Integer, TERM) -> O
    term2O   : TERM -> O
    rm2O     : (R, FG) -> O
    nocopy   : OP -> $

    1                   == makeop(1, 1)

    coerce(n:Integer):$ == n::R::$

    coerce(r:R):$       == (zero? r => 0; makeop(r, 1))

    coerce(op:OP):$     == nocopy copy op

    nocopy(op:OP):$     == makeop(1, op::FG)

    elt(x:$, r:M)       == +/[t.exp * termeval(t.gen, r) for t in terms x]

    rmeval(t, r)        == t.coef * monomeval(t.monom, r)

    termcopy t          == [[rm.coef, rm.monom] for rm in t]

    characteristic()    == characteristic()$R

    mkop(r, fg)         == [[r, fg]$RM]$TERM :: $

    evaluate(f, g)   == nocopy setProperty(retract(f)@OP,OPEVAL,g pretend None)

    if R has OrderedSet then

      makeop(r, fg) == (r >= 0 => mkop(r, fg); - mkop(-r, fg))

    else

      makeop(r, fg) == mkop(r, fg)

    inv(t:TERM):$ ==
      empty? t => 1
      c := first(t).coef
      m := first(t).monom
      inv(rest t) * makeop(1, inv m) * (recip(c)::R::$)

    x:$ ** i:Integer ==
      i = 0 => 1
      i > 0 => expt(x,i pretend PositiveInteger)$RepeatedSquaring($)
      (inv(retract(x)@TERM)) ** (-i)

    evaluateInverse(f, g) ==
      nocopy setProperty(retract(f)@OP, INVEVAL, g pretend None)

    coerce(x:$):O ==
      zero? x => (0$R)::O
      reduce(_+, [trm2O(t.exp, t.gen) for t in terms x])$List(O)

    trm2O(c, t) ==
      (c = 1) => term2O t
      c = -1 => - term2O t
      c::O * term2O t

    term2O t ==
      reduce(_*, [rm2O(rm.coef, rm.monom) for rm in t])$List(O)

    rm2O(c, m) ==
      (c = 1) => m::O
      (m = 1) => c::O
      c::O * m::O

    x:$ * y:$ ==
      +/[ +/[termprod0(t.exp * s.exp, t.gen, s.gen) for s in terms y]
          for t in terms x]

    termprod0(n, x, y) ==
      n >= 0 => termprod(n, x, y)::$
      - (termprod(-n, x, y)::$)

    termprod(n, x, y) ==
      lc := first(xx := termcopy x)
      lc.coef := n * lc.coef
      rm := last xx
      ((first(y).coef) = 1) =>
        rm.monom := rm.monom * first(y).monom
        concat_!(xx, termcopy rest y)
      ((rm.monom) = 1) =>
        rm.coef := rm.coef * first(y).coef
        rm.monom := first(y).monom
        concat_!(xx, termcopy rest y)
      concat_!(xx, termcopy y)

    if M has ExpressionSpace then

      opeval(op, r) ==
        (func := property(op, OPEVAL)) case "failed" => kernel(op, r)
        ((func::None) pretend (M -> M)) r

    else

      opeval(op, r) ==
        (func := property(op, OPEVAL)) case "failed" =>
          error "eval: operator has no evaluation function"
        ((func::None) pretend (M -> M)) r

    opInvEval(op, r) ==
      (func := property(op, INVEVAL)) case "failed" =>
         error "eval: operator has no inverse evaluation function"
      ((func::None) pretend (M -> M)) r

    termeval(t, r)  ==
      for rm in reverse t repeat r := rmeval(rm, r)
      r

    monomeval(m, r) ==
      for rec in reverse_! factors m repeat
        e := rec.exp
        g := rec.gen
        e > 0 =>
          for i in 1..e repeat r := opeval(g, r)
        e < 0 =>
          for i in 1..(-e) repeat r := opInvEval(g, r)
      r

    recip x ==
      (r := retractIfCan(x)@Union(R, "failed")) case "failed" => "failed"
      (r1 := recip(r::R)) case "failed" => "failed"
      r1::R::$

    retractIfCan(x:$):Union(R, "failed") ==
      (r:= retractIfCan(x)@Union(TERM,"failed")) case "failed" => "failed"
      empty?(t := r::TERM) => 0$R
      empty? rest t =>
        rm := first t
        (rm.monom = 1) => rm.coef
        "failed"
      "failed"

    retractIfCan(x:$):Union(OP, "failed") ==
      (r:= retractIfCan(x)@Union(TERM,"failed")) case "failed" => "failed"
      empty?(t := r::TERM) => "failed"
      empty? rest t =>
        rm := first t
        (rm.coef = 1) => retractIfCan(rm.monom)
        "failed"
      "failed"

    if R has CommutativeRing then
      termadj  : TERM -> $
      rmadj    : RM -> $
      monomadj : FG -> $
      opadj    : OP -> $

      r:R * x:$        == r::$ * x

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

      adjoint x        == +/[t.exp * termadj(t.gen) for t in terms x]

      rmadj t          == conjug(t.coef) * monomadj(t.monom)

      adjoint(op, adj) == nocopy setProperty(retract(op)@OP, OPADJ, adj::None)

      termadj t ==
        ans:$ := 1
        for rm in t repeat ans := rmadj(rm) * ans
        ans

      monomadj m ==
        ans:$ := 1
        for rec in factors m repeat ans := (opadj(rec.gen) ** rec.exp) * ans
        ans

      opadj op ==
        (adj := property(op, OPADJ)) case "failed" =>
           error "adjoint: operator does not have a defined adjoint"
        (adj::None) pretend $

      if R has conjugate:R -> R then

        conjug r == conjugate r else conjug r == r