Re: [Axiom-developer] [Bug?] "error in library; negative log"

[William Sit]

Karl Hegbloom wrote:
> 
> Attached is a session log showing an error that I receive while
> attempting to produce a sequence from an expression in Axiom.  Maxima
> seems to have no trouble with the similar expression, and computing the
> value of the expression by hand, as you can see, seems to work fine
> also.
> 
> Another problem I have is that taking the limit of an expression
> containing (-1)^n always returns "failed", where my TI-89 Titanium
> calculator will give a finite limit.  For example:
> 
>  limit( 2 + (-2/%pi)^n, n=%plusInfinity )  ===> "failed"
> 
> ... but the TI-89t returns 2.

Mathematically, the limit is 2 since (-2/%pi) has absolute value less than 1,
and hence (-2/%pi)^n converges to 0. So TI-89t is correct and Axiom is wrong.

> The TI-89t says that the limit of (-1)^n as n approaches infinity is -1,
> implying that it believes that infinity is an odd number.  That kind of
> makes sense to me, since if you divide infinity in half, you still have
> infinity, and you keep adding 1 to get to infinity, making it odd.  If
> infinity is even then the answer should be 1, and if we can't know if
> infinity is even or odd, then the answer is uncertain or undefined.

No, TI-89t is wrong to say that limit (-1)^n is -1. The limit does not exist
because the sequence (-1)^n oscillates between 1 and -1. There is no number L
(the assumed limit) such that given any epsilon > 0, there is a natural number N
such that |(-1)^n - L| < epsilon for all n > N.

> On the other hand, the TI-89t says that lim ( (-1)^n * (n + 1)/n ) is
> undefined. 

This is correct since (n + 1)/n = 1 + 1/n has a limiting value 1 as n approaches
infinity, but (-1)^n makes that alternating in sign, and hence the limit does
not exist.

> But it already told me that lim (-1)^n = -1, and that lim (n
> + 1)/n = 1.  If the limit of a product is the product of the limits of
> the factors, then lim ( (-1)^n * (n + 1)/n ) should be -1, right?

See above. The limit of a product is the product of the limits of the factors,
provided each  factor has its own limit. Here (-1)^n does not have a limit. 
Note that if you consider (-1)^n * (1/n), then the limit would exist and is 0,
and the converse does not hold (it is not necessary that each factor has a limit
for the product to have a limit).

> Karl Hegbloom <address@hidden>
> 
[coke and GCL ads snipped]
[AXIOM loading messages snipped, ditto below] 

To simplify the outputs and explanation, I have replaced your function f by
removing the initial addition of 2.0, which is not relevant to the discussion. I
have also replaced the floating point (-2.0/%pi::Float) by just a floating point
number -3.0. I have replaced the loop by a single eval. The follwwing transcript
is extracted from mine.

> (1) -> f == (-3.0)**n
                                                                   Type: Void
> (2) -> eval(f,n=3)

>    Compiling body of rule f to compute value of type Expression Float
>    >> Error detected within library code:
>    negative log

This is because you did not specify in the definition of f that n is an 
integer. 
Note that when you do not specify the type of arguments for a function, the
Interpreter may use any type specification as it sees fit at the time the
function is called. In this case, the Interpreter needs to figure out which eval
is to be used and it compiles f with target in Expression Float because n has
not been given a value and is treated as a symbol. The function f is technically
a constant function with value in Expression Float. 

> (3) -> f

               n
   (3)  (- 3.0)
                                                       Type: Expression Float

Note that you DID say f takes NO input argument. It would make a difference if
you had used:

> (2) ->  g n == (-3.0)**n
                                                                   Type: Void

> (3) ->  g 3
>   Compiling function g with type PositiveInteger -> Float
   (3)  - 27.0
                                                                  Type: Float
To see what happened in (2),

> )set mess bot on
> (4) -> f == (-3.0)**n
>   Compiled code for f has been cleared.
>   1 old definition(s) deleted for function or rule f
                                                                   Type: Void
> (5) -> eval(f,n=3)

[snipped]
> Function Selection for **
>     Arguments: (FLOAT,VARIABLE n)
>     Default target type: Expression Float

>[1]  signature:   (EXPR FLOAT,EXPR FLOAT) -> EXPR FLOAT
>     implemented: slot $$$ from EXPR FLOAT

>  Compiling body of rule f to compute value of type Expression Float
[snipped]
> Function Selection for eval
>      Arguments: (EXPR FLOAT,EQ POLY INT)
>      Default target type: Expression Float

> [1]  signature:   (EXPR FLOAT,EQ EXPR FLOAT) -> EXPR FLOAT
>      implemented: slot $$(Equation $) from EXPR FLOAT
[snipped]
>   >> Error detected within library code:
>   negative log

Now note that ** is taken from Expression Float, not Float and note the
signature of **. Axiom implements the exponentiation function in EXPR by
distinguishing whether the exponent is integer or otherwise.  The "otherwise"
cases are 
rather complicated, but suffices to say that a meaningful way to compute x**n
for a SYMBOL n with x:Expression Float in general is to use

   x**n == exp(n*log(x))

Now combine this with the fact you asked Axiom to FIRST compute f, THEN evaluate
by substituting n = 3, it should be clear why the computation produces the error
message above (in the special case x retracts to Float,
log is not defined for negative x). See float.spad:

   log x ==
      negative? x => error "negative log"
      zero? x => error "log 0 generated"

--- skip if you like
If anyone is inclined to get to the bottom of this, here is a summary of the
source:

  The implementation of ** in Expression (see expr.spad) in the case when a
symbol n is coerced to EXPR FLOAT is:

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

where CF is CombinatorialFunction(Float, Expression Float) package. We find in
combfunc.spad:

    x ** y               == 
      -- Do some basic simplifications
      is?(x,POWER) =>
        args : List F := argument first kernels x
        not(#args = 2) => error "Too many arguments to **"
        number?(first args) and number?(y) =>
          oppow [first(args)**y, second args]
        oppow [first args, (second args)* y]
      -- Generic case
      exp : Union(Record(val:F,exponent:Z),"failed") := isPower x
      exp case Record(val:F,exponent:Z) =>
        expr := exp::Record(val:F,exponent:Z)
        oppow [expr.val, (expr.exponent)*y]
      oppow [x, y]
and
    oppow   := operator(POWER::Symbol)$CommonOperators

>From CommonOperators:

    oppow   := operator(POWER, 2)

It is interesting to note what f is compiled into:

(6) -> f::InputForm

   (6)
   (/

     (+

       (+  (float 0 0 2)

         (*  (float 1 0 2)

           (**  (/ (float - 3 0 2) (float 1 0 2))

             (/  (+ (+ (float 0 0 2) (* (float 1 0 2) n)) (float 0 0 2))
              (float 1 0 2))
             )
           )
         )

      (float 0 0 2))

    (float 1 0 2))
                                                              Type: InputForm

Here (float a b c) means floating numbers using base c, with exponent b and
mantissa a, which is mathematically a*c^b.  So (float 0 0 2) is simply 0, (float
1 0 2) is 1, and (float - 3 0 2) is (-3.0).
Thus, one can see that the compiled Lisp code has a lot of cruft (adding 0 and
multiplying by 1 all over the place; however, this may be Lisp way of coercion
and I don't know how Lisp interacts with Axiom), but the compiled form for the
above f is equivalent to
           (** (float - 3 0 2) n)
(not surprising). Here n is a symbol and the expression belongs to the domain
Expression Float, and the only signature the Interpreter can find is what I
showed above.
---- 

Continuing with your output transcript:

> frame0 (5) -> -2.0/%pi::Float
> 
>    (5)  - 0.6366197723 6758134307 55351
>                                                                   Type: Float
> frame0 (6) -> % ** 2
> 
>    (6)  0.4052847345 6935108577 55179

                                                                    Type: Float
Here you actually supplied an integer 2 as exponent. This is implemented by
repeated multiplication (or something more clever such as Russian
multiplication). In Float, Axiom also implements ** by cases. See float.spad.

To do what you intend:
(7) -> g n == 2.0+(-2/%pi::Float)**n
   Compiled code for g has been cleared.
   1 old definition(s) deleted for function or rule g
                                                                   Type: Void
(8) -> [[i, g i] for i in 1..20]
   Compiling function g with type PositiveInteger -> Float

   (7)
   [[1.0,1.3633802276 3241865692], [2.0,2.4052847345 6935108578],
    [3.0,1.7419877245 3440408652], [4.0,2.1642557160 7494936303],
    [5.0,1.8954315634 2229166488], [6.0,2.0665703342 9093454695],
    [7.0,1.9576200089 3727145417], [8.0,2.0269799402 6329437858],
    [9.0,1.9828240365 7109059006], [10.0,2.0109345579 2830621047],
    [11.0,1.9930388442 2074156846], [12.0,2.0044316094 0760677601],
    [13.0,1.9971787498 2770734225], [14.0,2.0017960636 4247695142],
    [15.0,1.9988565903 7276863419], [16.0,2.0007279171 7661093314],
    [17.0,1.9995365935 3272349524], [18.0,2.0002950137 1971123349],
    [19.0,1.9998121884 3291212109], [20.0,2.0001195645 5708748421]]
                                                        Type: List List Float
   
William
 
William Sit
Department of Mathematics....Email: address@hidden
City College of New York................Tel: 212-650-5179
Convent Ave at West 138th Street........Fax: 212-862-0004
New York, NY 10031..Axiom, A Scientific Computation Sytem
USA............... http://www.nongnu.org/axiom/index.html