Subject: Numerical evaluation of the power function
From: Dieter Kaiser
Date: Thu, 28 Jan 2010 20:34:47 +0100
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.
...
;; Check for numerical evaluation of the sqrt.
((and (alike1 pot '((rat) 1 2))
(or (setq res (flonum-eval '%sqrt gr))
(and (not (member 'simp (car x) :test #'eq))
(setq res (big-float-eval '%sqrt gr)))))
(return res))
;; Numerical evaluation
((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)))))
...
> 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.
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
implementation of the numerical evaluation with rectform and float will
have similar problems:
(%i6) float(rectform((1.0+1.0*%i)^2));
(%o6) 2.0 %i
(%i7) float(rectform((1.0+1.0*%i)^2.0));
(%o7) 2.0 %i + 1.224500708041117E-16
> > /* 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.
> > /* 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.
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.
Dieter Kaiser