Numerical evaluation of the power function

On 1/28/10 2:34 PM, Dieter Kaiser wrote:
> Am Montag, den 25.01.2010, 09:25 -0500 schrieb Raymond Toy: 
>> On 1/23/10 8:48 AM, Dieter Kaiser wrote:
>>> Am Mittwoch, den 20.01.2010, 20:11 -0600 schrieb Barton Willis:
>>>   
>>>> Instead, would it be possible to use Raymond Toy's bigfloat package to do
>>>> floating point evaluation of ^?
>>>>     
>>> Yes, it is possible and perhaps a much better implementation, because it
>>> is faster. I have tried the following code:
>>>
>>> ((and (or (complex-float-numerical-eval-p gr pot)
>>>           (complex-bigfloat-numerical-eval-p gr pot))
>>>       (or (not (ratnump gr)) (not (ratnum pot))))
>>>  (let ((z (bigfloat:to gr))
>>>        (a (bigfloat:to pot)))
>>>    (return (to (bigfloat:expt z a)))))
>>>   
>> Can you tell me where exactly I should put the above code?
> 
> 
> I have put the new code just behind the numerical evaluation for the
> sqrt function. I have not cut out this code, because we have special
> routines to handle this case.

Thanks.  Now I can do tests you give.

>>> (%i1) float((1+%i)^2);
>>> (%o1) 2.0 %i + 1.224606353822378e-16

Curiously, clisp returns 2.0*%i.
> 
> 
> I think the only way to get a pure imaginary and most simple result is
> to have special code to handle exponents which represent an integer. An

What do you have in mind?  Are you saying x^y where y is equal to a
integer will compute the result by repeated squaring/multiplications?
If so, this also has some numerical issues, unless y is kept small.


> 
>>> /* This result seems to have only single float precision */
>>>
>>> (%i4) (1+%i)^2.0;
>>> (%o4) 2.000000003809306 %i - 8.742278017023905e-8
> 
>> I don't quite understand this one.  I think we should be calling
>> bigfloat:expt with 1+%i and 2.0.  The 1+%i becomes #c(1 1), and we
>> compute (cl:expt  #c(1 1) 2d0).  This ought to produce a double-float
>> result.  What lisp were you using for this?
> 
> 
> The result from above I have got with SBCL 1.0.29. But I can not
> reproduce it with CLISP 2.44. With CLISP I get 2.0*%i as a result. 

This seems like a bug in sbcl.  What does sbcl say for (cl:expt #c(1 1)
2d0)?

> OK, I see the problem. Usually, we convert all numbers to a float number
> before we do the numerical evaluation. Because the bigfloat package can
> handle all types of Maxima numbers this conversion is not done and the
> flag $numer no longer works as expected. Therefore, I think the expected
> conversion of numbers has to be build into the bigfloat package in
> addition.

I think that would be against the design principle of the bigfloat
package.  The intent is that if you have some Lisp code that works for
CL numbers (integer, single float, double float, whatever), then you can
just (essentially) compile the code in the bigfloat package to use
bigfloat operations.  Then the routine would continue to work as before,
but also be extended to work with bigfloats automatically.  So the
conversion that you want needs to be done outside of bigfloat.  Unless
we want to change the design principle for bigfloat.

Ray