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

)abbrev package EFSTRUC ElementaryFunctionStructurePackage
++ Author: Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 16 August 1995
++ Description:
++ ElementaryFunctionStructurePackage provides functions to test the
++ algebraic independence of various elementary functions, using the
++ Risch structure theorem (real and complex versions).
++ It also provides transformations on elementary functions
++ which are not considered simplifications.

ElementaryFunctionStructurePackage(R,F) : SIG == CODE where
  R : Join(IntegralDomain, OrderedSet, RetractableTo Integer,
           LinearlyExplicitRingOver Integer)
  F : Join(AlgebraicallyClosedField, TranscendentalFunctionCategory,
           FunctionSpace R)

  B   ==> Boolean
  N   ==> NonNegativeInteger
  Z   ==> Integer
  Q   ==> Fraction Z
  SY  ==> Symbol
  K   ==> Kernel F
  UP  ==> SparseUnivariatePolynomial F
  SMP ==> SparseMultivariatePolynomial(R, K)
  REC ==> Record(func:F, kers: List K, vals:List F)
  U   ==> Union(vec:Vector Q, func:F, fail: Boolean)
  POWER ==> "%power"::SY
  NTHR  ==> "nthRoot"::SY

  SIG ==> with

    normalize : F -> F
      ++ normalize(f) rewrites \spad{f} using the least possible number of
      ++ real algebraically independent kernels.

    normalize : (F, SY) -> F
      ++ normalize(f, x) rewrites \spad{f} using the least possible number of
      ++ real algebraically independent kernels involving \spad{x}.

    rischNormalize : (F, SY) -> REC
      ++ rischNormalize(f, x) returns \spad{[g, [k1,...,kn], [h1,...,hn]]}
      ++ such that \spad{g = normalize(f, x)} and each \spad{ki} was
      ++ rewritten as \spad{hi} during the normalization.

    realElementary : F -> F
      ++ realElementary(f) rewrites \spad{f} in terms of the 4 fundamental real
      ++ transcendental elementary functions: \spad{log, exp, tan, atan}.

    realElementary : (F, SY) -> F
      ++ realElementary(f,x) rewrites the kernels of \spad{f} involving
      ++ \spad{x} in terms of the 4 fundamental real
      ++ transcendental elementary functions: \spad{log, exp, tan, atan}.

    validExponential : (List K, F, SY) -> Union(F, "failed")
      ++ validExponential([k1,...,kn],f,x) returns \spad{g} if \spad{exp(f)=g}
      ++ and \spad{g} involves only \spad{k1...kn}, and "failed" otherwise.

    rootNormalize : (F, K) -> F
      ++ rootNormalize(f, k) returns \spad{f} rewriting either \spad{k} which
      ++ must be an nth-root in terms of radicals already in \spad{f}, or some
      ++ radicals in \spad{f} in terms of \spad{k}.

    tanQ : (Q, F) -> F
      ++ tanQ(q,a) is a local function with a conditional implementation.

  CODE ==> add

    import TangentExpansions F
    import IntegrationTools(R, F)
    import IntegerLinearDependence F
    import AlgebraicManipulations(R, F)
    import InnerCommonDenominator(Z, Q, Vector Z, Vector Q)

    k2Elem             : (K, List SY) -> F
    realElem           : (F, List SY) -> F
    smpElem            : (SMP, List SY) -> F
    deprel             : (List K, K, SY) -> U
    rootDep            : (List K, K)     -> U
    qdeprel            : (List F, F)     -> U
    factdeprel         : (List K, K)     -> U
    toR                : (List K, F) -> List K
    toY                : List K -> List F
    toZ                : List K -> List F
    toU                : List K -> List F
    toV                : List K -> List F
    ktoY               : K  -> F
    ktoZ               : K  -> F
    ktoU               : K  -> F
    ktoV               : K  -> F
    gdCoef?            : (Q, Vector Q) -> Boolean
    goodCoef           : (Vector Q, List K, SY) ->
                                 Union(Record(index:Z, ker:K), "failed")
    tanRN              : (Q, K) -> F
    localnorm          : F -> F
    rooteval           : (F, List K, K, Q) -> REC
    logeval            : (F, List K, K, Vector Q) -> REC
    expeval            : (F, List K, K, Vector Q) -> REC
    taneval            : (F, List K, K, Vector Q) -> REC
    ataneval           : (F, List K, K, Vector Q) -> REC
    depeval            : (F, List K, K, Vector Q) -> REC
    expnosimp          : (F, List K, K, Vector Q, List F, F) -> REC
    tannosimp          : (F, List K, K, Vector Q, List F, F) -> REC
    rtNormalize        : F -> F
    rootNormalize0     : F -> REC
    rootKernelNormalize: (F, List K, K) -> Union(REC, "failed")
    tanSum             : (F, List F) -> F

    comb?     := F has CombinatorialOpsCategory
    mpiover2:F := pi()$F / (-2::F)

    realElem(f, l)       == smpElem(numer f, l) / smpElem(denom f, l)

    realElementary(f, x) == realElem(f, [x])

    realElementary f     == realElem(f, variables f)

    toY ker              == [func for k in ker | (func := ktoY k) ^= 0]

    toZ ker              == [func for k in ker | (func := ktoZ k) ^= 0]

    toU ker              == [func for k in ker | (func := ktoU k) ^= 0]

    toV ker              == [func for k in ker | (func := ktoV k) ^= 0]

    rtNormalize f        == rootNormalize0(f).func

    toR(ker, x) == select(s+->is?(s, NTHR) and first argument(s) = x, ker)

    if R has GcdDomain then

      tanQ(c, x) ==
        tanNa(rootSimp zeroOf tanAn(x, denom(c)::PositiveInteger), numer c)

    else

      tanQ(c, x) ==
        tanNa(zeroOf tanAn(x, denom(c)::PositiveInteger), numer c)

    -- tanSum(c, [a1,...,an]) returns f(c, a1,...,an) such that
    -- if ai = tan(ui) then f(c, a1,...,an) = tan(c + u1 + ... + un).
    -- MUST BE CAREFUL FOR WHEN c IS AN ODD MULTIPLE of pi/2
    tanSum(c, l) ==
      k := c / mpiover2        -- k = - 2 c / pi, check for odd integer
                               -- tan((2n+1) pi/2 x) = - 1 / tan x
      (r := retractIfCan(k)@Union(Z, "failed")) case Z and odd?(r::Z) =>
           - inv tanSum l
      tanSum concat(tan c, l)

    rootNormalize0 f ==
      ker := select_!(s+->is?(s, NTHR) and empty? variables first argument s,
                      tower f)$List(K)
      empty? ker => [f, empty(), empty()]
      (n := (#ker)::Z - 1) < 1 => [f, empty(), empty()]
      for i in 1..n for kk in rest ker repeat
        (u := rootKernelNormalize(f, first(ker, i), kk)) case REC =>
          rec := u::REC
          rn  := rootNormalize0(rec.func)
          return [rn.func, concat(rec.kers,rn.kers), concat(rec.vals, rn.vals)]
      [f, empty(), empty()]

    deprel(ker, k, x) ==
      is?(k, "log"::SY) or is?(k, "exp"::SY) =>
        qdeprel([differentiate(g, x) for g in toY ker],
                 differentiate(ktoY k, x))
      is?(k, "atan"::SY) or is?(k, "tan"::SY) =>
        qdeprel([differentiate(g, x) for g in toU ker],
                 differentiate(ktoU k, x))
      is?(k, NTHR) => rootDep(ker, k)
      comb? and is?(k, "factorial"::SY) =>
        factdeprel([x for x in ker | is?(x,"factorial"::SY) and x^=k],k)
      [true]

    ktoY k ==
      is?(k, "log"::SY) => k::F
      is?(k, "exp"::SY) => first argument k
      0

    ktoZ k ==
      is?(k, "log"::SY) => first argument k
      is?(k, "exp"::SY) => k::F
      0

    ktoU k ==
      is?(k, "atan"::SY) => k::F
      is?(k,  "tan"::SY) => first argument k
      0

    ktoV k ==
      is?(k,  "tan"::SY) => k::F
      is?(k, "atan"::SY) => first argument k
      0

    smpElem(p, l) ==
      map(x+->k2Elem(x, l), y+->y::F, p)_
       $PolynomialCategoryLifting(IndexedExponents K, K, R, SMP, F)

    k2Elem(k, l) ==
      ez, iez, tz2: F
      kf := k::F
      not(empty? l) and empty? [v for v in variables kf | member?(v, l)] => kf
      empty?(args :List F := [realElem(a, l) for a in argument k]) => kf
      z := first args
      is?(k, POWER)       => (zero? z => 0; exp(last(args) * log z))
      is?(k, "cot"::SY)   => inv tan z
      is?(k, "acot"::SY)  => atan inv z
      is?(k, "asin"::SY)  => atan(z / sqrt(1 - z**2))
      is?(k, "acos"::SY)  => atan(sqrt(1 - z**2) / z)
      is?(k, "asec"::SY)  => atan sqrt(1 - z**2)
      is?(k, "acsc"::SY)  => atan inv sqrt(1 - z**2)
      is?(k, "asinh"::SY) => log(sqrt(1 + z**2) + z)
      is?(k, "acosh"::SY) => log(sqrt(z**2 - 1) + z)
      is?(k, "atanh"::SY) => log((z + 1) / (1 - z)) / (2::F)
      is?(k, "acoth"::SY) => log((z + 1) / (z - 1)) / (2::F)
      is?(k, "asech"::SY) => log((inv z) + sqrt(inv(z**2) - 1))
      is?(k, "acsch"::SY) => log((inv z) + sqrt(1 + inv(z**2)))
      is?(k, "%paren"::SY) or is?(k, "%box"::SY) =>
        empty? rest args => z
        kf
      if has?(op := operator k, "htrig") then iez  := inv(ez  := exp z)
      is?(k, "sinh"::SY)  => (ez - iez) / (2::F)
      is?(k, "cosh"::SY)  => (ez + iez) / (2::F)
      is?(k, "tanh"::SY)  => (ez - iez) / (ez + iez)
      is?(k, "coth"::SY)  => (ez + iez) / (ez - iez)
      is?(k, "sech"::SY)  => 2 * inv(ez + iez)
      is?(k, "csch"::SY)  => 2 * inv(ez - iez)
      if has?(op, "trig") then tz2  := tan(z / (2::F))
      is?(k, "sin"::SY)   => 2 * tz2 / (1 + tz2**2)
      is?(k, "cos"::SY)   => (1 - tz2**2) / (1 + tz2**2)
      is?(k, "sec"::SY)   => (1 + tz2**2) / (1 - tz2**2)
      is?(k, "csc"::SY)   => (1 + tz2**2) / (2 * tz2)
      op args

    --The next 5 functions are used by normalize, once a relation is found

    depeval(f, lk, k, v) ==
      is?(k, "log"::SY)  => logeval(f, lk, k, v)
      is?(k, "exp"::SY)  => expeval(f, lk, k, v)
      is?(k, "tan"::SY)  => taneval(f, lk, k, v)
      is?(k, "atan"::SY) => ataneval(f, lk, k, v)
      is?(k, NTHR) => rooteval(f, lk, k, v(minIndex v))
      [f, empty(), empty()]

    rooteval(f, lk, k, n) ==
      nv := nthRoot(x := first argument k, m := retract(n)@Z)
      l  := [r for r in concat(k, toR(lk, x)) |
             retract(second argument r)@Z ^= m]
      lv := [nv ** (n / (retract(second argument r)@Z::Q)) for r in l]
      [eval(f, l, lv), l, lv]

    ataneval(f, lk, k, v) ==
      w := first argument k
      s := tanSum [tanQ(qelt(v,i), x)
                   for i in minIndex v .. maxIndex v for x in toV lk]
      g := +/[qelt(v, i) * x for i in minIndex v .. maxIndex v for x in toU lk]
      h:F :=
        zero?(d := 1 + s * w) => mpiover2
        atan((w - s) / d)
      g := g + h
      [eval(f, [k], [g]), [k], [g]]

    gdCoef?(c, v) ==
      for i in minIndex v .. maxIndex v repeat
        retractIfCan(qelt(v, i) / c)@Union(Z, "failed") case "failed" =>
          return false
      true

    goodCoef(v, l, s) ==
      for i in minIndex v .. maxIndex v for k in l repeat
        is?(k, s) and
           ((r:=recip(qelt(v,i))) case Q) and
            (retractIfCan(r::Q)@Union(Z, "failed") case Z)
              and gdCoef?(qelt(v, i), v) => return([i, k])
      "failed"

    taneval(f, lk, k, v) ==
      u := first argument k
      fns := toU lk
      c := u - +/[qelt(v, i)*x for i in minIndex v .. maxIndex v for x in fns]
      (rec := goodCoef(v, lk, "tan"::SY)) case "failed" =>
          tannosimp(f, lk, k, v, fns, c)
      v0 := retract(inv qelt(v, rec.index))@Z
      lv := [qelt(v, i) for i in minIndex v .. maxIndex v |
                                                 i ^= rec.index]$List(Q)
      l  := [kk for kk in lk | kk ^= rec.ker]
      g := tanSum(-v0 * c, concat(tanNa(k::F, v0),
           [tanNa(x, - retract(a * v0)@Z) for a in lv for x in toV l]))
      [eval(f, [rec.ker], [g]), [rec.ker], [g]]

    tannosimp(f, lk, k, v, fns, c) ==
      every?(x+->is?(x, "tan"::SY), lk) =>
        dd := (d := (cd := splitDenominator v).den)::F
        newt := [tan(u / dd) for u in fns]$List(F)
        newtan := [tanNa(t, d) for t in newt]$List(F)
        h := tanSum(c, [tanNa(t, qelt(cd.num, i))
                        for i in minIndex v .. maxIndex v for t in newt])
        lk := concat(k, lk)
        newtan := concat(h, newtan)
        [eval(f, lk, newtan), lk, newtan]
      h := tanSum(c, [tanQ(qelt(v, i), x)
                      for i in minIndex v .. maxIndex v for x in toV lk])
      [eval(f, [k], [h]), [k], [h]]

    expnosimp(f, lk, k, v, fns, g) ==
      every?(x+->is?(x, "exp"::SY), lk) =>
        dd := (d := (cd := splitDenominator v).den)::F
        newe := [exp(y / dd) for y in fns]$List(F)
        newexp := [e ** d for e in newe]$List(F)
        h := */[e ** qelt(cd.num, i)
                for i in minIndex v .. maxIndex v for e in newe] * g
        lk := concat(k, lk)
        newexp := concat(h, newexp)
        [eval(f, lk, newexp), lk, newexp]
      h := */[exp(y) ** qelt(v, i)
                for i in minIndex v .. maxIndex v for y in fns] * g
      [eval(f, [k], [h]), [k], [h]]

    logeval(f, lk, k, v) ==
      z := first argument k
      c := z / (*/[x**qelt(v, i)
                   for x in toZ lk for i in minIndex v .. maxIndex v])
      -- CHANGED log ktoZ x TO ktoY x 
      -- SINCE WE WANT log exp f TO BE REPLACED BY f.
      g := +/[qelt(v, i) * x
              for i in minIndex v .. maxIndex v for x in toY lk] + log c
      [eval(f, [k], [g]), [k], [g]]

    rischNormalize(f, v) ==
      empty?(ker := varselect(tower f, v)) => [f, empty(), empty()]
      first(ker) ^= kernel(v)@K => error "Cannot happen"
      ker := rest ker
      (n := (#ker)::Z - 1) < 1 => [f, empty(), empty()]
      for i in 1..n for kk in rest ker repeat
        klist := first(ker, i)
        -- NO EVALUATION ON AN EMPTY VECTOR, WILL CAUSE INFINITE LOOP
        (c := deprel(klist, kk, v)) case vec and not empty?(c.vec) =>
          rec := depeval(f, klist, kk, c.vec)
          rn  := rischNormalize(rec.func, v)
          return [rn.func,
                   concat(rec.kers, rn.kers), concat(rec.vals, rn.vals)]
        c case func =>
          rn := rischNormalize(eval(f, [kk], [c.func]), v)
          return [rn.func, concat(kk, rn.kers), concat(c.func, rn.vals)]
      [f, empty(), empty()]

    rootNormalize(f, k) ==
      (u := rootKernelNormalize(f, toR(tower f, first argument k), k))
         case "failed" => f
      (u::REC).func

    rootKernelNormalize(f, l, k) ==
      (c := rootDep(l, k)) case vec =>
        rooteval(f, l, k, (c.vec)(minIndex(c.vec)))
      "failed"

    localnorm f ==
      for x in variables f repeat
        f := rischNormalize(f, x).func
      f

    validExponential(twr, eta, x) ==
      (c := solveLinearlyOverQ(construct([differentiate(g, x)
         for g in (fns := toY twr)]$List(F))@Vector(F),
           differentiate(eta, x))) case "failed" => "failed"
      v := c::Vector(Q)
      g := eta - +/[qelt(v, i) * yy
                        for i in minIndex v .. maxIndex v for yy in fns]
      */[exp(yy) ** qelt(v, i)
                for i in minIndex v .. maxIndex v for yy in fns] * exp g

    rootDep(ker, k) ==
      empty?(ker := toR(ker, first argument k)) => [true]
      [new(1,lcm(retract(second argument k)@Z,
       "lcm"/[retract(second argument r)@Z for r in ker])::Q)$Vector(Q)]

    qdeprel(l, v) ==
      (u := solveLinearlyOverQ(construct(l)@Vector(F), v))
        case Vector(Q) => [u::Vector(Q)]
      [true]

    expeval(f, lk, k, v) ==
      y   := first argument k
      fns := toY lk
      g:= y - +/[qelt(v, i) * z for i in minIndex v .. maxIndex v for z in fns]
      (rec := goodCoef(v, lk, "exp"::SY)) case "failed" =>
        expnosimp(f, lk, k, v, fns, exp g)
      v0 := retract(inv qelt(v, rec.index))@Z
      lv := [qelt(v, i) for i in minIndex v .. maxIndex v |
                                                 i ^= rec.index]$List(Q)
      l  := [kk for kk in lk | kk ^= rec.ker]
      h :F := */[exp(z) ** (- retract(a * v0)@Z) for a in lv for z in toY l]
      h := h * exp(-v0 * g) * (k::F) ** v0
      [eval(f, [rec.ker], [h]), [rec.ker], [h]]

    if F has CombinatorialOpsCategory then

      normalize f == rtNormalize localnorm factorials realElementary f

      normalize(f, x) ==
        rtNormalize(rischNormalize(factorials(realElementary(f,x),x),x).func)

      factdeprel(l, k) ==
        ((r := retractIfCan(n := first argument k)@Union(Z, "failed"))
          case Z) and (r::Z > 0) => [factorial(r::Z)::F]
        for x in l repeat
          m := first argument x
          ((r := retractIfCan(n - m)@Union(Z, "failed")) case Z) and
            (r::Z > 0) => return([*/[(m + i::F) for i in 1..r] * x::F])
        [true]

    else

      normalize f     == rtNormalize localnorm realElementary f

      normalize(f, x)== rtNormalize(rischNormalize(realElementary(f,x),x).func)