[Maxima-commits] CVS: maxima/src defint.lisp,1.27,1.28

[Raymond Toy <rtoy@us...>]

Update of /cvsroot/maxima/maxima/src
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv940

Modified Files:
	defint.lisp 
Log Message:
Add a few comments about what intsc1, intsc, intsc0, and scprod are
trying to do.


Index: defint.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/defint.lisp,v
retrieving revision 1.27
retrieving revision 1.28
diff -u -d -r1.27 -r1.28
--- defint.lisp	1 Sep 2006 18:16:10 -0000	1.27
+++ defint.lisp	1 Sep 2006 20:48:17 -0000	1.28
@@ -1726,10 +1726,13 @@
 		  (not (period %pi2 e var)))
 	      (return nil))
 	     ((not (equal a 0.))
+	      ;; Apply the substitution to make the lower limit 0.
 	      (setq e (maxima-substitute (m+ a var) var e))
 	      (setq a 0.)
 	      (setq b limit-diff)))
-;;;Multiples of 2*%pi in limits.
+       ;; At this point, we have an integral with limits 0 and b.
+       
+       ;; Multiples of 2*%pi in limits.
        (cond ((eq (ask-integer (setq d (let (($float nil))
 					 (m// limit-diff %pi2))) '$integer)
 		  '$yes)
@@ -1742,10 +1745,16 @@
 	      (cond ((setq ans (intsc e %pi2 var))
 		     (return (m* d ans)))
 		    (t (return nil)))))
+       ;; The integral is not over a full period (2*%pi) or multiple
+       ;; of a full period.  Need to do something special.
        (cond ((ratgreaterp %pi2 b)
+	      ;; Less than 1 full period, so intsc can integrate it.
 	      (return (intsc e b var)))
-	     (t (setq l a) 
-		(setq a 0.)))
+	     (t
+	      (setq l a) 
+	      ;; Why do we need this?  I think if we get here, a is
+	      ;; already 0.
+	      (setq a 0.)))
        (setq b (infr b))
        (cond ((null l) 
 	      (setq nzp2 (car b))
@@ -1768,6 +1777,8 @@
 			    limit-diff))
        (return ans))))
 
+;; integrate(sc, var, 0, b), where sc is f(sin(x), cos(x)).  I (rtoy)
+;; think this expects b to be less than 2*%pi.
 (defun intsc (sc b var)
   (cond ((eq ($sign b) '$neg)
 	 (setq b (m*t -1. b))
@@ -1776,42 +1787,61 @@
   (cond ((setq b (intsc0 (cdr sc) b var))
 	 (m* (resimplify (car sc)) b))))
 
+;; integrate(sc, var, 0, b), where sc is f(sin(x), cos(x)).
 (defun intsc0 (sc b var)
+  ;; Determine if sc is a product of sin's and cos's.
   (let ((nn* (scprod sc))
 	(dn* ()))
-    (cond (nn* (cond ((alike1 b half%pi)
-		      (bygamma (car nn*) (cadr nn*)))
-		     ((eq b '$%pi)
-		      (cond ((eq (real-branch (cadr nn*) -1.) '$yes)
-			     (m* (m+ 1. (m^ -1 (cadr nn*)))
-				 (bygamma (car nn*) (cadr nn*))))))
-		     ((alike1 b %pi2)
-		      (cond ((or (and (eq (ask-integer (car nn*) '$even)
-					  '$yes)
-				      (eq (ask-integer (cadr nn*) '$even)
-					  '$yes))
-				 (and (ratnump (car nn*))
-				      (eq (real-branch (car nn*) -1.)
-					  '$yes)
-				      (ratnump (cadr nn*))
-				      (eq (real-branch (cadr nn*) -1.)
-					  '$yes)))
-			     (m* 4.	(bygamma (car nn*) (cadr nn*))))
-			    ((or (eq (ask-integer (car nn*) '$odd) '$yes)
-				 (eq (ask-integer (cadr nn*) '$odd) '$yes))
-			     0.)
-			    (t nil)))
-		     ((alike1 b half%pi3)
-		      (m* (m+ 1. (m^ -1 (cadr nn*)) (m^ -1 (m+l nn*)))
-			  (bygamma (car nn*) (cadr nn*))))))
-	  (t (cond ((and (or (eq b '$%pi)
-			     (alike1 b %pi2)
-			     (alike1 b half%pi))
-			 (setq dn* (scrat sc b)))
-		    dn*)
-		   ((setq nn* (antideriv sc))
-		    (sin-cos-intsubs nn* var 0. b))
-		   (t ()))))))
+    (cond (nn*
+	   ;; We have a product of sin's and cos's.  We handle some
+	   ;; special cases here.
+	   (cond ((alike1 b half%pi)
+		  ;; Wang p. 110, formula (1):
+		  ;; integrate(sin(x)^m*cos(x)^n, x, 0, %pi/2) =
+		  ;;   gamma((m+1)/2)*gamma((n+1)/2)/2/gamma((n+m+2)/2)
+		  (bygamma (car nn*) (cadr nn*)))
+		 ((eq b '$%pi)
+		  ;; Wang p. 110, near the bottom, says
+		  ;;
+		  ;; int(f(sin(x),cos(x)), x, 0, %pi) =
+		  ;;   int(f(sin(x),cos(x)) + f(sin(x),-cos(x)),x,0,%pi/2)
+		  (cond ((eq (real-branch (cadr nn*) -1.) '$yes)
+			 (m* (m+ 1. (m^ -1 (cadr nn*)))
+			     (bygamma (car nn*) (cadr nn*))))))
+		 ((alike1 b %pi2)
+		  (cond ((or (and (eq (ask-integer (car nn*) '$even)
+				      '$yes)
+				  (eq (ask-integer (cadr nn*) '$even)
+				      '$yes))
+			     (and (ratnump (car nn*))
+				  (eq (real-branch (car nn*) -1.)
+				      '$yes)
+				  (ratnump (cadr nn*))
+				  (eq (real-branch (cadr nn*) -1.)
+				      '$yes)))
+			 (m* 4.	(bygamma (car nn*) (cadr nn*))))
+			((or (eq (ask-integer (car nn*) '$odd) '$yes)
+			     (eq (ask-integer (cadr nn*) '$odd) '$yes))
+			 0.)
+			(t nil)))
+		 ((alike1 b half%pi3)
+		  ;; Wang, p. 111 says
+		  ;;
+		  ;; int(f(sin(x),cos(x)),x,0,3*%pi/2) =
+		  ;;   int(f(sin(x),cos(x)),x,0,%pi)
+		  ;;   + int(f(-sin(x),-cos(x)),x,0,%pi/2)
+		  (m* (m+ 1. (m^ -1 (cadr nn*)) (m^ -1 (m+l nn*)))
+		      (bygamma (car nn*) (cadr nn*))))))
+	  (t
+	   ;; We don't have a product of sin's and cos's.
+	   (cond ((and (or (eq b '$%pi)
+			   (alike1 b %pi2)
+			   (alike1 b half%pi))
+		       (setq dn* (scrat sc b)))
+		  dn*)
+		 ((setq nn* (antideriv sc))
+		  (sin-cos-intsubs nn* var 0. b))
+		 (t ()))))))
 
 ;;;Is careful about substitution of limits where the denominator may be zero
 ;;;because of various assumptions made.
@@ -1840,6 +1870,8 @@
 	  (t (let (($%piargs ()))
 	       (intsubs exp ll ul))))))
 
+;; Determine whether E is of the form sin(x)^m*cos(x)^n and return the
+;; list (m n).
 (defun scprod (e)
   (let ((great-minus-1 #'(lambda (temp)
 			   (ratgreaterp temp -1)))

[Maxima-commits] CVS: maxima/tests rtest15.mac,1.42,1.43

[Raymond Toy <rtoy@us...>]

Update of /cvsroot/maxima/maxima/tests
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv5643

Modified Files:
	rtest15.mac 
Log Message:
Add test from Bug 1547769:  integrate(sqrt(x^3/(2*a-x)),x,0,2*a)


Index: rtest15.mac
===================================================================
RCS file: /cvsroot/maxima/maxima/tests/rtest15.mac,v
retrieving revision 1.42
retrieving revision 1.43
diff -u -d -r1.42 -r1.43
--- rtest15.mac	1 Sep 2006 17:27:49 -0000	1.42
+++ rtest15.mac	1 Sep 2006 18:21:04 -0000	1.43
@@ -513,3 +513,15 @@
 
 integrate(1/(sin(x)^2+1),x,0,20*%pi);
 10*sqrt(2)*%pi;
+
+/*
+ * Bug  1547769:  integrate(sqrt(x^3/(2*a-x)),x,0,2*a);
+ *
+ * We don't get an internal error, and we should be able to evaluate
+ * this.  defint.lisp, rev 1.27
+ */
+(assume(a > 0), 0);
+0;
+
+integrate(sqrt(x^3/(2*a-x)),x,0,2*a);
+3*%pi*a^2/2;

[Maxima-commits] CVS: maxima/src defint.lisp,1.26,1.27

[Raymond Toy <rtoy@us...>]

Update of /cvsroot/maxima/maxima/src
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv3479

Modified Files:
	defint.lisp 
Log Message:
Modify BATA0 to recognize f(x)^y, where f(x) has the desired form:
x^kk*(b*x^n+a)^l.  This allows maxima to evaluate the integral in Bug
[ 1547769 ] integrate(sqrt(x^3/(2*a-x)),x,0,2*a).


Index: defint.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/defint.lisp,v
retrieving revision 1.26
retrieving revision 1.27
diff -u -d -r1.26 -r1.27
--- defint.lisp	1 Sep 2006 17:21:01 -0000	1.26
+++ defint.lisp	1 Sep 2006 18:16:10 -0000	1.27
@@ -1920,6 +1920,25 @@
 (defun bata0 (e)
   (let (k c)
     (cond ((atom e) nil)
+	  ((mexptp e)
+	   ;; We have f(x)^y.  Look to see if f(x) has the desired
+	   ;; form.  Then f(x)^y has the desired form too, with
+	   ;; suitably modified values.
+	   ;;
+	   ;; XXX: Should we ask for the sign of f(x) if y is not an
+	   ;; integer?  This transformation we're going to do requires
+	   ;; that f(x)^y be real.
+	   (destructuring-bind (mexp base power)
+	       e
+	     (declare (ignore mexp))
+	     (multiple-value-bind (kk cc)
+		 (bata0 base)
+	       (when kk
+		 ;; Got a match.  Adjust kk and cc appropriately.
+		 (destructuring-bind (l a n b)
+		     cc
+		   (values (mul kk power)
+			   (list (mul l power) a n b)))))))
 	  ((and (mtimesp e)
 		(null (cdddr e))
 		(or (and (setq k (findp (cadr e)))

[Maxima-commits] CVS: maxima/tests rtest15.mac,1.41,1.42

[Raymond Toy <rtoy@us...>]

Update of /cvsroot/maxima/maxima/tests
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv15409

Modified Files:
	rtest15.mac 
Log Message:
Add a few tests to verify that INTSC1 handles multiple periods
correctly.


Index: rtest15.mac
===================================================================
RCS file: /cvsroot/maxima/maxima/tests/rtest15.mac,v
retrieving revision 1.41
retrieving revision 1.42
diff -u -d -r1.41 -r1.42
--- rtest15.mac	8 Jul 2006 18:01:40 -0000	1.41
+++ rtest15.mac	1 Sep 2006 17:27:49 -0000	1.42
@@ -498,3 +498,18 @@
 
 limit (abs(x), x, 0);
 0;
+
+/* Related to Bug [1044318] defint(1/(sin(x)^2+1),x,0,3*%pi)
+ *
+ * integrate(1/(sin(x)^2+1),x,0,n*2*%pi) should be
+ * n*integrate(1/(sin(x)^2+1),x,0,2*%pi), for positive integer n.  We
+ * were just returning the same value.  
+ */
+integrate(1/(sin(x)^2+1),x,0,2*%pi);
+sqrt(2)*%pi;
+
+integrate(1/(sin(x)^2+1),x,0,4*%pi);
+2*sqrt(2)*%pi;
+
+integrate(1/(sin(x)^2+1),x,0,20*%pi);
+10*sqrt(2)*%pi;

[Maxima-commits] CVS: maxima/src defint.lisp,1.25,1.26

[Raymond Toy <rtoy@us...>]

Update of /cvsroot/maxima/maxima/src
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv12542/src

Modified Files:
	defint.lisp 
Log Message:
o Change CV so that we call $SUBSTITUTE instead of SUBST if the
  integration limit is a number.  This causes some simplification to
  be done, which works around a bug in Bug 1504505 that shows up in
  cmucl.  I think it happens because the expression isn't simplified
  but has lots of terms like x^n*0.  This needs some more
  work/thought.  Doesn't fix Bug 1504505

o Add a comment for what PQR does.

o In INTSC1, fix a bug where we didn't multiply the answer by the
  number of periods in the integrand.  Before
  integrate(1/(sin(x)^2+1),x,0,n*2*%pi) returned the same as
  integrate(1/(sin(x)^2+1),x,0,2*%pi), instead of n times that.



Index: defint.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/defint.lisp,v
retrieving revision 1.25
retrieving revision 1.26
diff -u -d -r1.25 -r1.26
--- defint.lisp	6 May 2006 13:44:06 -0000	1.25
+++ defint.lisp	1 Sep 2006 17:21:01 -0000	1.26
@@ -657,8 +657,14 @@
       (let ((trans (intcv3 (m// (m+t 'll (m*t 'ul var))
 				(m+t 1. var))
 			   nil 'yx)))
-	(setf trans (subst ll 'll trans))
-	(setf trans (subst ul 'ul trans))
+	;; If the limit is a number, use $substitute so we simplify
+	;; the result.  Do we really want to do this?
+	(setf trans (if (mnump ll)
+			($substitute ll 'll trans)
+			(subst ll 'll trans)))
+	(setf trans (if (mnump ul)
+			($substitute ul 'ul trans)
+			(subst ul 'ul trans)))
 	(method-by-limits trans var 0. '$inf))
       ()))
 
@@ -669,6 +675,9 @@
 	  (t (m+t (eezz (car e) ll ul)
 		  (cv (m// (cdr e) dn*)))))))
 
+;; I think this takes a rational expression E, and finds the
+;; polynomial part.  A cons is returned.  The car is the quotient and
+;; the cdr is the remainder.
 (defun pqr (e)
   (let ((varlist (list var)))
     (newvar e)
@@ -1724,9 +1733,15 @@
        (cond ((eq (ask-integer (setq d (let (($float nil))
 					 (m// limit-diff %pi2))) '$integer)
 		  '$yes)
+	      ;; This looks wrong.  We never multiply by d because of
+	      ;; the return!
+	      #+nil
 	      (setq ans (m* d (cond ((setq ans (intsc e %pi2 var))
 				     (return ans))
-				    (t (return nil)))))))
+				    (t (return nil)))))
+	      (cond ((setq ans (intsc e %pi2 var))
+		     (return (m* d ans)))
+		    (t (return nil)))))
        (cond ((ratgreaterp %pi2 b)
 	      (return (intsc e b var)))
 	     (t (setq l a)