--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2009, Gabriel Dos Reis.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

-- SPAD files for the functional world should be compiled in the
-- following order:
--
--   op  kl  fspace  algfunc  ELEMNTRY  expr
)abbrev package EF ElementaryFunction
++ Author: Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 10 April 1995
++ Keywords: elementary, function, logarithm, exponential.
++ Examples:  )r EF INPUT
++ Description: Provides elementary functions over an integral domain.
ElementaryFunction(R, F): Exports == Implementation where
  R: IntegralDomain
  F: Join(FunctionSpace R, RadicalCategory)

  B   ==> Boolean
  L   ==> List
  Z   ==> Integer
  OP  ==> BasicOperator
  K   ==> Kernel F
  INV ==> error "Invalid argument"

  Exports ==> with
    exp     : F -> F
	++ exp(x) applies the exponential operator to x
    log     : F -> F
	++ log(x) applies the logarithm operator to x
    sin     : F -> F
	++ sin(x) applies the sine operator to x
    cos     : F -> F
	++ cos(x) applies the cosine operator to x 
    tan     : F -> F
	++ tan(x) applies the tangent operator to x
    cot     : F -> F
	++ cot(x) applies the cotangent operator to x
    sec     : F -> F
	++ sec(x) applies the secant operator to x
    csc     : F -> F
	++ csc(x) applies the cosecant operator to x
    asin    : F -> F
	++ asin(x) applies the inverse sine operator to x 
    acos    : F -> F
	++ acos(x) applies the inverse cosine operator to x
    atan    : F -> F
	++ atan(x) applies the inverse tangent operator to x
    acot    : F -> F
	++ acot(x) applies the inverse cotangent operator to x
    asec    : F -> F
	++ asec(x) applies the inverse secant operator to x
    acsc    : F -> F
	++ acsc(x) applies the inverse cosecant operator to x
    sinh    : F -> F
	++ sinh(x) applies the hyperbolic sine operator to x 
    cosh    : F -> F
	++ cosh(x) applies the hyperbolic cosine operator to x
    tanh    : F -> F
	++ tanh(x) applies the hyperbolic tangent operator to x
    coth    : F -> F
	++ coth(x) applies the hyperbolic cotangent operator to x
    sech    : F -> F
	++ sech(x) applies the hyperbolic secant operator to x
    csch    : F -> F
	++ csch(x) applies the hyperbolic cosecant operator to x
    asinh   : F -> F
	++ asinh(x) applies the inverse hyperbolic sine operator to x
    acosh   : F -> F
	++ acosh(x) applies the inverse hyperbolic cosine operator to x
    atanh   : F -> F
	++ atanh(x) applies the inverse hyperbolic tangent operator to x
    acoth   : F -> F
	++ acoth(x) applies the inverse hyperbolic cotangent operator to x
    asech   : F -> F
	++ asech(x) applies the	inverse hyperbolic secant operator to x
    acsch   : F -> F
	++ acsch(x) applies the inverse hyperbolic cosecant operator to x
    pi      : () -> F
	++ pi() returns the pi operator
    belong? : OP -> Boolean
	++ belong?(p) returns true if operator p is elementary
    operator: OP -> OP
	++ operator(p) returns an elementary operator with the same symbol as p
    -- the following should be local, but are conditional
    iisqrt2   : () -> F
	++ iisqrt2() should be local but conditional
    iisqrt3   : () -> F
	++ iisqrt3() should be local but conditional
    iiexp     : F -> F
	++ iiexp(x) should be local but conditional
    iilog     : F -> F
	++ iilog(x) should be local but conditional
    iisin     : F -> F
	++ iisin(x) should be local but conditional
    iicos     : F -> F
	++ iicos(x) should be local but conditional
    iitan     : F -> F
	++ iitan(x) should be local but conditional
    iicot     : F -> F
	++ iicot(x) should be local but conditional
    iisec     : F -> F
	++ iisec(x) should be local but conditional
    iicsc     : F -> F
	++ iicsc(x) should be local but conditional
    iiasin    : F -> F
	++ iiasin(x) should be local but conditional
    iiacos    : F -> F
	++ iiacos(x) should be local but conditional
    iiatan    : F -> F
	++ iiatan(x) should be local but conditional
    iiacot    : F -> F
	++ iiacot(x) should be local but conditional
    iiasec    : F -> F
	++ iiasec(x) should be local but conditional
    iiacsc    : F -> F
	++ iiacsc(x) should be local but conditional
    iisinh    : F -> F
	++ iisinh(x) should be local but conditional
    iicosh    : F -> F
	++ iicosh(x) should be local but conditional
    iitanh    : F -> F
	++ iitanh(x) should be local but conditional
    iicoth    : F -> F
	++ iicoth(x) should be local but conditional
    iisech    : F -> F
	++ iisech(x) should be local but conditional
    iicsch    : F -> F
	++ iicsch(x) should be local but conditional
    iiasinh   : F -> F
	++ iiasinh(x) should be local but conditional
    iiacosh   : F -> F
	++ iiacosh(x) should be local but conditional
    iiatanh   : F -> F
	++ iiatanh(x) should be local but conditional
    iiacoth   : F -> F
	++ iiacoth(x) should be local but conditional
    iiasech   : F -> F
	++ iiasech(x) should be local but conditional
    iiacsch   : F -> F
	++ iiacsch(x) should be local but conditional
    specialTrigs:(F, L Record(func:F,pole:B)) -> Union(F, "failed")
	++ specialTrigs(x,l) should be local but conditional
    localReal?: F -> Boolean
	++ localReal?(x) should be local but conditional

  Implementation ==> add
    ipi      : List F -> F
    iexp     : F -> F
    ilog     : F -> F
    iiilog   : F -> F
    isin     : F -> F
    icos     : F -> F
    itan     : F -> F
    icot     : F -> F
    isec     : F -> F
    icsc     : F -> F
    iasin    : F -> F
    iacos    : F -> F
    iatan    : F -> F
    iacot    : F -> F
    iasec    : F -> F
    iacsc    : F -> F
    isinh    : F -> F
    icosh    : F -> F
    itanh    : F -> F
    icoth    : F -> F
    isech    : F -> F
    icsch    : F -> F
    iasinh   : F -> F
    iacosh   : F -> F
    iatanh   : F -> F
    iacoth   : F -> F
    iasech   : F -> F
    iacsch   : F -> F
    dropfun  : F -> F
    kernel   : F -> K
    posrem   :(Z, Z) -> Z
    iisqrt1  : () -> F
    valueOrPole : Record(func:F, pole:B) -> F

    oppi  := operator('pi)$CommonOperators
    oplog := operator('log)$CommonOperators
    opexp := operator('exp)$CommonOperators
    opsin := operator('sin)$CommonOperators
    opcos := operator('cos)$CommonOperators
    optan := operator('tan)$CommonOperators
    opcot := operator('cot)$CommonOperators
    opsec := operator('sec)$CommonOperators
    opcsc := operator('csc)$CommonOperators
    opasin := operator('asin)$CommonOperators
    opacos := operator('acos)$CommonOperators
    opatan := operator('atan)$CommonOperators
    opacot := operator('acot)$CommonOperators
    opasec := operator('asec)$CommonOperators
    opacsc := operator('acsc)$CommonOperators
    opsinh := operator('sinh)$CommonOperators
    opcosh := operator('cosh)$CommonOperators
    optanh := operator('tanh)$CommonOperators
    opcoth := operator('coth)$CommonOperators
    opsech := operator('sech)$CommonOperators
    opcsch := operator('csch)$CommonOperators
    opasinh := operator('asinh)$CommonOperators
    opacosh := operator('acosh)$CommonOperators
    opatanh := operator('atanh)$CommonOperators
    opacoth := operator('acoth)$CommonOperators
    opasech := operator('asech)$CommonOperators
    opacsch := operator('acsch)$CommonOperators

    -- Pi is a domain...
    Pie, isqrt1, isqrt2, isqrt3: F

    -- following code is conditionalized on arbitraryPrecesion to recompute in
    -- case user changes the precision

    if R has TranscendentalFunctionCategory then
      Pie := pi()$R :: F
    else
      Pie := kernel(oppi, nil()$List(F))

    if R has TranscendentalFunctionCategory and R has arbitraryPrecision then
      pi() == pi()$R :: F
    else
      pi() == Pie

    if R has imaginary: () -> R then
      isqrt1 := imaginary()$R :: F
    else isqrt1 := sqrt(-1::F)

    if R has RadicalCategory then
      isqrt2 := sqrt(2::R)::F
      isqrt3 := sqrt(3::R)::F
    else
      isqrt2 := sqrt(2::F)
      isqrt3 := sqrt(3::F)

    iisqrt1() == isqrt1
    if R has RadicalCategory and R has arbitraryPrecision then
      iisqrt2() == sqrt(2::R)::F
      iisqrt3() == sqrt(3::R)::F
    else
      iisqrt2() == isqrt2
      iisqrt3() == isqrt3

    ipi l == pi()
    log x == oplog x
    exp x == opexp x
    sin x == opsin x
    cos x == opcos x
    tan x == optan x
    cot x == opcot x
    sec x == opsec x
    csc x == opcsc x
    asin x == opasin x
    acos x == opacos x
    atan x == opatan x
    acot x == opacot x
    asec x == opasec x
    acsc x == opacsc x
    sinh x == opsinh x
    cosh x == opcosh x
    tanh x == optanh x
    coth x == opcoth x
    sech x == opsech x
    csch x == opcsch x
    asinh x == opasinh x
    acosh x == opacosh x
    atanh x == opatanh x
    acoth x == opacoth x
    asech x == opasech x
    acsch x == opacsch x
    kernel x == retract(x)@K

    posrem(n, m)    == (negative?(r := n rem m) => r + m; r)
    valueOrPole rec == (rec.pole => INV; rec.func)
    belong? op      == has?(op, 'elem)

    operator op ==
      is?(op,'pi)    => oppi
      is?(op,'log)   => oplog
      is?(op,'exp)   => opexp
      is?(op,'sin)   => opsin
      is?(op,'cos)   => opcos
      is?(op,'tan)   => optan
      is?(op,'cot)   => opcot
      is?(op,'sec)   => opsec
      is?(op,'csc)   => opcsc
      is?(op,'asin)  => opasin
      is?(op,'acos)  => opacos
      is?(op,'atan)  => opatan
      is?(op,'acot)  => opacot
      is?(op,'asec)  => opasec
      is?(op,'acsc)  => opacsc
      is?(op,'sinh)  => opsinh
      is?(op,'cosh)  => opcosh
      is?(op,'tanh)  => optanh
      is?(op,'coth)  => opcoth
      is?(op,'sech)  => opsech
      is?(op,'csch)  => opcsch
      is?(op,'asinh) => opasinh
      is?(op,'acosh) => opacosh
      is?(op,'atanh) => opatanh
      is?(op,'acoth) => opacoth
      is?(op,'asech) => opasech
      is?(op,'acsch) => opacsch
      error "Not an elementary operator"

    dropfun x ==
      ((k := retractIfCan(x)@Union(K, "failed")) case "failed") or
        empty?(argument(k::K)) => 0
      first argument(k::K)

    if R has RetractableTo Z then
      specialTrigs(x, values) ==
        (r := retractIfCan(y := x/pi())@Union(Fraction Z, "failed"))
          case "failed" => "failed"
        q := r::Fraction(Integer)
        m := minIndex values
        (n := retractIfCan(q)@Union(Z, "failed")) case Z =>
          even?(n::Z) => valueOrPole(values.m)
          valueOrPole(values.(m+1))
        (n := retractIfCan(2*q)@Union(Z, "failed")) case Z =>
          one?(s := posrem(n::Z, 4)) => valueOrPole(values.(m+2))
          valueOrPole(values.(m+3))
        (n := retractIfCan(3*q)@Union(Z, "failed")) case Z =>
          one?(s := posrem(n::Z, 6)) => valueOrPole(values.(m+4))
          s = 2 => valueOrPole(values.(m+5))
          s = 4 => valueOrPole(values.(m+6))
          valueOrPole(values.(m+7))
        (n := retractIfCan(4*q)@Union(Z, "failed")) case Z =>
          one?(s := posrem(n::Z, 8)) => valueOrPole(values.(m+8))
          s = 3 => valueOrPole(values.(m+9))
          s = 5 => valueOrPole(values.(m+10))
          valueOrPole(values.(m+11))
        (n := retractIfCan(6*q)@Union(Z, "failed")) case Z =>
          one?(s := posrem(n::Z, 12)) => valueOrPole(values.(m+12))
          s = 5 => valueOrPole(values.(m+13))
          s = 7 => valueOrPole(values.(m+14))
          valueOrPole(values.(m+15))
        "failed"

    else specialTrigs(x, values) == "failed"

    isin x ==
      zero? x => 0
      y := dropfun x
      is?(x, opasin) => y
      is?(x, opacos) => sqrt(1 - y**2)
      is?(x, opatan) => y / sqrt(1 + y**2)
      is?(x, opacot) => inv sqrt(1 + y**2)
      is?(x, opasec) => sqrt(y**2 - 1) / y
      is?(x, opacsc) => inv y
      h  := inv(2::F)
      s2 := h * iisqrt2()
      s3 := h * iisqrt3()
      u  := specialTrigs(x, [[0,false], [0,false], [1,false], [-1,false],
                         [s3,false], [s3,false], [-s3,false], [-s3,false],
                          [s2,false], [s2,false], [-s2,false], [-s2,false],
                           [h,false], [h,false], [-h,false], [-h,false]])
      u case F => u :: F
      kernel(opsin, x)

    icos x ==
      zero? x => 1
      y := dropfun x
      is?(x, opasin) => sqrt(1 - y**2)
      is?(x, opacos) => y
      is?(x, opatan) => inv sqrt(1 + y**2)
      is?(x, opacot) => y / sqrt(1 + y**2)
      is?(x, opasec) => inv y
      is?(x, opacsc) => sqrt(y**2 - 1) / y
      h  := inv(2::F)
      s2 := h * iisqrt2()
      s3 := h * iisqrt3()
      u  := specialTrigs(x, [[1,false],[-1,false], [0,false], [0,false],
                             [h,false],[-h,false],[-h,false],[h,false],
                              [s2,false],[-s2,false],[-s2,false],[s2,false],
                               [s3,false], [-s3,false],[-s3,false],[s3,false]])
      u case F => u :: F
      kernel(opcos, x)

    itan x ==
      zero? x => 0
      y := dropfun x
      is?(x, opasin) => y / sqrt(1 - y**2)
      is?(x, opacos) => sqrt(1 - y**2) / y
      is?(x, opatan) => y
      is?(x, opacot) => inv y
      is?(x, opasec) => sqrt(y**2 - 1)
      is?(x, opacsc) => inv sqrt(y**2 - 1)
      s33 := (s3 := iisqrt3()) / (3::F)
      u := specialTrigs(x, [[0,false], [0,false], [0,true], [0,true],
                      [s3,false], [-s3,false], [s3,false], [-s3,false],
                       [1,false], [-1,false], [1,false], [-1,false],
                        [s33,false], [-s33, false], [s33,false], [-s33, false]])
      u case F => u :: F
      kernel(optan, x)

    icot x ==
      zero? x => INV
      y := dropfun x
      is?(x, opasin) => sqrt(1 - y**2) / y
      is?(x, opacos) => y / sqrt(1 - y**2)
      is?(x, opatan) => inv y
      is?(x, opacot) => y
      is?(x, opasec) => inv sqrt(y**2 - 1)
      is?(x, opacsc) => sqrt(y**2 - 1)
      s33 := (s3 := iisqrt3()) / (3::F)
      u := specialTrigs(x, [[0,true], [0,true], [0,false], [0,false],
                         [s33,false], [-s33,false], [s33,false], [-s33,false],
                          [1,false], [-1,false], [1,false], [-1,false],
                           [s3,false], [-s3, false], [s3,false], [-s3, false]])
      u case F => u :: F
      kernel(opcot, x)

    isec x ==
      zero? x => 1
      y := dropfun x
      is?(x, opasin) => inv sqrt(1 - y**2)
      is?(x, opacos) => inv y
      is?(x, opatan) => sqrt(1 + y**2)
      is?(x, opacot) => sqrt(1 + y**2) / y
      is?(x, opasec) => y
      is?(x, opacsc) => y / sqrt(y**2 - 1)
      s2 := iisqrt2()
      s3 := 2 * iisqrt3() / (3::F)
      h  := 2::F
      u  := specialTrigs(x, [[1,false],[-1,false],[0,true],[0,true],
                           [h,false], [-h,false], [-h,false], [h,false],
                            [s2,false], [-s2,false], [-s2,false], [s2,false],
                             [s3,false], [-s3,false], [-s3,false], [s3,false]])
      u case F => u :: F
      kernel(opsec, x)

    icsc x ==
      zero? x => INV
      y := dropfun x
      is?(x, opasin) => inv y
      is?(x, opacos) => inv sqrt(1 - y**2)
      is?(x, opatan) => sqrt(1 + y**2) / y
      is?(x, opacot) => sqrt(1 + y**2)
      is?(x, opasec) => y / sqrt(y**2 - 1)
      is?(x, opacsc) => y
      s2 := iisqrt2()
      s3 := 2 * iisqrt3() / (3::F)
      h  := 2::F
      u  := specialTrigs(x, [[0,true], [0,true], [1,false], [-1,false],
                            [s3,false], [s3,false], [-s3,false], [-s3,false],
                              [s2,false], [s2,false], [-s2,false], [-s2,false],
                                 [h,false], [h,false], [-h,false], [-h,false]])
      u case F => u :: F
      kernel(opcsc, x)

    iasin x ==
      zero? x => 0
      one? x =>   pi() / (2::F)
      x = -1 => - pi() / (2::F)
      y := dropfun x
      is?(x, opsin) => y
      is?(x, opcos) => pi() / (2::F) - y
      kernel(opasin, x)

    iacos x ==
      zero? x => pi() / (2::F)
      one? x => 0
      x = -1 => pi()
      y := dropfun x
      is?(x, opsin) => pi() / (2::F) - y
      is?(x, opcos) => y
      kernel(opacos, x)

    iatan x ==
      zero? x => 0
      one? x =>   pi() / (4::F)
      x = -1 => - pi() / (4::F)
      x = (r3:=iisqrt3()) => pi() / (3::F)
      one?(x*r3)          => pi() / (6::F)
      y := dropfun x
      is?(x, optan) => y
      is?(x, opcot) => pi() / (2::F) - y
      kernel(opatan, x)

    iacot x ==
      zero? x =>   pi() / (2::F)
      one? x  =>   pi() / (4::F)
      x = -1  =>   3 * pi() / (4::F)
      x = (r3:=iisqrt3())  =>  pi() / (6::F)
      x = -r3              =>  5 * pi() / (6::F)
      one?(xx:=x*r3)       =>  pi() / (3::F)
      xx = -1           =>     2* pi() / (3::F)
      y := dropfun x
      is?(x, optan) => pi() / (2::F) - y
      is?(x, opcot) => y
      kernel(opacot, x)

    iasec x ==
      zero? x => INV
      one? x => 0
      x = -1 => pi()
      y := dropfun x
      is?(x, opsec) => y
      is?(x, opcsc) => pi() / (2::F) - y
      kernel(opasec, x)

    iacsc x ==
      zero? x => INV
      one? x =>   pi() / (2::F)
      x = -1 => - pi() / (2::F)
      y := dropfun x
      is?(x, opsec) => pi() / (2::F) - y
      is?(x, opcsc) => y
      kernel(opacsc, x)

    isinh x ==
      zero? x => 0
      y := dropfun x
      is?(x, opasinh) => y
      is?(x, opacosh) => sqrt(y**2 - 1)
      is?(x, opatanh) => y / sqrt(1 - y**2)
      is?(x, opacoth) => - inv sqrt(y**2 - 1)
      is?(x, opasech) => sqrt(1 - y**2) / y
      is?(x, opacsch) => inv y
      kernel(opsinh, x)

    icosh x ==
      zero? x => 1
      y := dropfun x
      is?(x, opasinh) => sqrt(y**2 + 1)
      is?(x, opacosh) => y
      is?(x, opatanh) => inv sqrt(1 - y**2)
      is?(x, opacoth) => y / sqrt(y**2 - 1)
      is?(x, opasech) => inv y
      is?(x, opacsch) => sqrt(y**2 + 1) / y
      kernel(opcosh, x)

    itanh x ==
      zero? x => 0
      y := dropfun x
      is?(x, opasinh) => y / sqrt(y**2 + 1)
      is?(x, opacosh) => sqrt(y**2 - 1) / y
      is?(x, opatanh) => y
      is?(x, opacoth) => inv y
      is?(x, opasech) => sqrt(1 - y**2)
      is?(x, opacsch) => inv sqrt(y**2 + 1)
      kernel(optanh, x)

    icoth x ==
      zero? x => INV
      y := dropfun x
      is?(x, opasinh) => sqrt(y**2 + 1) / y
      is?(x, opacosh) => y / sqrt(y**2 - 1)
      is?(x, opatanh) => inv y
      is?(x, opacoth) => y
      is?(x, opasech) => inv sqrt(1 - y**2)
      is?(x, opacsch) => sqrt(y**2 + 1)
      kernel(opcoth, x)

    isech x ==
      zero? x => 1
      y := dropfun x
      is?(x, opasinh) => inv sqrt(y**2 + 1)
      is?(x, opacosh) => inv y
      is?(x, opatanh) => sqrt(1 - y**2)
      is?(x, opacoth) => sqrt(y**2 - 1) / y
      is?(x, opasech) => y
      is?(x, opacsch) => y / sqrt(y**2 + 1)
      kernel(opsech, x)

    icsch x ==
      zero? x => INV
      y := dropfun x
      is?(x, opasinh) => inv y
      is?(x, opacosh) => inv sqrt(y**2 - 1)
      is?(x, opatanh) => sqrt(1 - y**2) / y
      is?(x, opacoth) => - sqrt(y**2 - 1)
      is?(x, opasech) => y / sqrt(1 - y**2)
      is?(x, opacsch) => y
      kernel(opcsch, x)

    iasinh x ==
      is?(x, opsinh) => first argument kernel x
      kernel(opasinh, x)

    iacosh x ==
      is?(x, opcosh) => first argument kernel x
      kernel(opacosh, x)

    iatanh x ==
      is?(x, optanh) => first argument kernel x
      kernel(opatanh, x)

    iacoth x ==
      is?(x, opcoth) => first argument kernel x
      kernel(opacoth, x)

    iasech x ==
      is?(x, opsech) => first argument kernel x
      kernel(opasech, x)

    iacsch x ==
      is?(x, opcsch) => first argument kernel x
      kernel(opacsch, x)

    iexp x ==
      zero? x => 1
      is?(x, oplog) => first argument kernel x
      before?(x,0) and empty? variables x => inv iexp(-x)
      h  := inv(2::F)
      i  := iisqrt1()
      s2 := h * iisqrt2()
      s3 := h * iisqrt3()
      u  := specialTrigs(x / i, [[1,false],[-1,false], [i,false], [-i,false],
            [h + i * s3,false], [-h + i * s3, false], [-h - i * s3, false],
             [h - i * s3, false], [s2 + i * s2, false], [-s2 + i * s2, false],
              [-s2 - i * s2, false], [s2 - i * s2, false], [s3 + i * h, false],
               [-s3 + i * h, false], [-s3 - i * h, false], [s3 - i * h, false]])
      u case F => u :: F
      kernel(opexp, x)

-- THIS DETERMINES WHEN TO PERFORM THE log exp f -> f SIMPLIFICATION
-- CURRENT BEHAVIOR:
--     IF R IS COMPLEX(S) THEN ONLY ELEMENTS WHICH ARE RETRACTABLE TO R
--     AND EQUAL TO THEIR CONJUGATES ARE DEEMED REAL (OVERRESTRICTIVE FOR NOW)
--     OTHERWISE (e.g. R = INT OR FRAC INT), ALL THE ELEMENTS ARE DEEMED REAL

    if (R has imaginary:() -> R) and (R has conjugate: R -> R) then
         localReal? x ==
            (u := retractIfCan(x)@Union(R, "failed")) case R
               and (u::R) = conjugate(u::R)

    else localReal? x == true

    iiilog x ==
      zero? x => INV
      one? x => 0
      (u := isExpt(x, opexp)) case Record(var:K, exponent:Integer) =>
           rec := u::Record(var:K, exponent:Integer)
           arg := first argument(rec.var);
           localReal? arg => rec.exponent * first argument(rec.var);
           ilog x
      ilog x

    ilog x ==
      ((num1 := one?(num := numer x)) or num = -1) and not one?(den := denom x)
        and empty? variables x => - kernel(oplog, (num1 => den; -den)::F)
      kernel(oplog, x)

    if R has ElementaryFunctionCategory then
      iilog x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iiilog x
        log(r::R)::F

      iiexp x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iexp x
        exp(r::R)::F

    else
      iilog x == iiilog x
      iiexp x == iexp x

    if R has TrigonometricFunctionCategory then
      iisin x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => isin x
        sin(r::R)::F

      iicos x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => icos x
        cos(r::R)::F

      iitan x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => itan x
        tan(r::R)::F

      iicot x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => icot x
        cot(r::R)::F

      iisec x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => isec x
        sec(r::R)::F

      iicsc x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => icsc x
        csc(r::R)::F

    else
      iisin x == isin x
      iicos x == icos x
      iitan x == itan x
      iicot x == icot x
      iisec x == isec x
      iicsc x == icsc x

    if R has ArcTrigonometricFunctionCategory then
      iiasin x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iasin x
        asin(r::R)::F

      iiacos x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iacos x
        acos(r::R)::F

      iiatan x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iatan x
        atan(r::R)::F

      iiacot x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iacot x
        acot(r::R)::F

      iiasec x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iasec x
        asec(r::R)::F

      iiacsc x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iacsc x
        acsc(r::R)::F

    else
      iiasin x == iasin x
      iiacos x == iacos x
      iiatan x == iatan x
      iiacot x == iacot x
      iiasec x == iasec x
      iiacsc x == iacsc x

    if R has HyperbolicFunctionCategory then
      iisinh x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => isinh x
        sinh(r::R)::F

      iicosh x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => icosh x
        cosh(r::R)::F

      iitanh x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => itanh x
        tanh(r::R)::F

      iicoth x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => icoth x
        coth(r::R)::F

      iisech x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => isech x
        sech(r::R)::F

      iicsch x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => icsch x
        csch(r::R)::F

    else
      iisinh x == isinh x
      iicosh x == icosh x
      iitanh x == itanh x
      iicoth x == icoth x
      iisech x == isech x
      iicsch x == icsch x

    if R has ArcHyperbolicFunctionCategory then
      iiasinh x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iasinh x
        asinh(r::R)::F

      iiacosh x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iacosh x
        acosh(r::R)::F

      iiatanh x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iatanh x
        atanh(r::R)::F

      iiacoth x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iacoth x
        acoth(r::R)::F

      iiasech x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iasech x
        asech(r::R)::F

      iiacsch x ==
        (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iacsch x
        acsch(r::R)::F

    else
      iiasinh x == iasinh x
      iiacosh x == iacosh x
      iiatanh x == iatanh x
      iiacoth x == iacoth x
      iiasech x == iasech x
      iiacsch x == iacsch x

    evaluate(oppi, ipi)$BasicOperatorFunctions1(F)
    evaluate(oplog, iilog)
    evaluate(opexp, iiexp)
    evaluate(opsin, iisin)
    evaluate(opcos, iicos)
    evaluate(optan, iitan)
    evaluate(opcot, iicot)
    evaluate(opsec, iisec)
    evaluate(opcsc, iicsc)
    evaluate(opasin, iiasin)
    evaluate(opacos, iiacos)
    evaluate(opatan, iiatan)
    evaluate(opacot, iiacot)
    evaluate(opasec, iiasec)
    evaluate(opacsc, iiacsc)
    evaluate(opsinh, iisinh)
    evaluate(opcosh, iicosh)
    evaluate(optanh, iitanh)
    evaluate(opcoth, iicoth)
    evaluate(opsech, iisech)
    evaluate(opcsch, iicsch)
    evaluate(opasinh, iiasinh)
    evaluate(opacosh, iiacosh)
    evaluate(opatanh, iiatanh)
    evaluate(opacoth, iiacoth)
    evaluate(opasech, iiasech)
    evaluate(opacsch, iiacsch)
    derivative(opexp, exp)
    derivative(oplog, inv)
    derivative(opsin, cos)
    derivative(opcos, - sin #1)
    derivative(optan, 1 + tan(#1)**2)
    derivative(opcot, - 1 - cot(#1)**2)
    derivative(opsec, tan(#1) * sec(#1))
    derivative(opcsc, - cot(#1) * csc(#1))
    derivative(opasin, inv sqrt(1 - #1**2))
    derivative(opacos, - inv sqrt(1 - #1**2))
    derivative(opatan, inv(1 + #1**2))
    derivative(opacot, - inv(1 + #1**2))
    derivative(opasec, inv(#1 * sqrt(#1**2 - 1)))
    derivative(opacsc, - inv(#1 * sqrt(#1**2 - 1)))
    derivative(opsinh, cosh)
    derivative(opcosh, sinh)
    derivative(optanh, 1 - tanh(#1)**2)
    derivative(opcoth, 1 - coth(#1)**2)
    derivative(opsech, - tanh(#1) * sech(#1))
    derivative(opcsch, - coth(#1) * csch(#1))
    derivative(opasinh, inv sqrt(1 + #1**2))
    derivative(opacosh, inv sqrt(#1**2 - 1))
    derivative(opatanh, inv(1 - #1**2))
    derivative(opacoth, inv(1 - #1**2))
    derivative(opasech, - inv(#1 * sqrt(1 - #1**2)))
    derivative(opacsch, - inv(#1 * sqrt(1 + #1**2)))

