new defint.lisp bug - was: new defint.lisp and radexpand:false?
Subject: new defint.lisp bug - was: new defint.lisp and radexpand:false?
From: Edwin Woollett
Date: Wed, 11 Jan 2012 10:36:21 -0800
On Jan. 11, 2012, I wrote:
--------------------------
>Here is an example which starts with the current
>5.25.1 defint.lisp, resulting in a wrong answer,
-------------------------
To demonstrate that the wrong answer is not due to
insufficient precision in the numerical work, I use
mfloat(expr, ndigits), previously discussed on the
list, which goes for an answer accurate to ndigits
using bfloat methods.
--------------------------
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
(%i1) mfloat1(ex,ndigits) :=
block([fpprec:16,old,new,kk,small,goodans:false],
small:bfloat(10^(-(ndigits+1))),
old : expand(bfloat(ex)),
for kk:2 thru 10 do
(new:realpart(block([fpprec:kk*10],expand(bfloat(ex)))),
if abs(new-old) < small then
(goodans:true,return()),
old:new),
if not goodans then false
else if not numberp(new) then false
else float(new))$
(%i2) mfloat(expr,mdigits):=
block([er,ei,err],
er : realpart(expr),
err : mfloat1(er,mdigits),
if not complex_number_p(err) then
(print(" realpart not evaluated"),
return(false)),
if expand(expr - er) = 0 then err
else
(ei : imagpart(expr),
ei : mfloat1(ei,mdigits),
if not complex_number_p(ei) then
(print(" imagpart not evaluated"),
return(false)),
err + %i*ei))$
(%i3) quad_qags(exp(x^5),x,1,2);
(%o3) [1.0132394896940184E+12,0.66155076101722,147,0]
(%i4) i1: integrate(exp(x^5),x,1,2);
(%o4) (gamma_incomplete(1/5,-32)-gamma_incomplete(1/5,-1))/5
(%i5) mfloat(i1,15);
(%o5) -5.9556722908249524E+11*%i-8.1972796653426025E+11
----------------------------
Ted