Fix for Larry Bates's integration bug



"quotient is not exact" is a message from the polynomial
GCD/ reduction program.  Basically the program determined
that A/B  appears and that GCD(A,B) = G.  It then
needs to divide A/G  and B/G  and proceed with the reduced
fraction.  But what it found was that G did not divide
A  or B.

One possibility is there is a bug in the GCD routine, in
which case using one of the others may help.  If you look
at the setting of GCD, and it is SPMOD,  you might try
doing the same calculation after doing

GCD:SUBRES;

If this works, then the error is in the GCD routine.

It is much less likely that there is a bug in the subresultant
GCD.


I suppose there are other possible ways of provoking this
message other than the GCD.  If it is possible to do a
back-trace and see what functions were active when this
message occurred, that might give a hint.

RJF


TJG


Wolfgang Jenkner wrote:

> Camm Maguire <camm@enhanced.com> writes:
> 
> 
>>(C1) f(a):=integrate(sin(a)/(1-cos(a)^2*sin(t)^2), t,0,2*%Pi);
>>
>>                                      SIN(a)
>>(D1)           f(a) := INTEGRATE(-------------------, t, 0, 2 %PI)
>>                                        2       2
>>                                 1 - COS (a) SIN (t)
>>(C2) f(%Pi/2);
>>
>>(D2)                      2 %PI
>>(C3) f(3.14159/2);
>>
> 
> [...]
> 
> 
>>(D3)                        0
>>
>>The correct answer is 2Pi, and in fact does not depend on the choice
>>of the variable a if 0<a<Pi.
>>
>>Larry Bates
>>
>>
> 
> The origin of the bug is revealed by tracing the function PERIOD (in
> defint), which tests if some expression is periodic with respect to a
> given period and variable. I trace some more functions to tell the
> whole story.
> 
> (C1) trace(?intsc1,?period,?sinint);
> (D1) 			   [INTSC1, PERIOD, SININT]
> (C2) f(a):=integrate(sin(a)/(1-cos(a)^2*sin(t)^2), t,0,2*%pi);
> 				      SIN(a)
> (D2) 	      f(a) := INTEGRATE(-------------------, t, 0, 2 %PI)
> 				       2       2
> 				1 - COS (a) SIN (t)
> (C3) f(1.5);
> 
> RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
> 
> RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
> 					 1
> 1 Enter ?INTSC1 [0, 2 %PI, -----------------------------]
> 						    2
> 			   1 - .0050037516997773 SIN (t)
> 				       1
>  1 Enter ?PERIOD [2 %PI, -----------------------------, t]
> 						  2
> 			 1 - .0050037516997773 SIN (t)
> 
> RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
>  1 Exit  ?PERIOD FALSE
> 1 Exit  ?INTSC1 FALSE
> 
> RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
> 			       1
> 1 Enter ?SININT [-----------------------------, t]
> 					  2
> 		 1 - .0050037516997773 SIN (t)
> 
> RAT replaced -.005003751699777 by -1107//221234 = -.005003751683738
> 			      220127 TAN(t)
> 		221234 ATAN(-----------------)
> 			    SQRT(48699576718)
> 1 Exit  ?SININT ------------------------------
> 		      SQRT(48699576718)
> (D3) 				       0
> (C4) untrace(?intsc1,?period,?sinint);
> (D4) 			   [INTSC1, PERIOD, SININT]
> 
> The bug in PERIOD (which results from comparing expressions which have
> been simplified in different ways) has the consequence that defint
> thinks that INTSC1 can't be used, so SININT is tried to find an
> indefinite integral; the resulting expression, however, is valid only
> in (open) intervals where tan has no pole...
> 
> After applying the patch, I get the result
> 
> (D7) 		  2.235435669225693E-6 SQRT(800452634953) %PI
> (C8) %/%pi,numer;
> (D8) 			       1.999999999999074
> 
> BTW, I replaced your 3.14159/2 by 1.5 in order to get slightly less
> ugly numbers and also since otherwise apparently some numbers become
> too small (in absolute value) and that vile `LOG(0) has been
> generated' error turns up somehow.
> 
> This seems to completely fix the case where a is a float.  Sadly, in
> other cases we run into some unrelated bug.
> 
> (C9) f(3/2);
> quotient is not exact
> 
> I haven't fixed this one.  Note, however, that other methods like
> doing a partial fraction decomposition wrt sin(t) beforehand will
> work.
> 
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> Index: defint.lisp
> ===================================================================
> RCS file: /cvsroot/maxima/maxima/src/defint.lisp,v
> retrieving revision 1.3
> diff -C2 -r1.3 defint.lisp
> *** defint.lisp	11 Jul 2001 05:06:22 -0000	1.3
> --- defint.lisp	8 Jul 2002 15:35:43 -0000
> ***************
> *** 1492,1496 ****
>   
>   (DEFUN PERIOD (P E VAR)
> !        (ALIKE1 E (NO-ERR-SUB (m+ P VAR) E))) 
>   
>   (DEFUN INFR (A)
> --- 1492,1498 ----
>   
>   (DEFUN PERIOD (P E VAR)
> !   (and (ALIKE1 (no-err-sub var E) (setq e (NO-ERR-SUB (m+ P VAR) E)))
> !        ;; means there was no error
> !        (not (eq e t))))
>   
>   (DEFUN INFR (A)
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cut ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> 
> Wolfgang
>