Opened 11 years ago
Closed 17 months ago
#11590 closed defect (fixed)
Integrating the sgn() function can produce incorrect results
Reported by:  mjo  Owned by:  burcin 

Priority:  major  Milestone:  sage9.3 
Component:  calculus  Keywords:  
Cc:  nbruin, chapoton  Merged in:  
Authors:  Michael Orlitzky  Reviewers:  Frédéric Chapoton 
Report Upstream:  Fixed upstream, in a later stable release.  Work issues:  
Branch:  bf1b50b (Commits, GitHub, GitLab)  Commit:  bf1b50b7d9264e5695ec8cb900de253e6706136b 
Dependencies:  Stopgaps:  #12731 
Description
Actual result:
sage: integrate(x * sgn(x^2  1/4), x, 1, 0) 1/2
Since the argument to sgn() has only one root, 1/2, on (1, 0), there are only two subintervals on which sgn() can have different values. In particular,
sage: sgn(x^2  1/4)(x = 0.75) 1 sage: sgn(x^2  1/4)(x = 0.25) 1
Now, the original, actual result should be equivalent to the sum of the following:
sage: integrate(x, x, 1, 0.5) 0.375 sage: integrate(x, x, 0.5, 0) 0.125
So, something went wrong during the initial integration.
Change History (24)
comment:1 Changed 11 years ago by
 Type changed from PLEASE CHANGE to defect
comment:2 Changed 11 years ago by
 Report Upstream changed from N/A to Reported upstream. Little or no feedback.
comment:3 followup: ↓ 13 Changed 10 years ago by
 Report Upstream changed from Reported upstream. Little or no feedback. to Reported upstream. Developers acknowledge bug.
Even with #11483, this isn't working right. It should, though  see Barton's post at the Maxima ticket:
By the way: (%i4) load(abs_integrate)$ Correct antiderivative: (%i5) 'integrate(x*signum(x^21/4),x); (%o5) abs(x^21/4)/2 Correct definite integral (%i6) 'integrate(x*signum(x^21/4),x,1,0); (%o6) 1/4
so I'm not sure why this is still returning the "wrong" thing. Probably something about the integration code...
comment:4 followup: ↓ 8 Changed 10 years ago by
I'm stumped on this one. We get the correct antiderivative:
sage: integrate(x*sgn(x^21/4),x) 1/2*abs(x^2  1/4)
And ECL gives us the right answer, so it's not an environment setting:
sage: from sage.interfaces.maxima_lib import ecl_eval sage: ecl_eval("#$'integrate(x*signum(x^21/4),x,1,0);$") <ECL: ((RAT SIMP) 1 4)>
But going through maxima_eval
is still trouble:
sage: integrate(x*sgn(x^21/4),x,1,0) 1/2 sage: from sage.interfaces.maxima_lib import maxima_eval sage: a = '($INTEGRATE ((MTIMES SIMP) $X ((%SIGNUM SIMP) ((MPLUS SIMP) ((RAT SIMP) ( 1) 4) ((MEXPT SIMP) $X 2))) ) $X 1 0)' sage: maxima_eval(a) <ECL: ((RAT SIMP) 1 2)>
comment:5 followup: ↓ 6 Changed 10 years ago by
 Cc nbruin added
Well, I'm glad it wasn't just me!
Could it be that because we get an answer without abs_integrate
it just returns the 'regular' answer? But I don't recall the integrate code doing that, I think once we turn on abs_integrate
it should just 'be on'...
comment:6 in reply to: ↑ 5 Changed 10 years ago by
Replying to kcrisman:
Well, I'm glad it wasn't just me!
Could it be that because we get an answer without
abs_integrate
it just returns the 'regular' answer? But I don't recall the integrate code doing that, I think once we turn onabs_integrate
it should just 'be on'...
It's possible. The abs_integrate code looks like it defines extra definite integration methods, and then loops through them until 'integrate
is gone from the expression. Within maxima, the extra methods obviously get tried first, because we get the right answer. But I suppose something in the ECL integration could be trying the default definite integration first.
comment:7 Changed 10 years ago by
Something else is weird here. In standalone maxima5.26, we don't get back the wrong answer. But through sage, without abs_integrate, we get 1/2
:
sage: integrate(x * sgn(x^2  1/4), x, 1, 0) 1/2 sage: maxima.console() ;;; Loading #P"/home/mjo/src/sage5.0.beta1/local/lib/ecl/sbbsdsockets.fas" ;;; Loading #P"/home/mjo/src/sage5.0.beta1/local/lib/ecl/sockets.fas" ;;; Loading #P"/home/mjo/src/sage5.0.beta1/local/lib/ecl/defsystem.fas" ;;; Loading #P"/home/mjo/src/sage5.0.beta1/local/lib/ecl/cmp.fas" Maxima 5.26.0 http://maxima.sourceforge.net using Lisp ECL 11.1.1 Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) display2d:false; (%o1) false (%i2) 'integrate(x*signum(x^21/4),x); (%o2) 'integrate(x*signum(x^21/4),x) (%i3) 'integrate(x*signum(x^21/4),x,1,0); (%o3) 'integrate(x*signum(x^21/4),x,1,0)
Intriguing!
comment:8 in reply to: ↑ 4 Changed 10 years ago by
Some random observations that may or may not be relevant. First, let's separate the "Reader" (parser) and the "Evaluator".
sage: from sage.libs.ecl import * sage: I1=EclObject("#$integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car() sage: I2=EclObject("#$'integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car() sage: I1 <ECL: (($INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: I2 <ECL: ((%INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)>
Note the subtle difference between $INTEGRATE
(function) %INTEGRATE
(inert integral form). Sage's sr_integral
produces the $INTEGRATE
expression, so the two alternatives tested are *not* equivalent. I suspect those trigger different codepaths. You could try changing
interfaces/maxima_lib.py:205:max_integrate=EclObject("$INTEGRATE")
to max_integrate=EclObject("%INTEGRATE")
instead and see if that solves more than it messes up.
Slightly less likely to make a difference: If you look at the full expression resulting for EclObject("#$1+2$")
, you'll notice an MEVAL*
, whereas maxima_eval is defined with MEVAL
. I've never worked out which one is the proper one to use. So, you should test all combinations (these skip setting up the proper maximaerrorcatching environment):
sage: EclObject(["meval*",["quote", I1]]).eval() sage: EclObject(["meval*",["quote", I2]]).eval() sage: EclObject(["meval",["quote", I1]]).eval() sage: EclObject(["meval",["quote", I2]]).eval()
comment:9 followup: ↓ 10 Changed 10 years ago by
I took a much stupider approach, but I did essentially try replacing $INTEGRATE
with %INTEGRATE
. Once I discovered the (#$$)
syntax, I tried replacing the call to maxima_eval
in sr_integral
with a simple string substitution. It causes more problems than it solves, I think.
If I replace meval
with meval*
in a few places, I'm able to break some additional things but not fix any =)
I guess the Right Thing to do would be to go read the maxima source...
comment:10 in reply to: ↑ 9 Changed 10 years ago by
Replying to mjo:
I tried replacing the call to
maxima_eval
insr_integral
with a simple string substitution. It causes more problems than it solves, I think.
Yes, that would make things much worse (not to mention a major step back in other ways). The sr_integral
interface tries to avoid stringbased conversion, but maxima_lib in general is still capable of it in a much better way than crafting your own stringmangling based on #$...$
. For instance:
sage.calculus.calculus.maxima(sage.calculus.calculus.dummy_integrate(x * sgn(x^2  1/4), x, 1, 0))
should work and compares nicely to the pexpectbased
maxima(sage.calculus.calculus.dummy_integrate(x * sgn(x^2  1/4), x, 1, 0))
You should really try the commands I suggested earlier to see where the problem is, though (no changes to the source required!). I don't have ready access to a sage build with the right patches, but you obviously do. Even on a clearly insufficiently patched sage I observe different behaviour:
sage: maxima_eval(I1) <ECL: ((RAT SIMP) 1 2)> sage: maxima_eval(I2) RuntimeError: ECL says: Error executing code in Maxima: first: empty argument.
comment:11 followup: ↓ 12 Changed 10 years ago by
Sorry, I should have been more clear: I did try everything you suggested. Replacing $INTEGRATE
with %INTEGRATE
causes the same test failures in maxima_lib.py that I got when I swapped out the whole thing for string substitution. It does fix this particular answer, though.
Here's a session with just the abs_integrate patch applied, maxima5.26. I left the error at the top intact, for whatever reason it doesn't work at all until after I've integrated something.
sage: from sage.libs.ecl import * sage: I1=EclObject("#$integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car()  TypeError Traceback (most recent call last) /home/mjo/src/sage5.0.beta1/devel/sagedevel/<ipython console> in <module>() /home/mjo/src/sage5.0.beta1/local/lib/python2.7/sitepackages/sage/libs/ecl.so in sage.libs.ecl.EclObject.cadr (sage/libs/ecl.c:5528)() TypeError: cadr can only be applied to a cons sage: integrate(x*sgn(x^21/4),x,1,0) 1/2 sage: I1=EclObject("#$integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car() sage: I2=EclObject("#$'integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car() sage: I1 <ECL: (($INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: I2 <ECL: ((%INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: EclObject(["meval*",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval*",["quote", I2]]).eval() <ECL: ((RAT SIMP) 1 4)> sage: EclObject(["meval",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval",["quote", I2]]).eval() <ECL: ((RAT SIMP) 1 4)>
And here's a session with $INTEGRATE
switched to %INTEGRATE
:
sage: from sage.libs.ecl import * sage: I1=EclObject("#$integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car()  TypeError Traceback (most recent call last) /home/mjo/src/sage5.0.beta1/devel/sagedevel/<ipython console> in <module>() /home/mjo/src/sage5.0.beta1/local/lib/python2.7/sitepackages/sage/libs/ecl.so in sage.libs.ecl.EclObject.cadr (sage/libs/ecl.c:5528)() TypeError: cadr can only be applied to a cons sage: integrate(x*sgn(x^21/4),x,1,0) 1/4 sage: I1=EclObject("#$integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car() sage: I2=EclObject("#$'integrate(x*signum(x^21/4),x,1,0);$").cadr().cdr().car() sage: I1 <ECL: (($INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: I2 <ECL: ((%INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: EclObject(["meval*",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval*",["quote", I2]]).eval() <ECL: ((RAT SIMP) 1 4)> sage: EclObject(["meval",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval",["quote", I2]]).eval() <ECL: ((RAT SIMP) 1 4)>
The same thing, except the original integration actually works. Problem is, other tests start to fail with %INTEGRATE
. Here's an example:
$ sage t sage/interfaces/maxima_lib.py sage t "devel/sagedevel/sage/interfaces/maxima_lib.py" ********************************************************************** File "/home/mjo/src/sage5.0.beta1/devel/sagedevel/sage/interfaces/maxima_lib.py", line 661: sage: integral(x^n,x) Expected: Traceback (most recent call last): ... ValueError: Computation failed since Maxima requested additional constraints; using the 'assume' command before integral evaluation *may* help (example of legal syntax is 'assume(n+1>0)', see `assume?` for more details) Is n+1 zero or nonzero? Got: x^(n + 1)/(n + 1) ********************************************************************** File "/home/mjo/src/sage5.0.beta1/devel/sagedevel/sage/interfaces/maxima_lib.py", line 685: sage: integrate(sgn(x)  sgn(1x), x) Expected: abs(x  1) + abs(x) Got: (x  1)*sgn(x  1) + x*sgn(x) ********************************************************************** File "/home/mjo/src/sage5.0.beta1/devel/sagedevel/sage/interfaces/maxima_lib.py", line 700: sage: integrate(cos(x + abs(x)), x) Expected: 1/4*(sgn(x) + 1)*sin(2*x)  1/2*x*sgn(x) + 1/2*x Got: 1/4*(2*x  sin(2*x))*sgn(x) + 1/2*x + 1/4*sin(2*x) ********************************************************************** File "/home/mjo/src/sage5.0.beta1/devel/sagedevel/sage/interfaces/maxima_lib.py", line 711: sage: integral(abs(cos(x))*sin(x),(x,pi/2,pi)) Expected: 1/2 Got: 1/2
comment:12 in reply to: ↑ 11 Changed 10 years ago by
Replying to mjo:
I left the error at the top intact, for whatever reason it doesn't work at all until after I've integrated something.
sage: from sage.libs.ecl import *
This doesn't load maxima inside ecl yet, but once you evaluate an integral, maxima is loaded and the #$...$
macro works. Sorry about that.
comment:13 in reply to: ↑ 3 ; followup: ↓ 14 Changed 10 years ago by
Replying to kcrisman:
Even with #11483, this isn't working right. It should, though  see Barton's post at the Maxima ticket:
By the way: (%i4) load(abs_integrate)$ Correct antiderivative: (%i5) 'integrate(x*signum(x^21/4),x); (%o5) abs(x^21/4)/2 Correct definite integral (%i6) 'integrate(x*signum(x^21/4),x,1,0); (%o6) 1/4so I'm not sure why this is still returning the "wrong" thing.
Did you try it with integrate
rather than 'integrate
? Given what we've seen elsewhere in this ticket, I suspect it still gives the wrong answer. In sage, we are interfacing with integrate
. If that is still broken, then the bug is not fixed as far as sage is concerned.
If maxima's position is that we should use a different integration routine (i.e., 'integrate
) then we need a heuristic on when to use what routine ... isn't it maxima's job to figure this out?
comment:14 in reply to: ↑ 13 Changed 10 years ago by
Did you try it with
integrate
rather than'integrate
?
(%i5) display2d:false; (%o5) false (%i6) integrate(x*signum(x^21/4),x,1,0); (%o6) 1/2 (%i7) 'integrate(x*signum(x^21/4),x,1,0); (%o7) 1/4
Note that this was just Barton Willis (Maxima dev) example.
Given what we've seen elsewhere in this ticket, I suspect it still gives the wrong answer. In sage, we are interfacing with
integrate
. If that is still broken, then the bug is not fixed as far as sage is concerned.
I don't think they claimed it was fixed.
If maxima's position is that we should use a different integration routine (i.e.,
'integrate
) then we need a heuristic on when to use what routine ... isn't it maxima's job to figure this out?
Yes, I am a little stumped on this. I thought that the apostrophe just meant "don't evaluate nounform", but apparently I was mistaken. I just reported this example on the Maxima ticket, but I wouldn't hold my breath waiting for it, since there is a way to get it to work correctly there.
comment:15 Changed 10 years ago by
 Stopgaps set to #12731
comment:17 Changed 7 years ago by
 Stopgaps set to #12731
comment:18 Changed 6 years ago by
FYI fixed in Maxima by commit 5a300aa, which I have cherrypicked to branch5_38.
Related bug reports: https://sourceforge.net/p/maxima/bugs/2242 https://sourceforge.net/p/maxima/bugs/3123
comment:19 Changed 6 years ago by
 Report Upstream changed from Reported upstream. Developers acknowledge bug. to Fixed upstream, in a later stable release.
comment:20 Changed 3 years ago by
This now works in 8.9.b7. We need to add the doctest
sage: integrate(x * sgn(x^2  1/4), x, 1, 0) 1/4
comment:21 Changed 17 months ago by
 Branch set to u/mjo/ticket/11590
 Cc chapoton added
 Commit set to bf1b50b7d9264e5695ec8cb900de253e6706136b
 Status changed from new to needs_review
New commits:
bf1b50b  Trac #11590: add a doctest for the integral reported on this ticket.

comment:22 Changed 17 months ago by
 Reviewers set to Frédéric Chapoton
 Status changed from needs_review to positive_review
ok, this is now handled by giac..
comment:23 Changed 17 months ago by
 Milestone set to sage9.3
comment:24 Changed 17 months ago by
 Branch changed from u/mjo/ticket/11590 to bf1b50b7d9264e5695ec8cb900de253e6706136b
 Resolution set to fixed
 Status changed from positive_review to closed
This appears to be in Maxima, and is reported at their bug tracker.