Numerical evaluation of the power function

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?

In the meantime, I just created the simple function:

(defun $bfexp (z a)
  (let ((zz (bigfloat:to z))
    (aa (bigfloat::to a)))
    (to (bigfloat::expt zz aa))))

for testing.
> But I have observed some inconsistencies which might be due to the
> implementation of the expt function in the numerical package:
>
> /* Better would be 2.0*%i */
>
> (%i1) float((1+%i)^2);
> (%o1) 2.0 %i + 1.224606353822378e-16
>   
IIRC, float((1+%i)^2) floats the parts and the numerically evaluates
them.  If so, bigfloat:expt calls cl:expt with the arguments #C(1.0,
1.0) and 2.0, I think.  So that's the source of error here.
>
> /* 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?
> /* Compare the result form above with this one,
>    both results should be the same. */
>
> (%i5) (1.0+%i)^2.0;
> (%o5) 2.0 %i + 1.224606353822377e-16
>
> /* I think the result should be 2.0*%i 
>    and it should be identical to the first example */
>
> (%i18) (1+%i)^2,numer;
> (%o18) 2 %i
>
> There might be some more small inconsistencies.
>   
Not exactly sure yet on how to remove these.

Ray