From: Rupert Swarbrick <rswarbrick@us...>  20140729 19:46:12

This is an automated email from the git hooks/postreceive script. It was generated because a ref change was pushed to the repository containing the project "Maxima CAS". The branch, master has been updated via a7ccd8746ee3f5a6a0e4f32ebd6cfce0bd32fc76 (commit) via 417c1b7fba7e3c690677cdb7ad4e24f8b6a6d439 (commit) via f008c59cbb580bcb4b1e9b4744096a9420bd994d (commit) from 6db8907fce2cff0e9ef01fe660d3a6d647ade567 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below.  Log  commit a7ccd8746ee3f5a6a0e4f32ebd6cfce0bd32fc76 Author: Rupert Swarbrick <rswarbrick@...> Date: Tue Jul 29 20:45:34 2014 +0100 Correctly expand definite integrals in SP2INTEG2 The nontrivial code in SP2INTEG2 is for the case when we are trying to expand integrate (f(u), u, a(x), b(x)) in x. The previous code was just incorrect, as far as I can tell. The new version expands a(x) and b(x) in x. If the results are not monomials in x then we give up rather than trying to substitute arbitrary stuff into a power series. If they are monomials, we then expand the integrand (around zero: not sure about that bit). If that works, integrating should be easy, since it's just a power series. Next, substitute in the end points and subtract the results. Rather than adding another test, this fixes the test I added yesterday. Oops... diff git a/src/series.lisp b/src/series.lisp index 9643f79..89c7188 100644  a/src/series.lisp +++ b/src/series.lisp @@ 515,15 +515,58 @@ (t (sinint exp var)))) (t (sinint exp var))))) +;; Expand a definite integral, integrate(exp, v, lo, hi). (defun sp2integ2 (exp v lo hi)  (if (eq v var) (setq v (gensym) exp (subst v var exp)))  (cond ((and (free lo var) (free hi var))  (cond ((free exp var)  (list '(%integrate) (subst var v exp) var lo hi))  (t (sp3form (sp2expand exp)  (list '(%integrate) v lo hi)))))  (t (m+ (sp2sub (setq exp (sp2expand (subst var v exp))) hi)  (m* 1 (sp2sub exp lo)))))) + ;; Deal with the case where the integration variable is also the expansion + ;; variable. Something's a bit odd if we're doing this, but we assume that the + ;; aliasing was accidental and that the (bound) integration variable should be + ;; called something else. + (when (eq v var) + (setq v (gensym) + exp (subst v var exp))) + + (when (boundp '*sp2integrecursionguard*) + (throw 'psex + (list 'err '(mtext) + "Recursion when trying to expand the definite integral " + (outof (symbolvalue '*sp2integrecursionguard*))))) + + (cond + ((not (and (free lo var) (free hi var))) + ;; Suppose one or both of the endpoints involves VAR. We'll give up unless + ;; they are monomials in VAR (because substituting a nonmonomial into a + ;; power series is more difficult). If they *are* monomials in VAR, then we + ;; try to expand the integral as a power series and find an antiderivative. + (let ((loexp (sp2expand lo)) + (hiexp (sp2expand hi)) + (*sp2integrecursionguard* exp)) + (declare (special *sp2integrecursionguard*)) + + (unless (and (smono loexp var) (smono hiexp var)) + (throw 'psex + (list 'err '(mtext) + "Endpoints of definite integral " (outof exp) + " are not monomials in the expansion variable."))) + + ;; Since this is a formal powerseries calculation, we needn't concern + ;; ourselves with radii of convergence here. Just expand the integrand + ;; about zero. (Is there something cleverer we could do?) + (let ((antiderivative (sinint ($powerseries exp v 0) v))) + ;; If the expansion had failed, we'd have thrown an error to toplevel, + ;; so we can assume that ANTIDERIVATIVE is a sum of monomials in + ;; V. Substituting in LOEXP and HIEXP, we'll get two sums of + ;; monomials in VAR. The difference of the two expressions isn't quite + ;; of the right form, but hopefully it's close enough for any other + ;; code that uses it. + (m (maximasubstitute hiexp v antiderivative) + (maximasubstitute loexp v antiderivative))))) + + ((free exp var) + (list '(%integrate) (subst var v exp) var lo hi)) + + (t + (sp3form (sp2expand exp) + (list '(%integrate) v lo hi))))) ;;************************************************************************************ ;; phase three miscellaneous garbage and final simplification diff git a/tests/rtest16.mac b/tests/rtest16.mac index 224c4ac..d032bb2 100644  a/tests/rtest16.mac +++ b/tests/rtest16.mac @@ 2024,6 +2024,8 @@ sumcontract(intosum(powerseries(1+ (1x)^(a),x,0)  powerseries((1x)^(a),x,0))) /* #2775: Don't expand a powerseries with log(ab)=log(a)+log(b) if a is negative + + (and also #2767) */ (gensumnum : 0, powerseries (log(2x), x, 0)); sum(2^(i11)*x^(i1+1)/(i1+1), i1, 0, inf)$ +log(2)  sum (2^(i21)*x^(i2+1)/(i2+1), i2, 0, inf)$ commit 417c1b7fba7e3c690677cdb7ad4e24f8b6a6d439 Author: Rupert Swarbrick <rswarbrick@...> Date: Tue Jul 29 20:39:06 2014 +0100 Modify SP1LOG2 to generate a definite integral The previous code wrongly assumed that integrate(f'(x)/f(x), x) = log(f(x)) which is, of course, missing a constant of integration. Instead of generating an indefinite integral to be expanded, generate a definite integral. It turns out that SP2INTEG2, which is responsible for expanding definite integrals, is also broken, but that will be addressed in the next patch. diff git a/src/trgred.lisp b/src/trgred.lisp index db261ed..81e3f8d 100644  a/src/trgred.lisp +++ b/src/trgred.lisp @@ 510,24 +510,35 @@ ;; tries again after expressing E as integrate(diff(e)/e). ((sp1log2 e)))) +;; We didn't manage to expand the expression, so make use of the fact that +;; diff(log(f(x)), x) = f'(x)/f(x) and return integrate(f'(x)/f(x), x), hoping +;; that a later stage will be able to do something useful with it. +;; +;; We have to be a little bit careful because an indefinite integral might have +;; the wrong constant term. Instead, rewrite as +;; +;; log(f(x0+h)) = log(f(x0+h))  log(f(x0)) + log(f(x0)) +;; = integrate(diff(log(f(x0+k)), k), k, 0, h) + log(f(x0)) +;; = integrate(diff(f(x0+k))/f(x0+k), k, 0, h) + log(f(x0)) +;; +;; The "x0" about which we expand is always zero (see the code in $powerseries) (defun sp1log2 (e)  (and $verbose  (prog2  (mtell (intl:gettext "trigreduce: failed to expand.~%~%"))  (showexp (list '(%log) e))  (mtell (intl:gettext "trigreduce: try again after applying rule:~2%~M~%~%")  (list '(mlabel) nil  (outof  (list '(mequal)  (list '(%log) e)  (list '(%integrate)  (list '(mquotient)  (list '(%derivative) e var 1)  e)  var)))))))  (list '(%integrate)  (sp1 ($ratsimp (list '(mtimes) (sdiff e var) (list '(mexpt) e 1))))  var)) + (when $verbose + (mtell (intl:gettext "trigreduce: failed to expand.~%~%")) + (showexp (list '(%log) e)) + (mtell (intl:gettext "trigreduce: try again after applying rule:~2%~M~%~%") + (list '(mlabel) nil + (outof + `((mequal) + ((%log) ,e) + ((%integrate) + ((mquotient) ((%derivative) ,e ,var 1) ,e) ,var)))))) + (let* ((dummysym ($gensym))) + (m+ (list '(%log) ($limit e var 0)) + (list '(%integrate) + (maximasubstitute dummysym var + (sp1 (m// (sdiff e var) e))) + dummysym 0 var)))) (defun sp1trig (e) (cond ((atom (cadr e)) (simplify e)) commit f008c59cbb580bcb4b1e9b4744096a9420bd994d Author: Rupert Swarbrick <rswarbrick@...> Date: Mon Jul 28 23:54:05 2014 +0100 Tidy up expansions in EXPANDTRIGOFSUM It seemed a bit ugly having all the m+, m and m* written out four times, so I factored it out. Compiled code size goes from 2774 > 1108 bytes with x86_64 SBCL, so this might even be a tiny performance win but that wasn't really why I did it... diff git a/src/trgred.lisp b/src/trgred.lisp index b847f5e..db261ed 100644  a/src/trgred.lisp +++ b/src/trgred.lisp @@ 541,27 +541,15 @@ ;; Return the expansion of ((trigfun) ((mplus) a b)). For example sin(a+b) = ;; sin(a)cos(b) + cos(a)sin(b). (defun expandtrigofsum (trigfun a b)  (ecase trigfun  (%sin  (m+ (m* (sp1trig (list '(%sin) a))  (sp1trig (list '(%cos) b)))  (m* (sp1trig (list '(%cos) a))  (sp1trig (list '(%sin) b)))))  (%cos  (m (m* (sp1trig (list '(%cos) a))  (sp1trig (list '(%cos) b)))  (m* (sp1trig (list '(%sin) a))  (sp1trig (list '(%sin) b)))))  (%sinh  (m+ (m* (sp1trig (list '(%sinh) a))  (sp1trig (list '(%cosh) b)))  (m* (sp1trig (list '(%cosh) a))  (sp1trig (list '(%sinh) b)))))  (%cosh  (m+ (m* (sp1trig (list '(%cosh) a))  (sp1trig (list '(%cosh) b)))  (m* (sp1trig (list '(%sinh) a))  (sp1trig (list '(%sinh) b))))))) + (flet ((expandit (op f1 f2 f3 f4) + (funcall op + (m* (sp1trig (list f1 a)) (sp1trig (list f2 b))) + (m* (sp1trig (list f3 a)) (sp1trig (list f4 b)))))) + (ecase trigfun + (%sin (expandit #'add2* '(%sin) '(%cos) '(%cos) '(%sin))) + (%cos (expandit #'sub* '(%cos) '(%cos) '(%sin) '(%sin))) + (%sinh (expandit #'add2* '(%sinh) '(%cosh) '(%cosh) '(%sinh))) + (%cosh (expandit #'sub* '(%cosh) '(%cosh) '(%sinh) '(%sinh)))))) ;; Try to expand f(a+b) where f is sin, cos, sinh or cosh. (defun sp1trigex (e)  Summary of changes: src/series.lisp  59 ++++++++++++++++++++++++++++++++++++ src/trgred.lisp  75 ++++++++++++++++++++++++++ tests/rtest16.mac  4 ++ 3 files changed, 91 insertions(+), 47 deletions() hooks/postreceive  Maxima CAS 