new defint.lisp bug - was: new defint.lisp and radexpand:false?



On Jan. 11, 2012, Richard Fateman wrote:
----------------------------------
>> (%i3) integrate(exp(x^5),x,1,2);
>> (%o3) (gamma_incomplete(1/5,-32)-gamma_incomplete(1/5,-1))/5
>
>for what it is worth, Mathematica 7 returns
>
>1/5 (-1)^(4/5) (Gamma[1/5, -32] - Gamma[1/5, -1])}
>
>
>Notice the factor of (-1)^(4/5).
-------------------------------------------
This is a very important observation, and explains
a difference between Maxima and Mathematica results
for these type of integrals.  Maxima evaluates the
factor (-1)^(4/5) numerically to 1.0, which differs
from Mathematica's evaluation (see below).

Here I consider a different (but related via a simple
change of variables)  integral, and compare in detail
the results from Maxima and Mathematica (using the
Wolfram alpha site).

-----------------------------------------------------
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)


(%i1) integrate(exp(y)/y^(4/5),y,1,2);
(%o1) gamma_incomplete(1/5,-2)-gamma_incomplete(1/5,-1)

/*  MMa NIntegrate[Exp[y]/y^(4/5),{y,1,2}] yields
  
    (-1)^(4/5)*(Gamma[1/5,-2] - Gamma[1/5,-1])    
    
    does not agree with Maxima, since in Mathematica
    (-1)^(4/5) float value is not 1.0, whereas
    is does have value 1.0 in Maxima.  */


(%i2) expand(float(%o1));
(%o2) -1.951024982232477*%i-2.685355511932373

/* this Maxima float value agrees with the first 14 digits of
     
     Mma's   N[Gamma[1/5,-2] - Gamma[1/5,-1]] which gives 
  
    -2.685355511932382... -1.951024982232483... i
    
    which translates into
    
    -2.685355511932382 - 1.951024982232483*%i     */
    
 /* the approximate (real) value of the integral
    is given by quad_qags  */     
    

(%i3) quad_qags(exp(y)/y^(4/5),y,1,2);
(%o3) [3.319281956502171,3.6851432533315483E-14,21,0]

   /* Maxima gives value 1.0 to (-1)^(4/5)  */

(%i4) expand(float(-1)^(4/5));
(%o4) 1.0


  /* whereas Mathematica gives
  
    N[(-1)^(4/5)]  --->  
    
      -0.8090169943749474... + 0.5877852522924731... i

      which translates into
      
      -0.8090169943749474 + 0.5877852522924731*%i     
      
 If we multiply this factor by %o2 and expand we get the
   correct numerical answer to within floating point roundoff
   errors:    */
   
(%i5) expand((-0.8090169943749474 + 0.5877852522924731*%i)*%o2);

(%o5) 4.4408920985006262E-16*%i+3.319281956502161


/* the above was with 5.25.1 defint.lisp.
   The updated version makes no difference here  */
   
   
(%i6) load("defint-new.lisp");
(%o6) "c:/work2/defint-new.lisp"

(%i7) integrate(exp(y)/y^(4/5),y,1,2);
(%o7) gamma_incomplete(1/5,-2)-gamma_incomplete(1/5,-1)

(%i8) expand(float(%));
(%o8) -1.951024982232477*%i-2.685355511932373

--------------------------------------
Ted