Subject: Numerical evaluation of the power function
From: Dieter Kaiser
Date: Sat, 23 Jan 2010 14:48:21 +0100
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
(%i2) (1+%i)^2,bfloat;
(%o2) 2.0b0 %i - 3.216245299353273b-16
(%i3) (1.0+%i)^2;
(%o3) 2.0 %i + 1.224606353822378e-16
/* This result seems to have only single float precision */
(%i4) (1+%i)^2.0;
(%o4) 2.000000003809306 %i - 8.742278017023905e-8
/* 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.
Perhaps, we can remove this inconsistencies and implement the numerical
evaluation of the power function with the bigfloat package.
Dieter Kaiser