Subject: Numerical evaluation of the power function
From: Raymond Toy
Date: Sat, 23 Jan 2010 10:23:02 -0500
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)))))
>
> 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
>
bigfloat:expt is not very smart and uses exp and log to compute the
result. Hence, there's a fair amount of roundoff. This needs more work
and is the function that I've spent the least amount of time on.
>
> /* This result seems to have only single float precision */
>
> (%i4) (1+%i)^2.0;
> (%o4) 2.000000003809306 %i - 8.742278017023905e-8
>
Oops. I'll fix that asap.
Ray