#### 7.87 bug in integration in Maple 6 and ﬁx (16.5.01)

##### 7.87.1 Josef Hekrdla

I am really horriﬁed! 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 ﬁrst-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 oﬃcial 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 ﬁx 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);
14   if d = 1 then
15     const := const*num/coeff(den,_X,1);
16     if n = 1 then
else
end if
elif n = 1 and degree(num,_X) = d-1 and rem(num,diff(den,_X),_X,q) = 0 then
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);
^^^^^^^^^^^^^ 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);
^^^^^^^^^^^^^



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
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
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 ﬁx, 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 "diﬀ" instead of "lcoeﬀ". 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 ﬁx works in all cases. Be wary when doing integrals with a quadratic denominator.