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

)abbrev domain EXPR Expression
++ Author: Manuel Bronstein
++ Date Created: 19 July 1988
++ Date Last Updated: October 1993 (P.Gianni), February 1995 (MB)
++ Description: 
++ Top-level mathematical expressions involving symbolic functions.

Expression(R) : SIG == CODE where
  R : OrderedSet

  Q   ==> Fraction Integer
  K   ==> Kernel %
  MP  ==> SparseMultivariatePolynomial(R, K)
  AF  ==> AlgebraicFunction(R, %)
  EF  ==> ElementaryFunction(R, %)
  CF  ==> CombinatorialFunction(R, %)
  LF  ==> LiouvillianFunction(R, %)
  AN  ==> AlgebraicNumber
  KAN ==> Kernel AN
  FSF ==> FunctionalSpecialFunction(R, %)
  ESD ==> ExpressionSpace_&(%)
  FSD ==> FunctionSpace_&(%, R)
  SYMBOL ==> "%symbol"
  ALGOP  ==> "%alg"
  POWER  ==> "%power"::Symbol
  SUP    ==> SparseUnivariatePolynomial

  SIG ==> FunctionSpace R with

    if R has IntegralDomain then

      AlgebraicallyClosedFunctionSpace R

      TranscendentalFunctionCategory

      CombinatorialOpsCategory

      LiouvillianFunctionCategory

      SpecialFunctionCategory

      reduce: % -> %
        ++ reduce(f) simplifies all the unreduced algebraic quantities
        ++ present in f by applying their defining relations.

      number?: % -> Boolean
        ++ number?(f) tests if f is rational

      simplifyPower: (%,Integer) -> %
        ++ simplifyPower(f,n) is not documented

      if R has GcdDomain then

        factorPolynomial : SUP  % -> Factored SUP %
          ++ factorPolynomial(p) is not documented

        squareFreePolynomial : SUP % -> Factored SUP %
          ++ squareFreePolynomial(p) is not documented

      if R has RetractableTo Integer then RetractableTo AN

  CODE ==> add

    import KernelFunctions2(R, %)

    retNotUnit     : % -> R

    retNotUnitIfCan: % -> Union(R, "failed")

    belong? op == true

    retNotUnit x ==
      (u := constantIfCan(k := retract(x)@K)) case R => u::R
      error "Not retractable"

    retNotUnitIfCan x ==
      (r := retractIfCan(x)@Union(K,"failed")) case "failed" => "failed"
      constantIfCan(r::K)

    if R has IntegralDomain then

      reduc  : (%, List Kernel %) -> %

      commonk   : (%, %) -> List K

      commonk0  : (List K, List K) -> List K

      toprat    : % -> %

      algkernels: List K -> List K

      evl       : (MP, K, SparseUnivariatePolynomial %) -> Fraction MP

      evl0      : (MP, K) -> SparseUnivariatePolynomial Fraction MP

      Rep := Fraction MP

      0                == 0$Rep

      1                == 1$Rep

      one? x           == (x = 1)$Rep

      zero? x          == zero?(x)$Rep

      - x:%            == -$Rep x

      n:Integer * x:%  == n *$Rep x

      coerce(n:Integer) ==  coerce(n)$Rep@Rep::%

      x:% * y:%        == reduc(x *$Rep y, commonk(x, y))

      x:% + y:%        == reduc(x +$Rep y, commonk(x, y))

      (x:% - y:%):%    == reduc(x -$Rep y, commonk(x, y))

      x:% / y:%        == reduc(x /$Rep y, commonk(x, y))

      number?(x:%):Boolean ==
        if R has RetractableTo(Integer) then
          ground?(x) or ((retractIfCan(x)@Union(Q,"failed")) case Q)
        else
          ground?(x)

      simplifyPower(x:%,n:Integer):% ==
        k : List K := kernels x
        is?(x,POWER) =>
          -- Look for a power of a number in case we can do a simplification
          args : List % := argument first k
          not(#args = 2) => error "Too many arguments to **"
          number?(args.1) =>
             reduc((args.1) **$Rep n, algkernels kernels (args.1))**(args.2)
          (first args)**(n*second(args))
        reduc(x **$Rep n, algkernels k)

      x:% ** n:NonNegativeInteger ==
        n = 0 => 1$%
        n = 1 => x
        simplifyPower(numerator x,n pretend Integer) / 
           simplifyPower(denominator x,n pretend Integer)

      x:% ** n:Integer ==
        n = 0 => 1$%
        n = 1 => x
        n = -1 => 1/x
        simplifyPower(numerator x,n) / 
           simplifyPower(denominator x,n)

      x:% ** n:PositiveInteger == 
        n = 1 => x
        simplifyPower(numerator x,n pretend Integer) / 
           simplifyPower(denominator x,n pretend Integer)

      x:% < y:%        == x <$Rep y

      x:% = y:%        == x =$Rep y

      numer x          == numer(x)$Rep

      denom x          == denom(x)$Rep

      coerce(p:MP):%   == coerce(p)$Rep

      reduce x         == reduc(x, algkernels kernels x)

      commonk(x, y)    == commonk0(algkernels kernels x, algkernels kernels y)

      algkernels l     == select_!(x +-> has?(operator x, ALGOP), l)

      toprat f == ratDenom(f,algkernels kernels f)$AlgebraicManipulations(R, %)

      x:MP / y:MP ==
       reduc(x /$Rep y,commonk0(algkernels variables x,algkernels variables y))

-- since we use the reduction from FRAC SMP which asssumes that the
-- variables are independent, we must remove algebraic from the denominators
      reducedSystem(m:Matrix %):Matrix(R) ==
        mm:Matrix(MP) := reducedSystem(map(toprat, m))$Rep
        reducedSystem(mm)$MP

-- since we use the reduction from FRAC SMP which asssumes that the
-- variables are independent, we must remove algebraic from the denominators
      reducedSystem(m:Matrix %, v:Vector %):
       Record(mat:Matrix R, vec:Vector R) ==
        r:Record(mat:Matrix MP, vec:Vector MP) :=
          reducedSystem(map(toprat, m), map(toprat, v))$Rep
        reducedSystem(r.mat, r.vec)$MP

-- The result MUST be left sorted deepest first   MB 3/90
      commonk0(x, y) ==
        ans := empty()$List(K)
        for k in reverse_! x repeat if member?(k, y) then ans := concat(k, ans)
        ans

      rootOf(x:SparseUnivariatePolynomial %, v:Symbol) == rootOf(x,v)$AF

      pi()                      == pi()$EF

      exp x                     == exp(x)$EF

      log x                     == log(x)$EF

      sin x                     == sin(x)$EF

      cos x                     == cos(x)$EF

      tan x                     == tan(x)$EF

      cot x                     == cot(x)$EF

      sec x                     == sec(x)$EF

      csc x                     == csc(x)$EF

      asin x                    == asin(x)$EF

      acos x                    == acos(x)$EF

      atan x                    == atan(x)$EF

      acot x                    == acot(x)$EF

      asec x                    == asec(x)$EF

      acsc x                    == acsc(x)$EF

      sinh x                    == sinh(x)$EF

      cosh x                    == cosh(x)$EF

      tanh x                    == tanh(x)$EF

      coth x                    == coth(x)$EF

      sech x                    == sech(x)$EF

      csch x                    == csch(x)$EF

      asinh x                   == asinh(x)$EF

      acosh x                   == acosh(x)$EF

      atanh x                   == atanh(x)$EF

      acoth x                   == acoth(x)$EF

      asech x                   == asech(x)$EF

      acsch x                   == acsch(x)$EF

      abs x                     == abs(x)$FSF

      Gamma x                   == Gamma(x)$FSF

      Gamma(a, x)               == Gamma(a, x)$FSF

      Beta(x,y)                 == Beta(x,y)$FSF

      digamma x                 == digamma(x)$FSF

      polygamma(k,x)            == polygamma(k,x)$FSF

      besselJ(v,x)              == besselJ(v,x)$FSF

      besselY(v,x)              == besselY(v,x)$FSF

      besselI(v,x)              == besselI(v,x)$FSF

      besselK(v,x)              == besselK(v,x)$FSF

      airyAi x                  == airyAi(x)$FSF

      airyBi x                  == airyBi(x)$FSF

      x:% ** y:%                == x **$CF y

      factorial x               == factorial(x)$CF

      binomial(n, m)            == binomial(n, m)$CF

      permutation(n, m)         == permutation(n, m)$CF

      factorials x              == factorials(x)$CF

      factorials(x, n)          == factorials(x, n)$CF

      summation(x:%, n:Symbol)           == summation(x, n)$CF

      summation(x:%, s:SegmentBinding %) == summation(x, s)$CF

      product(x:%, n:Symbol)             == product(x, n)$CF

      product(x:%, s:SegmentBinding %)   == product(x, s)$CF

      erf x                              == erf(x)$LF

      Ei x                               == Ei(x)$LF

      Si x                               == Si(x)$LF

      Ci x                               == Ci(x)$LF

      li x                               == li(x)$LF

      dilog x                            == dilog(x)$LF

      fresnelS x                         == fresnelS(x)$LF

      fresnelC x                         == fresnelC(x)$LF

      integral(x:%, n:Symbol)            == integral(x, n)$LF

      integral(x:%, s:SegmentBinding %)  == integral(x, s)$LF

      operator op ==
        belong?(op)$AF  => operator(op)$AF
        belong?(op)$EF  => operator(op)$EF
        belong?(op)$CF  => operator(op)$CF
        belong?(op)$LF  => operator(op)$LF
        belong?(op)$FSF => operator(op)$FSF
        belong?(op)$FSD => operator(op)$FSD
        belong?(op)$ESD => operator(op)$ESD
        nullary? op and has?(op, SYMBOL) => operator(kernel(name op)$K)
        (n := arity op) case "failed" => operator name op
        operator(name op, n::NonNegativeInteger)

      reduc(x, l) ==
        for k in l repeat
          p := minPoly k
          x := evl(numer x, k, p) /$Rep evl(denom x, k, p)
        x

      evl0(p, k) ==
        numer univariate(p::Fraction(MP),
                     k)$PolynomialCategoryQuotientFunctions(IndexedExponents K,
                                                            K,R,MP,Fraction MP)

      -- uses some operations from Rep instead of % in order not to
      -- reduce recursively during those operations.
      evl(p, k, m) ==
        degree(p, k) < degree m => p::Fraction(MP)
        (((evl0(p, k) pretend SparseUnivariatePolynomial($)) rem m)
          pretend SparseUnivariatePolynomial Fraction MP) (k::MP::Fraction(MP))

      if R has GcdDomain then

        noalg?: SUP % -> Boolean

        noalg? p ==
          while p ^= 0 repeat
            not empty? algkernels kernels leadingCoefficient p => return false
            p := reductum p
          true

        gcdPolynomial(p:SUP %, q:SUP %) ==
          noalg? p and noalg? q => gcdPolynomial(p, q)$Rep
          gcdPolynomial(p, q)$GcdDomain_&(%)

        factorPolynomial(x:SUP %) : Factored SUP % ==
          uf:= factor(x pretend SUP(Rep))$SupFractionFactorizer(
                                          IndexedExponents K,K,R,MP)
          uf pretend Factored SUP %

        squareFreePolynomial(x:SUP %) : Factored SUP % ==
          uf:= squareFree(x pretend SUP(Rep))$SupFractionFactorizer(
                                          IndexedExponents K,K,R,MP)
          uf pretend Factored SUP %

      if R is AN then
        -- this is to force the coercion R -> EXPR R to be used
        -- instead of the coercioon AN -> EXPR R which loops.
        -- simpler looking code will fail! MB 10/91
        coerce(x:AN):% == (monomial(x, 0$IndexedExponents(K))$MP)::%

      if (R has RetractableTo Integer) then

        x:% ** r:Q                           == x **$AF r

        minPoly k                            == minPoly(k)$AF

        definingPolynomial x                 == definingPolynomial(x)$AF

        retract(x:%):Q                       == retract(x)$Rep

        retractIfCan(x:%):Union(Q, "failed") == retractIfCan(x)$Rep

        if not(R is AN) then

          k2expr  : KAN -> %

          smp2expr: SparseMultivariatePolynomial(Integer, KAN) -> %

          R2AN    : R  -> Union(AN, "failed")

          k2an    : K  -> Union(AN, "failed")

          smp2an  : MP -> Union(AN, "failed")


          coerce(x:AN):% == smp2expr(numer x) / smp2expr(denom x)

          k2expr k       == map(x+->x::%, k)$ExpressionSpaceFunctions2(AN, %)

          smp2expr p ==
            map(k2expr,x+->x::%,p)_
              $PolynomialCategoryLifting(IndexedExponents KAN,
                   KAN, Integer, SparseMultivariatePolynomial(Integer, KAN), %)

          retractIfCan(x:%):Union(AN, "failed") ==
            ((n:= smp2an numer x) case AN) and ((d:= smp2an denom x) case AN)
                 => (n::AN) / (d::AN)
            "failed"

          R2AN r ==
            (u := retractIfCan(r::%)@Union(Q, "failed")) case Q => u::Q::AN
            "failed"

          k2an k ==
            not(belong?(op := operator k)$AN) => "failed"
            arg:List(AN) := empty()
            for x in argument k repeat
              if (a := retractIfCan(x)@Union(AN, "failed")) case "failed" then
                return "failed"
              else arg := concat(a::AN, arg)
            (operator(op)$AN) reverse_!(arg)

          smp2an p ==
            (x1 := mainVariable p) case "failed" => R2AN leadingCoefficient p
            up := univariate(p, k := x1::K)
            (t  := k2an k) case "failed" => "failed"
            ans:AN := 0
            while not ground? up repeat
              (c:=smp2an leadingCoefficient up) case "failed" _
                => return "failed"
              ans := ans + (c::AN) * (t::AN) ** (degree up)
              up  := reductum up
            (c := smp2an leadingCoefficient up) case "failed" => "failed"
            ans + c::AN

      if R has ConvertibleTo InputForm then

        convert(x:%):InputForm == convert(x)$Rep

        import MakeUnaryCompiledFunction(%, %, %)

        eval(f:%, op: BasicOperator, g:%, x:Symbol):% == 
          eval(f,[op],[g],x)

        eval(f:%, ls:List BasicOperator, lg:List %, x:Symbol) ==
          -- handle subsrcipted symbols by renaming -> eval -> renaming back
          llsym:List List Symbol:=[variables g for g in lg]
          lsym:List Symbol:= removeDuplicates concat llsym
          lsd:List Symbol:=select (scripted?,lsym)
          empty? lsd=> eval(f,ls,[compiledFunction(g, x) for g in lg])
          ns:List Symbol:=[new()$Symbol for i in lsd]
          lforwardSubs:List Equation % := _
            [(i::%)= (j::%) for i in lsd for j in ns]
          lbackwardSubs:List Equation % := _
            [(j::%)= (i::%) for i in lsd for j in ns]
          nlg:List % :=[subst(g,lforwardSubs) for g in lg]
          res:% :=eval(f, ls, [compiledFunction(g, x) for g in nlg])
          subst(res,lbackwardSubs)

      if R has PatternMatchable Integer then

        patternMatch(x:%, p:Pattern Integer,
         l:PatternMatchResult(Integer, %)) ==
          patternMatch(x, p, l)$PatternMatchFunctionSpace(Integer, R, %)

      if R has PatternMatchable Float then

        patternMatch(x:%, p:Pattern Float,
         l:PatternMatchResult(Float, %)) ==
          patternMatch(x, p, l)$PatternMatchFunctionSpace(Float, R, %)

    else  -- R is not an integral domain

      operator op ==
        belong?(op)$FSD => operator(op)$FSD
        belong?(op)$ESD => operator(op)$ESD
        nullary? op and has?(op, SYMBOL) => operator(kernel(name op)$K)
        (n := arity op) case "failed" => operator name op
        operator(name op, n::NonNegativeInteger)

      if R has Ring then

        Rep := MP

        0              == 0$Rep

        1              == 1$Rep

        - x:%          == -$Rep x

        n:Integer *x:% == n *$Rep x

        x:% * y:%      == x *$Rep y

        x:% + y:%      == x +$Rep y

        x:% = y:%      == x =$Rep y

        x:% < y:%      == x <$Rep y

        numer x        == x@Rep

        coerce(p:MP):% == p

        reducedSystem(m:Matrix %):Matrix(R) ==
          reducedSystem(m)$Rep

        reducedSystem(m:Matrix %, v:Vector %):
         Record(mat:Matrix R, vec:Vector R) ==
          reducedSystem(m, v)$Rep

        if R has ConvertibleTo InputForm then

          convert(x:%):InputForm == convert(x)$Rep

        if R has PatternMatchable Integer then

          kintmatch: (K,Pattern Integer,PatternMatchResult(Integer,Rep))
                                     -> PatternMatchResult(Integer, Rep)

          kintmatch(k, p, l) ==
            patternMatch(k, p, l pretend PatternMatchResult(Integer, %)
              )$PatternMatchKernel(Integer, %)
                pretend PatternMatchResult(Integer, Rep)

          patternMatch(x:%, p:Pattern Integer,
           l:PatternMatchResult(Integer, %)) ==
            patternMatch(x@Rep, p,
                         l pretend PatternMatchResult(Integer, Rep),
                          kintmatch
                           )$PatternMatchPolynomialCategory(Integer,
                            IndexedExponents K, K, R, Rep)
                              pretend PatternMatchResult(Integer, %)

        if R has PatternMatchable Float then

          kfltmatch: (K, Pattern Float, PatternMatchResult(Float, Rep))
                                     -> PatternMatchResult(Float, Rep)

          kfltmatch(k, p, l) ==
            patternMatch(k, p, l pretend PatternMatchResult(Float, %)
              )$PatternMatchKernel(Float, %)
                pretend PatternMatchResult(Float, Rep)

          patternMatch(x:%, p:Pattern Float,
           l:PatternMatchResult(Float, %)) ==
            patternMatch(x@Rep, p,
                         l pretend PatternMatchResult(Float, Rep),
                          kfltmatch
                           )$PatternMatchPolynomialCategory(Float,
                            IndexedExponents K, K, R, Rep)
                              pretend PatternMatchResult(Float, %)

      else   -- R is not even a ring

        if R has AbelianMonoid then

          import ListToMap(K, %)

          kereval        : (K, List K, List %) -> %

          subeval        : (K, List K, List %) -> %

          Rep := FreeAbelianGroup K

          0              == 0$Rep

          x:% + y:%      == x +$Rep y

          x:% = y:%      == x =$Rep y

          x:% < y:%      == x <$Rep y

          coerce(k:K):%  == coerce(k)$Rep

          kernels x      == [f.gen for f in terms x]

          coerce(x:R):%  == (zero? x => 0; constantKernel(x)::%)

          retract(x:%):R == (zero? x => 0; retNotUnit x)

          coerce(x:%):OutputForm == coerce(x)$Rep

          kereval(k, lk, lv) == 
           match(lk, lv, k, (x2:K):% +-> map(x1 +-> eval(x1, lk, lv), x2))

          subeval(k, lk, lv) ==
            match(lk, lv, k,
             (x:K):% +->
               kernel(operator x, [subst(a, lk, lv) for a in argument x]))

          isPlus x ==
            empty?(l := terms x) or empty? rest l => "failed"
            [t.exp *$Rep t.gen for t in l]$List(%)

          isMult x ==
            empty?(l := terms x) or not empty? rest l => "failed"
            t := first l
            [t.exp, t.gen]

          eval(x:%, lk:List K, lv:List %) ==
            _+/[t.exp * kereval(t.gen, lk, lv) for t in terms x]

          subst(x:%, lk:List K, lv:List %) ==
            _+/[t.exp * subeval(t.gen, lk, lv) for t in terms x]

          retractIfCan(x:%):Union(R, "failed") ==
            zero? x => 0
            retNotUnitIfCan x

          if R has AbelianGroup then -(x:%) == -$Rep x

        else   -- R is nothing

            import ListToMap(K, %)

            Rep := K

            x:% < y:%      == x <$Rep y

            x:% = y:%      == x =$Rep y

            coerce(k:K):%  == k

            kernels x      == [x pretend K]

            coerce(x:R):%  == constantKernel x

            retract(x:%):R == retNotUnit x

            retractIfCan(x:%):Union(R, "failed") == retNotUnitIfCan x

            coerce(x:%):OutputForm               == coerce(x)$Rep

            eval(x:%, lk:List K, lv:List %) ==
              match(lk, lv, x pretend K, 
               (x1:K):% +-> map(x2 +-> eval(x2, lk, lv), x1))

            subst(x, lk, lv) ==
              match(lk, lv, x pretend K,
               (x1:K):% +-> 
                 kernel(operator x1, [subst(a, lk, lv) for a in argument x1]))

            if R has ConvertibleTo InputForm then
              convert(x:%):InputForm == convert(x)$Rep