7.87 bug in integration in Maple 6 and fix (16.5.01)

7.87.1 Josef Hekrdla

I am really horrified! If you integrate elementary function 1/(9x^2+6x+2), you obtain the wrong result arctan(1+3x) instead of the correct 1/3*arctan(1+3x)! It is only for Maple6. Maple5 returns the correct result.

It is corrected with Maple 7. (U. Klein)

7.87.2 Carl DeVore(18.5.01)

This bug results from trying to program a great variety of simple cases rather than just using a handful of algorithms for the general cases. In this problem, the special case that is handled in procedure `int/ratpoly/ratpoly` is when the denominator is exactly one more than a perfect square of first-degree binomial. In other words,

   > int(c/((a*x+b)^2+1), x)
 

will give the wrong answer. Change that 1 to anything else and you’ll get the right answer. This case is so special – why do they even bother with it? It’s hard to imagine that it is worth the bother to treat it separately.

The general case algorithms need to be there anyway. Does it really save that much time to sort out and handle the simple cases? When there are so many special cases, it makes it that much more likely that a trivial bug will go undetected for a long time. Does Maple have an official programming philosophy?

In this case, the bug is truly trivial – essentially a typo – and easy for the user to correct. The fact that someone did not notice this bug sooner is a sign that it is not worth the trouble to treat this case separately. I’ve traced it to lines 21 and 23 of procedure `int/ratpoly/ratpoly`. I show the fix after the code below.

> restart; 
> showstat(`int/ratpoly/ratpoly`); 
 
`int/ratpoly/ratpoly` := proc(f) 
local ans, den, num, d, n, k, const, q, r, rest, g, const2, A, B, answer; 
   1   if type(f,polynom(anything,_X)) then 
   2     RETURN(`int/polynom`(args)) 
       end if; 
   3   if type(f,`+`) then 
   4     return map(procname,f) 
       end if; 
   5   num := 1; 
   6   const := 1; 
   7   if type(f,`*`) then 
   8     rest, const := selectremove(has,f,_X); 
   9     num, den := numer(rest), denom(rest) 
       else 
  10     const, num, den := 1, 1, denom(f) 
       end if; 
  11   den, n := op(`if`(den::`^`,den,[den, 1])); 
  12   d := degree(den,_X); 
  13   answer := FAIL; 
  14   if d = 1 then 
  15     const := const*num/coeff(den,_X,1); 
  16     if n = 1 then 
  17       answer := const*ln(den) 
         else 
  18       answer := -const/(n-1)/(den^(n-1)) 
         end if 
       elif n = 1 and degree(num,_X) = d-1 and rem(num,diff(den,_X),_X,q) = 0 then 
  19     answer := const*q*ln(den) 
       elif n = 1 and d = 2+2*degree(num,_X) and type(den,polynom(polynom(rational),_X)) 
            and psqrt(den-1) <> _NOSQRT and rem(diff(den-1,_X),num,_X) = 0 then 
  20     r := psqrt(den-1); 
  21     answer := const/content(r,_X)*arctan(r) 
^^^^^^^^^^^^^ The bug is that they divide by the common factor (the "content") of the binomial rather than by its leading coefficient. 
 
       elif n = 1 and d = 2+2*degree(num,_X) and type(den,polynom(polynom(rational),_X)) 
            and psqrt(den+1) <> _NOSQRT and rem(diff(den+1,_X),num,_X) = 0 then 
  22     r := psqrt(den+1); 
  23     answer := const/content(r,_X)*arctanh(r) 
                         ^^^^^^^^^^^^^
 

Same bug here.

       elif d = 2 then 
  24     answer := const*`int/ratpoly/quadratic`(num,den,n-1) 
       elif d = 3 then 
  25     r := `int/ratpoly/cubic`(num,den,n); 
  26     answer := `if`(r = FAIL,FAIL,const*r) 
       elif 0 < ldegree(num,_X) and `int/ratpoly/subs`(num,den,'k') then 
  27     num := collect(num/(_X^(k-1)),_X); 
  28     r := `int/indef1`(subs(_X = _X^(1/k),num/(den^n))); 
  29     answer := 1/k*const*subs(_X = _X^k,r) 
       elif irem(d,2,'q') = 0 and degree(num,_X) = q-1 and 
            `int/ratpoly/arctan`(num,den^(n+1),'q') then 
  30     answer := const*q 
       elif d = 4 then 
  31     r := `int/ratpoly/quartic`(num,den,n); 
  32     answer := `if`(r = FAIL,FAIL,const*r) 
       elif d = 5 then 
  33     r := `int/ratpoly/quintic`(num,den,n-1); 
  34     answer := `if`(r = FAIL,FAIL,const*r) 
       elif d = 6 then 
  35     r := `int/ratpoly/sextic`(num,den,n-1); 
  36     answer := `if`(r = FAIL,FAIL,const*r) 
       end if; 
  37   if answer = FAIL then 
  38     if type(f,'ratpoly(rational,_X)') then 
  39       const*`int/risch/ratpoly`(num/(den^n),_X) 
         else 
  40       g, const2, A, B := `int/ratpoly/horowitz`(num/(den^n),_X); 
  41       if A = 0 then 
  42         const*g 
           else 
  43         const*(g+const2*sum(subs(_X = _R,A/diff(B,_X))*ln(_X-_R),_R = RootOf(B,_X))) 
           end if 
         end if 
       else 
  44     answer 
       end if 
end proc
 

An example of the bug in action:

> int(c/((a*x+b)^2+1), x); 
 
                          c arctan(a x + b)
 

Before you test my fix, make sure to do a restart to clear the erroneous answers from the remember tables.

> restart;
 

Since the "content" command is used only on the erroneous lines, you can easily correct this bug:

> `int/ratpoly/ratpoly`:= subs(content= lcoeff, eval(`int/ratpoly/ratpoly`)):
 

You could also use "diff" instead of "lcoeff". That would be more natural.

Test it:

> int(c/((a*x+b)^2+1), x); 
 
                          c arctan(a x + b) 
                          ----------------- 
                                  a 
 
> int(1/(9*x^2+6*x+2), x); 
 
                         1/3 arctan(1 + 3 x)
 

In release 5, these special cases are not treated. They are merely passed on to `int/ratpoly/quadratic` as in line 24 above. Why were these cases singled out in release 6?

I am not absolutely sure that my fix works in all cases. Be wary when doing integrals with a quadratic denominator.