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

)abbrev package STTF StreamTranscendentalFunctions
++ Author: William Burge, Clifton J. Williamson
++ Date Created: 1986
++ Date Last Updated: 6 March 1995
++ Description:
++ StreamTranscendentalFunctions implements transcendental functions on
++ Taylor series, where a Taylor series is represented by a stream of
++ its coefficients.

StreamTranscendentalFunctions(Coef) : SIG == CODE where
  Coef : Algebra Fraction Integer

  L   ==> List
  I   ==> Integer
  RN  ==> Fraction Integer
  SG  ==> String
  ST  ==> Stream Coef
  STT ==> StreamTaylorSeriesOperations Coef
  YS  ==> Y$ParadoxicalCombinatorsForStreams(Coef)

  SIG ==> with

    exp : ST -> ST
      ++ exp(st) computes the exponential of a power series st.

    log : ST -> ST
      ++ log(st) computes the log of a power series.

    "**" : (ST,ST) -> ST
      ++ st1 ** st2 computes the power of a power series st1 by another
      ++ power series st2.

    sincos  : ST -> Record(sin:ST, cos:ST)
      ++ sincos(st) returns a record containing the sine and cosine
      ++ of a power series st.

    sin : ST -> ST
      ++ sin(st) computes sine of a power series st.

    cos : ST -> ST
      ++ cos(st) computes cosine of a power series st.

    tan : ST -> ST
      ++ tan(st) computes tangent of a power series st.

    cot : ST -> ST
      ++ cot(st) computes cotangent of a power series st.

    sec : ST -> ST
      ++ sec(st) computes secant of a power series st.

    csc : ST -> ST
      ++ csc(st) computes cosecant of a power series st.

    asin : ST -> ST
      ++ asin(st) computes arcsine of a power series st.

    acos : ST -> ST
      ++ acos(st) computes arccosine of a power series st.

    atan : ST -> ST
      ++ atan(st) computes arctangent of a power series st.

    acot : ST -> ST
      ++ acot(st) computes arccotangent of a power series st.

    asec : ST -> ST
      ++ asec(st) computes arcsecant of a power series st.

    acsc : ST -> ST
      ++ acsc(st) computes arccosecant of a power series st.

    sinhcosh: ST -> Record(sinh:ST, cosh:ST)
      ++ sinhcosh(st) returns a record containing
      ++ the hyperbolic sine and cosine
      ++ of a power series st.
    sinh    : ST -> ST
      ++ sinh(st) computes the hyperbolic sine of a power series st.
    cosh    : ST -> ST
      ++ cosh(st) computes the hyperbolic cosine of a power series st.
    tanh    : ST -> ST
      ++ tanh(st) computes the hyperbolic tangent of a power series st.
    coth    : ST -> ST
      ++ coth(st) computes the hyperbolic cotangent of a power series st.
    sech    : ST -> ST
      ++ sech(st) computes the hyperbolic secant of a power series st.
    csch    : ST -> ST
      ++ csch(st) computes the hyperbolic cosecant of a power series st.
    asinh   : ST -> ST
      ++ asinh(st) computes the inverse hyperbolic sine of a power series st.
    acosh   : ST -> ST
      ++ acosh(st) computes the inverse hyperbolic cosine
      ++ of a power series st.
    atanh   : ST -> ST
      ++ atanh(st) computes the inverse hyperbolic tangent
      ++ of a power series st.
    acoth   : ST -> ST
      ++ acoth(st) computes the inverse hyperbolic
      ++ cotangent of a power series st.
    asech   : ST -> ST
      ++ asech(st) computes the inverse hyperbolic secant of a
      ++ power series st.
    acsch   : ST -> ST
      ++ acsch(st) computes the inverse hyperbolic
      ++ cosecant of a power series st.

  CODE ==> add

    import StreamTaylorSeriesOperations Coef

    TRANSFCN : Boolean := Coef has TranscendentalFunctionCategory

--% Error Reporting

    TRCONST : SG := "series expansion involves transcendental constants"

    NPOWERS : SG := "series expansion has terms of negative degree"

    FPOWERS : SG := "series expansion has terms of fractional degree"

    MAYFPOW : SG := "series expansion may have terms of fractional degree"

    LOGS : SG := "series expansion has logarithmic term"

    NPOWLOG : SG :=
       "series expansion has terms of negative degree or logarithmic term"

    FPOWLOG : SG :=
       "series expansion has terms of fractional degree or logarithmic term"

    NOTINV : SG := "leading coefficient not invertible"

--% Exponentials and Logarithms

    expre:(Coef,ST,ST) -> ST
    expre(r,e,dx) == lazyIntegrate(r,e*dx)

    exp z ==
      empty? z => 1 :: ST
      (coef := frst z) = 0 => YS(y +-> expre(1,y,deriv z))
      TRANSFCN => YS(y +-> expre(exp coef,y,deriv z))
      error concat("exp: ",TRCONST)

    log z ==
      empty? z => error "log: constant coefficient should not be 0"
      (coef := frst z) = 0 => error "log: constant coefficient should not be 0"
      coef = 1 => lazyIntegrate(0,deriv z/z)
      TRANSFCN => lazyIntegrate(log coef,deriv z/z)
      error concat("log: ",TRCONST)

    z1:ST ** z2:ST == exp(z2 * log z1)

--% Trigonometric Functions

    sincosre:(Coef,Coef,L ST,ST,Coef) -> L ST
    sincosre(rs,rc,sc,dx,sign) ==
      [lazyIntegrate(rs,(second sc)*dx),lazyIntegrate(rc,sign*(first sc)*dx)]

    sincos z ==
      empty? z => [0 :: ST,1 :: ST]
      l :=
        (coef := frst z) = 0 => YS(y +-> sincosre(0,1,y,deriv z,-1),2)
        TRANSFCN => YS(y +-> sincosre(sin coef,cos coef,y,deriv z,-1),2)
        error concat("sincos: ",TRCONST)
      [first l,second l]

    sin z == sincos(z).sin
    cos z == sincos(z).cos

    tanre:(Coef,ST,ST,Coef) -> ST
    tanre(r,t,dx,sign) == lazyIntegrate(r,((1 :: ST) + sign*t*t)*dx)

    tan z ==
      empty? z => 0 :: ST
      (coef := frst z) = 0 => YS(y +-> tanre(0,y,deriv z,1))
      TRANSFCN => YS(y +-> tanre(tan coef,y,deriv z,1))
      error concat("tan: ",TRCONST)

    cotre:(Coef,ST,ST) -> ST
    cotre(r,t,dx) == lazyIntegrate(r,-((1 :: ST) + t*t)*dx)

    cot z ==
      empty? z => error "cot: cot(0) is undefined"
      (coef := frst z) = 0 => error concat("cot: ",NPOWERS)
      TRANSFCN => YS(y +-> cotre(cot coef,y,deriv z))
      error concat("cot: ",TRCONST)

    sec z ==
      empty? z => 1 :: ST
      frst z = 0 => recip(cos z) :: ST
      TRANSFCN =>
        cosz := cos z
        first cosz = 0 => error concat("sec: ",NPOWERS)
        recip(cosz) :: ST
      error concat("sec: ",TRCONST)

    csc z ==
      empty? z => error "csc: csc(0) is undefined"
      TRANSFCN =>
        sinz := sin z
        first sinz = 0 => error concat("csc: ",NPOWERS)
        recip(sinz) :: ST
      error concat("csc: ",TRCONST)

    orderOrFailed : ST -> Union(I,"failed")
    orderOrFailed x ==
    -- returns the order of x or "failed"
    -- if -1 is returned, the series is identically zero
      for n in 0..1000 repeat
        empty? x => return -1
        not zero? frst x => return n :: I
        x := rst x
      "failed"

    asin z ==
      empty? z => 0 :: ST
      (coef := frst z) = 0 =>
        integrate(0,powern(-1/2,(1 :: ST) - z*z) * (deriv z))
      TRANSFCN =>
        coef = 1 or coef = -1 =>
          x := (1 :: ST) - z*z
          -- compute order of 'x'
          (ord := orderOrFailed x) case "failed" =>
            error concat("asin: ",MAYFPOW)
          (order := ord :: I) = -1 => return asin(coef) :: ST
          odd? order => error concat("asin: ",FPOWERS)
          squirt := powern(1/2,x)
          (quot := (deriv z) exquo squirt) case "failed" =>
             error concat("asin: ",NOTINV)
          integrate(asin coef,quot :: ST)
        integrate(asin coef,powern(-1/2,(1 :: ST) - z*z) * (deriv z))
      error concat("asin: ",TRCONST)

    acos z ==
      empty? z =>
        TRANSFCN => acos(0)$Coef :: ST
        error concat("acos: ",TRCONST)
      TRANSFCN =>
        coef := frst z
        coef = 1 or coef = -1 =>
          x := (1 :: ST) - z*z
          -- compute order of 'x'
          (ord := orderOrFailed x) case "failed" =>
            error concat("acos: ",MAYFPOW)
          (order := ord :: I) = -1 => return acos(coef) :: ST
          odd? order => error concat("acos: ",FPOWERS)
          squirt := powern(1/2,x)
          (quot := (-deriv z) exquo squirt) case "failed" =>
             error concat("acos: ",NOTINV)
          integrate(acos coef,quot :: ST)
        integrate(acos coef,-powern(-1/2,(1 :: ST) - z*z) * (deriv z))
      error concat("acos: ",TRCONST)

    atan z ==
      empty? z => 0 :: ST
      (coef := frst z) = 0 =>
        integrate(0,(recip((1 :: ST) + z*z) :: ST) * (deriv z))
      TRANSFCN =>
        (y := recip((1 :: ST) + z*z)) case "failed" =>
          error concat("atan: ",LOGS)
        integrate(atan coef,(y :: ST) * (deriv z))
      error concat("atan: ",TRCONST)

    acot z ==
      empty? z =>
        TRANSFCN => acot(0)$Coef :: ST
        error concat("acot: ",TRCONST)
      TRANSFCN =>
        (y := recip((1 :: ST) + z*z)) case "failed" =>
          error concat("acot: ",LOGS)
        integrate(acot frst z,-(y :: ST) * (deriv z))
      error concat("acot: ",TRCONST)

    asec z ==
      empty? z => error "asec: constant coefficient should not be 0"
      TRANSFCN =>
        (coef := frst z) = 0 =>
          error "asec: constant coefficient should not be 0"
        coef = 1 or coef = -1 =>
          x := z*z - (1 :: ST)
          -- compute order of 'x'
          (ord := orderOrFailed x) case "failed" =>
            error concat("asec: ",MAYFPOW)
          (order := ord :: I) = -1 => return asec(coef) :: ST
          odd? order => error concat("asec: ",FPOWERS)
          squirt := powern(1/2,x)
          (quot := (deriv z) exquo squirt) case "failed" =>
            error concat("asec: ",NOTINV)
          (quot2 := (quot :: ST) exquo z) case "failed" =>
            error concat("asec: ",NOTINV)
          integrate(asec coef,quot2 :: ST)
        integrate(asec coef,(powern(-1/2,z*z-(1::ST))*(deriv z)) / z)
      error concat("asec: ",TRCONST)

    acsc z ==
      empty? z => error "acsc: constant coefficient should not be zero"
      TRANSFCN =>
        (coef := frst z) = 0 =>
          error "acsc: constant coefficient should not be zero"
        coef = 1 or coef = -1 =>
          x := z*z - (1 :: ST)
          -- compute order of 'x'
          (ord := orderOrFailed x) case "failed" =>
            error concat("acsc: ",MAYFPOW)
          (order := ord :: I) = -1 => return acsc(coef) :: ST
          odd? order => error concat("acsc: ",FPOWERS)
          squirt := powern(1/2,x)
          (quot := (-deriv z) exquo squirt) case "failed" =>
            error concat("acsc: ",NOTINV)
          (quot2 := (quot :: ST) exquo z) case "failed" =>
            error concat("acsc: ",NOTINV)
          integrate(acsc coef,quot2 :: ST)
        integrate(acsc coef,-(powern(-1/2,z*z-(1::ST))*(deriv z)) / z)
      error concat("acsc: ",TRCONST)

--% Hyperbolic Trigonometric Functions

    sinhcosh z ==
      empty? z => [0 :: ST,1 :: ST]
      l :=
        (coef := frst z) = 0 => YS(y +-> sincosre(0,1,y,deriv z,1),2)
        TRANSFCN => YS(y +-> sincosre(sinh coef,cosh coef,y,deriv z,1),2)
        error concat("sinhcosh: ",TRCONST)
      [first l,second l]

    sinh z == sinhcosh(z).sinh

    cosh z == sinhcosh(z).cosh

    tanh z ==
      empty? z => 0 :: ST
      (coef := frst z) = 0 => YS(y +-> tanre(0,y,deriv z,-1))
      TRANSFCN => YS(y +-> tanre(tanh coef,y,deriv z,-1))
      error concat("tanh: ",TRCONST)

    coth z ==
      tanhz := tanh z
      empty? tanhz => error "coth: coth(0) is undefined"
      (frst tanhz) = 0 => error concat("coth: ",NPOWERS)
      recip(tanhz) :: ST

    sech z ==
      coshz := cosh z
      (empty? coshz) or (frst coshz = 0) => error concat("sech: ",NPOWERS)
      recip(coshz) :: ST

    csch z ==
      sinhz := sinh z
      (empty? sinhz) or (frst sinhz = 0) => error concat("csch: ",NPOWERS)
      recip(sinhz) :: ST

    asinh z ==
      empty? z => 0 :: ST
      (coef := frst z) = 0 => log(z + powern(1/2,(1 :: ST) + z*z))
      TRANSFCN =>
        x := (1 :: ST) + z*z
        -- compute order of 'x', in case coefficient(z,0) = +- %i
        (ord := orderOrFailed x) case "failed" =>
          error concat("asinh: ",MAYFPOW)
        (order := ord :: I) = -1 => return asinh(coef) :: ST
        odd? order => error concat("asinh: ",FPOWERS)
        -- the argument to 'log' must have a non-zero constant term
        log(z + powern(1/2,x))
      error concat("asinh: ",TRCONST)

    acosh z ==
      empty? z =>
        TRANSFCN => acosh(0)$Coef :: ST
        error concat("acosh: ",TRCONST)
      TRANSFCN =>
        coef := frst z
        coef = 1 or coef = -1 =>
          x := z*z - (1 :: ST)
          -- compute order of 'x'
          (ord := orderOrFailed x) case "failed" =>
            error concat("acosh: ",MAYFPOW)
          (order := ord :: I) = -1 => return acosh(coef) :: ST
          odd? order => error concat("acosh: ",FPOWERS)
          -- the argument to 'log' must have a non-zero constant term
          log(z + powern(1/2,x))
        log(z + powern(1/2,z*z - (1 :: ST)))
      error concat("acosh: ",TRCONST)

    atanh z ==
      empty? z => 0 :: ST
      (coef := frst z) = 0 =>
        (inv(2::RN)::Coef) * log(((1 :: ST) + z)/((1 :: ST) - z))
      TRANSFCN =>
        coef = 1 or coef = -1 => error concat("atanh: ",LOGS)
        (inv(2::RN)::Coef) * log(((1 :: ST) + z)/((1 :: ST) - z))
      error concat("atanh: ",TRCONST)

    acoth z ==
      empty? z =>
        TRANSFCN => acoth(0)$Coef :: ST
        error concat("acoth: ",TRCONST)
      TRANSFCN =>
        frst z = 1 or frst z = -1 => error concat("acoth: ",LOGS)
        (inv(2::RN)::Coef) * log((z + (1 :: ST))/(z - (1 :: ST)))
      error concat("acoth: ",TRCONST)

    asech z ==
      empty? z => error "asech: asech(0) is undefined"
      TRANSFCN =>
        (coef := frst z) = 0 => error concat("asech: ",NPOWLOG)
        coef = 1 or coef = -1 =>
          x := (1 :: ST) - z*z
          -- compute order of 'x'
          (ord := orderOrFailed x) case "failed" =>
            error concat("asech: ",MAYFPOW)
          (order := ord :: I) = -1 => return asech(coef) :: ST
          odd? order => error concat("asech: ",FPOWERS)
          log(((1 :: ST) + powern(1/2,x))/z)
        log(((1 :: ST) + powern(1/2,(1 :: ST) - z*z))/z)
      error concat("asech: ",TRCONST)

    acsch z ==
      empty? z => error "acsch: acsch(0) is undefined"
      TRANSFCN =>
        frst z = 0 => error concat("acsch: ",NPOWLOG)
        x := z*z + (1 :: ST)
        -- compute order of 'x'
        (ord := orderOrFailed x) case "failed" =>
          error concat("acsc: ",MAYFPOW)
        (order := ord :: I) = -1 => return acsch(frst z) :: ST
        odd? order => error concat("acsch: ",FPOWERS)
        log(((1 :: ST) + powern(1/2,x))/z)
      error concat("acsch: ",TRCONST)