Opened 5 years ago

Closed 4 years ago

## #10682 closed defect (fixed)

Reported by: Owned by: fmaltey burcin critical sage-5.0 symbolics maxima 5.26.0 binomial sum mjo sage-5.0.beta8 Dmitrii Pasechnik Jean-Pierre Flori, Nils Bruin N/A

I test

```sage: var ('n,k')
sage: sum (binomial(n,k)*k^2, k, 0, n)  # is right
sage: sum (binomial(n,k)*k^2, k, 1, n)  # is right n(n+1)2^(n-2)
sage: sum (binomial(n,k)*k^2, k, 2, n)  # is false : I get 0
```

This works correctly on Maxima 5.26 - we need to upgrade! The new spkg is at

Install the spkg and apply:

### comment:1 Changed 5 years ago by mboratko

• Description modified (diff)

### comment:2 Changed 4 years ago by was

• Priority changed from major to critical

This is not some general issue where sum fails on all inputs with lower bound != 0, 1. For example, this is fine:

```sage: var('n,k')
(n, k)
sage: sum(k, k, 1, n)
1/2*n^2 + 1/2*n
sage: sum(k, k, 2, n)
1/2*n^2 + 1/2*n - 1
```

Starting at 2, I guess this should be the answer:

```sage: (sum (binomial(n,k)*k^2, k, 1, n) - sum (binomial(n,k)*k^2, k, 1, 1))
(n^2 + n)*2^(n - 2) - n
```

The answer seems fine for starting at 3,4, by the way:

```sage: a = (sum (binomial(n,k)*k^2, k, 1, n) - sum (binomial(n,k)*k^2, k, 1, 2))
sage: b = sum (binomial(n,k)*k^2, k, 3, n)
sage: bool(a==b)
True
sage: a = (sum (binomial(n,k)*k^2, k, 1, n) - sum (binomial(n,k)*k^2, k, 1, 3))
sage: b = sum (binomial(n,k)*k^2, k, 4, n)
sage: bool(a==b)
True
```

This is probably a serious bug in Maxima, since the real work is done by the function

```sage.calculus.calculus.maxima.sr_sum
```

which just calls maxima in some complicated way.

I am raising the priority since this is a subtle and very serious mathematically incorrect result. Also, somebody needs to isolate this to a pure-maxima session exhibiting the bug (if it is in maxima!) and report upstream. I tried for a minute and failed.

### comment:3 Changed 4 years ago by nbruin

```Maxima 5.23.2 http://maxima.sourceforge.net
using Lisp ECL 11.1.1
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%o1) /usr/local/sage/4.7.1/local/share/maxima/5.23.2/share/contrib/solve_rec/simplify_sum.mac
(%i2) simplify_sum(sum(binomial(n,k)*k^2,k,2,n));
(%o2)                                  0
```

### comment:4 follow-up: ↓ 5 Changed 4 years ago by nbruin

Possibly the following ticket is relevant:

If so, it might be fixed. I don't have access to maxima 5.26 to test.

### comment:5 in reply to: ↑ 4 ; follow-up: ↓ 7 Changed 4 years ago by dimpase

• Status changed from new to needs_info

Possibly the following ticket is relevant:

If so, it might be fixed. I don't have access to maxima 5.26 to test.

It is fixed. (What's a big deal of having access to maxima 5.26? It compiles with ECL supplied with Sage 4.8):

```\$ ./maxima-local
Maxima 5.26.0_41_gbc6210e http://maxima.sourceforge.net
using Lisp ECL 11.1.1
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%o1)        /tmp/maxima/share/contrib/solve_rec/simplify_sum.mac
(%i2) simplify_sum(sum(binomial(n,k)*k^2,k,2,n));

2       n
(n  + n) 2  - 4 n
(%o2)                          -----------------
4
(%i3)
```

Thus, we need to update maxima spkg to 5.26.

By the way, I never understood why sum() is hardcoded this way in Sage. Sometimes I'd rather use Python's sum(), but I don't know how to do this.

### comment:6 Changed 4 years ago by dimpase

• Description modified (diff)
• Summary changed from sum fails with lower bound != 0 or 1 to sum fails with lower bound != 0 or 1 (upgrade maxima to 5.26)

### comment:7 in reply to: ↑ 5 Changed 4 years ago by nbruin

Sometimes I'd rather use Python's sum(), but I don't know how to do this.

This is off-topic but a good question. It shows people don't think of checking the python documentation for something like this, so perhaps a more prominent pointer is necessary

```sage: import __builtin__
sage: __builtin__.sum?
...
Docstring:
sum(sequence[, start]) -> value
```

(What's a big deal of having access to maxima 5.26? It compiles with ECL supplied with Sage 4.8)

I tried on sage.math.washington.edu with the 4.8 precompiled tarball there and it didn't work for me. Luckily there are other people more skilled at compiling maxima.

### comment:8 Changed 4 years ago by dimpase

• Description modified (diff)
• Status changed from needs_info to needs_review

### comment:9 Changed 4 years ago by dimpase

• Status changed from needs_review to needs_work

### comment:10 follow-ups: ↓ 12 ↓ 14 Changed 4 years ago by dimpase

• Description modified (diff)
• Work issues set to several doctests need to be patched due to changes in output format/term order

Several doctests need to be patched due to changes in output format/term order. See
the log.

### comment:12 in reply to: ↑ 10 ; follow-up: ↓ 13 Changed 4 years ago by dimpase

Several doctests need to be patched due to changes in output format/term order. See
the log.

The first test failure, in devel/sage/sage/matrix/matrix2.pyx, actually shows that 5.26 has improved in quality of such computations. The error between the exact and .n() answers has decreased by several orders of magnitude.

### comment:13 in reply to: ↑ 12 Changed 4 years ago by dimpase

The first test failure, in devel/sage/sage/matrix/matrix2.pyx, actually shows that 5.26 has improved in quality of such computations. The error between the exact and .n() answers has decreased by several orders of magnitude.

Here is the evidence (a is the matrix from the matrix/matrix2.pyx, line 10369, test, b is the exact matrix). The following is with the new spkg:

```sage: b=matrix([[1,pi],[100,1/100]])
sage: a=matrix(RR,[[1.,pi.n()],[100.,.01]])
sage: (a.exp()-b.exp()).norm()
9.32525865157e-07
```

The same computation on Sage 4.8 gives (a.exp()-b.exp()).norm() equal to 0.373383533611

### comment:14 in reply to: ↑ 10 ; follow-up: ↓ 15 Changed 4 years ago by mjo

Several doctests need to be patched due to changes in output format/term order. See
the log.

The failure in calculus/desolvers.py is what kept me from upgrading straight to 5.26 (I settled for 5.24 instead).

The issue is somewhere in our maxima interface... if you try the commands in a maxima console, you get the correct (simple) answer. If you try them with maxima._eval_line(), you get the mess from the doctest.

### comment:15 in reply to: ↑ 14 Changed 4 years ago by dimpase

Several doctests need to be patched due to changes in output format/term order. See
the log.

The failure in calculus/desolvers.py is what kept me from upgrading straight to 5.26 (I settled for 5.24 instead).

The issue is somewhere in our maxima interface... if you try the commands in a maxima console, you get the correct (simple) answer. If you try them with maxima._eval_line(), you get the mess from the doctest.

OK, this will need to be sorted out, and few other things. I am preparing the patch for most other, "trivial", failures, and will post it shortly.

### comment:16 follow-up: ↓ 20 Changed 4 years ago by nbruin

Thanks for preparing the spkg, Dima! My problem was probably that I didn't rebuild ecl first. It works very well now. I think I found the culprit for the desolve problem:

```Maxima 5.26.0 http://maxima.sourceforge.net
using Lisp ECL 11.1.1
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) display2d: false;

(%o1) false

(%o2) "/mnt/usb1/scratch/nbruin/5.0/local/share/maxima/5.26.0/share/contrib/diffequations/contrib_ode.mac"
(%i3) domain : complex;

(%o3) complex
(%i4) assume(x>0);

(%o4) [x > 0]
(%i5) assume(y>0);

(%o5) [y > 0]
(%i6) contrib_ode(x*'diff(y,x,1)-x*sqrt(y^2+x^2)-y,y,x);

(%o6) [(log((2*sqrt(4/x^2)*x^2*sqrt((4*y^4+4*x^2*y^2)/x^2)+8*y^2+4*x^2)/x^2)
-sqrt(4/x^2)*x*asinh(2*y^2/(x*sqrt(4*y^2)))
-sqrt(4/x^2)*x*asinh(2*y/sqrt(4*x^2))+sqrt(4/x^2)*x^2)
/(sqrt(4/x^2)*x)
= %c]
(%i7) domain: real;

(%o7) real
(%i8) contrib_ode(x*'diff(y,x,1)-x*sqrt(y^2+x^2)-y,y,x);

(%o8) [x-asinh(y/x) = %c]
```

so the problem is that we set "domain: complex" and apparently the behaviour of something there has changed. The assumptions x>0 and y>0 imply that x,y are real (and maxima's assume facility agrees with that), but as we are well aware, assumptions don't get used everywhere. At the very least, it seems that assumptions on variables do not affect all the things that "domain: complex" affects.

### comment:17 Changed 4 years ago by dimpase

there still is a bug in

```sage/functions/special.py", line 1464:
sage: elliptic_e(z, 1)
Expected:
sin(z)
Got:
sin(z) + 2*round(z/pi)
```

Note that BOTH expected and the received values are wrong!
The right answer is abs(sin(z)) + 2*round(z/pi)

Do you guys know how to report this upstream?

Dima

### comment:18 Changed 4 years ago by dimpase

to fix the problem in sage/interfaces/maxima_lib.py, printing extra (false), I added the patch mentioned
here to the spkg.

It fixes the problem. (I didn't want to mess with the unstable git version of Maxima, so we are applying the patch directly to
the stable Maxima 5.26.0).

### comment:19 Changed 4 years ago by dimpase

with the attached patch, and (updated) spkg, only one failure remains when running make ptestlong:

```sage -t  --long -force_lib devel/sage/sage/calculus/desolvers.py
**********************************************************************
File "/mnt/usb1/scratch/dima/sage-5.0.beta5/devel/sage-main/sage/calculus/desolvers.py", line 358:
sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True)
Expected:
[x - arcsinh(y(x)/x) == c]
Got:
[1/2*(2*x^2*sqrt(x^(-2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)/sqrt(x^2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x)) + log(4*(2*x^2*sqrt((x^2*y(x)^2 + y(x)^4)/x^2)*sqrt(x^(-2)) + x^2 + 2*y(x)^2)/x^2))/(x*sqrt(x^(-2))) == c]
**********************************************************************
1 of  59 in __main__.example_1
***Test Failed*** 1 failures.
```

This ought to be fixed by that complex/real domain thing described above by nbruin, but I don't know where this should be done.

### comment:20 in reply to: ↑ 16 Changed 4 years ago by dimpase

Thanks for preparing the spkg, Dima! My problem was probably that I didn't rebuild ecl first. It works very well now. I think I found the culprit for the desolve problem:

so the problem is that we set "domain: complex" and apparently the behaviour of something there has changed. The assumptions x>0 and y>0 imply that x,y are real (and maxima's assume facility agrees with that), but as we are well aware, assumptions don't get used everywhere. At the very least, it seems that assumptions on variables do not affect all the things that "domain: complex" affects.

Should we provide an option to do "domain: real" on the fly?
I don't know the details of the maxima interface, but I imagine this is doable...

Dima

### comment:21 Changed 4 years ago by dimpase

Here is the patch I propose for calculus/desolvers.py:

```68c68
< def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False, domain='complex'):
---
> def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False):
113,115d112
<     - ``domain`` - (optional) specifies the domain to use to solve the ODE.
<       By default it is complex numbers; setting domain='real' will make it real numbers.
<       It corresponds to Maxima's 'domain : complex;' or, respectively, 'domain : real;'
361,365d357
<         sage: desolve(x*diff(y,x)-x*sqrt(y^2+x^2)-y == 0, y, contrib_ode=True, domain='real')
<         [x - arcsinh(y(x)/x) == c]
<
<     Trac #10682 updated Maxima to 5.26, and it started to show a different solution in the complex domain for the ODE above::
<
367,368c359,360
<         [1/2*(2*x^2*sqrt(x^(-2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)/sqrt(x^2)) - 2*x*sqrt(x^(-2))*arcsinh(y(x)^2/(sqrt(y(x)^2)*x)) + log(4*(2*x^2*sqrt((x^2*y(x)^2 + y(x)^4)/x^2)*sqrt(x^(-2)) + x^2 + 2*y(x)^2)/x^2))/(x*sqrt(x^(-2))) == c]
<
---
>         [x - arcsinh(y(x)/x) == c]
>
428d419
<     P("domain : "+domain)
```

It gives deslove() one more optional parameter domain, which takes care of this new discrepancy with this ODE.
I'll update the whole patch shortly, test, and then expect the ticket to be ready for review.

### comment:22 Changed 4 years ago by dimpase

• Description modified (diff)
• Status changed from needs_work to needs_review

### comment:23 follow-up: ↓ 25 Changed 4 years ago by jpflori

Could you please post a diff between the new and old spkg's here?

Just for review of course.

### comment:24 follow-up: ↓ 26 Changed 4 years ago by jpflori

• Status changed from needs_review to needs_work

As far as the Sage patch is concerned, it would be nice to limit the length of code lines to 79 characters.

For example, in sage/symbolic/relation.py you corrected the output but put everything on one line, which is not very readable in 80 char terminals.

The error in f.nintegal() should definitely raise an Error. Letting it output of tuple of reals is wrong. I'll try to locate where this should take place if you have no time or no idea.

### comment:25 in reply to: ↑ 23 Changed 4 years ago by dimpase

Could you please post a diff between the new and old spkg's here?

Just for review of course.

Here it is (and of course I took the current maxima release 5.26.0-source from sf.net and replaced src/ with it, and then
ran the script spkg-dist)

```\$ hg diff -r tip -r 35
diff --git a/SPKG.txt b/SPKG.txt
--- a/SPKG.txt
+++ b/SPKG.txt
@@ -40,11 +40,6 @@

== Changelog ==

-=== maxima-5.26.0 (Dima Pasechnik, February 28th 2012) ===
- * upgrading to version 5.26.0 to take care of #10682.
- * added patch/comm.patch to fix Maxima bug #3484414.
-   (false) precedes the display output if display2d is set to false.
-
=== maxima-5.24.0.p0 (Michael Orlitzky, January 29th 2012) ===
* Trac #12094: Version bump to prevent a regression with the
abs_integrate package.
diff --git a/patches/comm.patch b/patches/comm.patch
deleted file mode 100644
--- a/patches/comm.patch
+++ /dev/null
@@ -1,12 +0,0 @@
---- src/src/comm.lisp
-+++ src/src/comm.lisp
-@@ -753,7 +753,9 @@
-        (setq ans (list '(mequal simp) (disp2 l) ans)))
-     (if lablist (nconc lablist (cons (elabel ans) nil)))
-     (setq tim (get-internal-run-time))
--    (displa (list '(mlable) (if lablist linelable) ans))
-+    (let ((*display-labels-p* nil))
-+      (declare (special *display-labels-p*))
-+      (displa (list '(mlable) (if lablist linelable) ans)))
-     (mterpri)
-     (timeorg tim)))
```

### comment:26 in reply to: ↑ 24 ; follow-up: ↓ 27 Changed 4 years ago by dimpase

As far as the Sage patch is concerned, it would be nice to limit the length of code lines to 79 characters.

For example, in sage/symbolic/relation.py you corrected the output but put everything on one line, which is not very readable in 80 char terminals.

oops, OK, I somehow had trouble splitting lines, and thought that the doctests framework doesn't quite allow this.
What is the syntax? Will it just ignore '\n', anywhere ?

The error in f.nintegal() should definitely raise an Error. Letting it output of tuple of reals is wrong. I'll try to locate where this should take place if you have no time or no idea.

Are you saying that this "error code"=6 is a new addition in Maxima, and is not properly handled by the Sage interface?

### comment:27 in reply to: ↑ 26 ; follow-up: ↓ 29 Changed 4 years ago by jpflori

oops, OK, I somehow had trouble splitting lines, and thought that the doctests framework doesn't quite allow this. What is the syntax? Will it just ignore '\n', anywhere ?

I think that you can more or less split where there is a space.

You can just copy what used to be done for the previous examples.

Are you saying that this "error code"=6 is a new addition in Maxima, and is not properly handled by the Sage interface?

I'm not sure that this is new, but:

• Sage used to return a ValueError? so I guess that it used to catch a Maxima error and translate it
• So it seems that Maxima's behavior has changed even though the "6" case used to be there, maybe a Maxima expert could tell us why this was decided. Anyhow, your fix might be the right one then...

A last comment, in your modification to desolve() you modify the domain var of Maxima.

So if I use domain=real with your code, then all of Maxima will use domain=real afterwards.

Maybe we'd better set back domain to whatever it was worth before at the end of the computation.

### comment:28 Changed 4 years ago by jpflori

By the way the upstream tarball does not seem to include the latest changelogs.

### comment:29 in reply to: ↑ 27 Changed 4 years ago by kcrisman

A last comment, in your modification to desolve() you modify the domain var of Maxima.

So if I use domain=real with your code, then all of Maxima will use domain=real afterwards.

Maybe we'd better set back domain to whatever it was worth before at the end of the computation.

Yes, definitely. I'd be surprised if this didn't cause doctest failures otherwise, since we have domain:complex because it made other things work in the past. Otherwise this is a good solution, I think, and gives the user more access to Maxima's flexibility without too much effort, which is a good thing.

### comment:30 Changed 4 years ago by jpflori

What I told was wrong for the "6" stuff.

In fact, Maxima used to return the expression unmodified when the precision was too high, now it behaves correctly, so I'm happy with your fix.

However, we could rid of some of the error handling in nintegral and mention (and before that find...) the corresponding commit in maxima git tree.

### comment:31 Changed 4 years ago by mjo

I think it would be nice if we handled this automatically rather than requiring the user to know about the internals of solution/simplification. I mean, none of us knew what the problem was until Nils pointed it out.

My understanding is that we default to domain: complex because it will make the fewest invalid simplifications and is compatible with pynac. In general, that's the right thing to do. But here? We know that all of the variables in the system are real, from the assumption that x,y>0.

Why not just set the domain to real in that case?

(I think we should add something like a set_domain_real method to maxima_lib either way to make this cleaner.)

### comment:32 follow-ups: ↓ 33 ↓ 40 Changed 4 years ago by jpflori

If we just rely on assumptions of the form 'x > 0', it should not be too complicated to check that all variables are supposed to be real.

Anyway, isn't there something in Maxima which "knows" that a variable is real after such assumptions ?

If that is the case, its current behavior could be considered quite strange.

### comment:33 in reply to: ↑ 32 Changed 4 years ago by kcrisman

If we just rely on assumptions of the form 'x > 0', it should not be too complicated to check that all variables are supposed to be real.

Anyway, isn't there something in Maxima which "knows" that a variable is real after such assumptions ?

If I understand correctly, Maxima and assumptions play weakly together. They are ignored in a lot of "dummy variable" contexts like solving and integration. So maybe here too.

Interestingly, our GSL numerical ODE solvers seem to only have real domains (see this ask.sagemath.org question, I can't verify this independently), so maybe we could instead just always change to domain:real for ODE and then convert back when done... these are all hacks, but some hacks are better than others.

### comment:34 Changed 4 years ago by mjo

For the assumptions, we don't necessarily need maxima to recognize them. If the user tells us that x>0, then we should be able to set domain: real whether or not maxima is going to use that information.

In simplification, I gather that maxima doesn't really know about real/complex. The docs simply say,

```Option variable: domain

Default value: real

When domain is set to complex, sqrt (x^2) will remain sqrt (x^2) instead of returning abs(x).
```

```The notion of a "domain" of simplification is still in its infancy, and controls little more than this at the moment.
```

I think we should leave the default complex to be on the safe side. We could strictly improve some things by checking the assumptions, though. There are a few easy cases that we should be able to handle safely:

1. assume(x, 'real')
2. assume(x > constant)
3. assume(x > y) (where y is known to be real by similar reasoning)
4. et cetera

If all variables in an expression are determined to be real, we can set domain: real before asking maxima to do anything with it.

### comment:35 follow-up: ↓ 36 Changed 4 years ago by kcrisman

Good points. That being said, if the DE is one that is more complicated, and some nasty expression plays the role that x does here, we don't want to try to figure out if that expression is real or not.

### comment:36 in reply to: ↑ 35 ; follow-up: ↓ 37 Changed 4 years ago by dimpase

Good points. That being said, if the DE is one that is more complicated, and some nasty expression plays the role that x does here, we don't want to try to figure out if that expression is real or not.

I thing we should not try to do our own processing of assume()'s at this point. On another ticket, maybe.

Regarding the domain switch, I simply could not find a place to hook it up, as an independent call, instead opting for desolve() to have an extra parameter. If you read desolve() code you see that it digs up the maxima instance to call, as the parent of the expression. I could not find an independent, "global" hook to call maxima, and it seems that it might not even be exposed. It seems that the original design aimed at possibly many independent maxima instances to be run.

By the way, the call of desolve with domain='complex', (or just with domain omitted) will reset the domain back to complex. As such, this gives one a way to control the domain, something that was missing in the maxima interface before (although it is a hack).

Any thoughts on this?

### comment:37 in reply to: ↑ 36 ; follow-up: ↓ 38 Changed 4 years ago by kcrisman

Regarding the domain switch, I simply could not find a place to hook it up, as an independent call, instead opting for desolve() to have an extra parameter. If you read desolve() code you see that it digs up the maxima instance to call, as the parent of the expression. I could not find an independent, "global" hook to call maxima, and it seems that it might not even be exposed. It seems that the original design aimed at possibly many independent maxima instances to be run.

Yes, and still is. But you should be able to import sage.calculus.calculus.maxima and the like.

```sage: sage.calculus.calculus.maxima('domain')
complex
sage: sage.calculus.calculus.maxima('domain:real')
real
sage: sage.calculus.calculus.maxima('domain')
real
sage: sage.calculus.calculus.maxima('domain:complex')
complex
sage: sage.calculus.calculus.maxima('domain')
complex
```

### comment:38 in reply to: ↑ 37 Changed 4 years ago by dimpase

Regarding the domain switch, I simply could not find a place to hook it up, as an independent call, instead opting for desolve() to have an extra parameter. If you read desolve() code you see that it digs up the maxima instance to call, as the parent of the expression. I could not find an independent, "global" hook to call maxima, and it seems that it might not even be exposed. It seems that the original design aimed at possibly many independent maxima instances to be run.

Yes, and still is. But you should be able to import sage.calculus.calculus.maxima and the like.

I think I experimented with something like this, and switching domain this way had no effect! But I might be wrong. I'll recheck. We'll be travelling NZ-Singapore in the coming 16 hours, so it would take time.

### comment:39 Changed 4 years ago by mjo

We might as well just add domain as a @property on maxima_lib. That way we can do,

```old_domain = maxima.domain
maxima.domain = 'real'
result = perform_things()
maxima.domain = old_domain
return result
```

### comment:40 in reply to: ↑ 32 ; follow-up: ↓ 41 Changed 4 years ago by nbruin

Anyway, isn't there something in Maxima which "knows" that a variable is real after such assumptions ?

It's just a little weak in its deductions:

```(%i1) domain: complex;
(%o1)                               complex
(%i2) display2d: false;
(%o2) false
(%i3) sqrt(x^2);
(%o3) sqrt(x^2)
(%i4) assume(x>0);
(%o4) [x > 0]
(%i5) sqrt(x^2);
(%o5) x
(%i6) sqrt(4/x^2);
(%o6) sqrt(4/x^2)
(%i7) domain:real;
(%o7) real
(%i8) sqrt(4/x^2);
(%o8) 2/x
```

It is ingrained now that the calculus maxima runs with domain: complex, so if you want to temporarily switch to domain: real, you should switch back afterwards. In principle this could play havoc with multiple threads, but maxima is not threadsafe anyway and our signal management isn't compatible with ECL's ideas of dealing with threads either. So, if in desolve (conditional on domain="real") we just do:

```P("load('contrib_ode)")
...
try:
P("domain: real")
soln = P(cmd)
finally:
P("domain: complex")
```

we should be fine (the "try ... finally" should guard that we end up with domain: complex). I don't think performance is a concern if we're doing a "load" every time anyway!.

### comment:41 in reply to: ↑ 40 Changed 4 years ago by jpflori

ideas of dealing with threads either. So, if in desolve (conditional on domain="real") we just do: P("load('contrib_ode)") ... try: P("domain: real") soln = P(cmd) finally: P("domain: complex") we should be fine (the "try ... finally" should guard that we end up with domain: complex). I don't think performance is a concern if we're doing a "load" every time anyway!.

That's the kind of solution I would propose.

As far as the Maxima instance is concerned, you get it through the ._maxima_().parent() call in nintegral(...) and it happens to be the unique library instance of Maxima (which is currently used for all the calculus stuff). In particular, a call to nintegral should not modify its settings.

To refine Nils answer, Id say that in try we set domain to the option passed and in "finally" we should restore the previous value of domain which could be "real"rather than "complex" if the user decided to change it globally.

To get the previous value something like str(P('domain').strip(' ') would do at the Sage level although this is kind of dirty (and surely slow)

### comment:42 Changed 4 years ago by dimpase

• Status changed from needs_work to needs_review

please review the revised patch. I have fixed the formatting problems and removed the domain parameter (which I added in the previous version of the patch) from desolve().

Instead, I simply call sage: sage.calculus.calculus.maxima('domain:real')
and sage: sage.calculus.calculus.maxima('domain:complex') in the docstest.
This works just as well.

### comment:43 follow-ups: ↓ 45 ↓ 52 Changed 4 years ago by jpflori

• Authors set to Dima Pasechnick
• Keywords maxima 5.26.0 binomial sum added
• Reviewers set to Jean-Pierre Flori,
• Status changed from needs_review to needs_work
• Work issues changed from several doctests need to be patched due to changes in output format/term order to domain issue, error handling in nintegral

I would rather follow Nils approach and at least include a try/finally statement in cas an error is raised.

Doing so, we're sure not to let the domain set to real after leaving the function because of an error.

Then, you'd better use "P" already present in the code rather than sage.calculus.calculs.maxima so that you're sure to use the same maxima as the one that will be used for solving the equation. Currently it's equivalent, but if something else changes in the future, such a choice will make everything easier.

Finally, on my side, I'd really prefer first to use the option approach you already implemented, and then to bakup the previous value of domain and set it back in the finally statement (rather than setting it to complex all the time).*

If you want to do it more cleanly, you can use

```def desolve(..., domain='real')
...
olddomain=P.get('domain')
P.set('domain',domain)
...
try:
...
finally:
P.set('domain',olddomain)
...

```

rather than P('domain :'+str)

And the error handling at the end of nintegral should definitely be reworked.

For example the test "if 'quad_qags' in str(v)" is useless now that Maxima behavior has changed.

Moreover, I'm not sure the "ERROR in srt(err)" can ever be reached now (was it before???).

I could not locate the commit corresponding to this change although I tried a little hard. Maxima source tree is a little obscure for me. This commit should be mentionned in the doc and would help to answer the two above remarks about error handling. Some expert of Maxima is needed for that.

A last question is to know what we do when nintegral returns an error code.

Do we let the situation as is? Or do we raise an error. I would agrre that this should be treated in another ticket, but I also don't think there is a better choice to make. Some discussion should happen.

### comment:44 Changed 4 years ago by mjo

The whole

```desolve(..., make_it_work=True)
```

approach is silly. If we don't get a solution without it, we should use contrib_ode. If we know the variables are real, we should set domain=real.

It might require a second ticket and some work on MaximaLib?/desolve, but we shouldn't just paper over the fact that this doesn't work (and is a regression in that regard). Requiring the user to manually set the domain is a regression.

### comment:45 in reply to: ↑ 43 ; follow-up: ↓ 46 Changed 4 years ago by nbruin

And the error handling at the end of nintegral should definitely be reworked.

For example the test "if 'quad_qags' in str(v)" is useless now that Maxima behavior has changed.

Moreover, I'm not sure the "ERROR in srt(err)" can ever be reached now (was it before???).

I could not locate the commit corresponding to this change although I tried a little hard. Maxima source tree is a little obscure for me. This commit should be mentionned in the doc and would help to answer the two above remarks about error handling. Some expert of Maxima is needed for that.

I did find this commit in the history of src/numerical/slatec/quadpack.lisp. It seems to indicate that something about quadpack error reporting has changed and may have become more configurable. Perhaps we can configure it in a way that the reported errors are more suitable and detectable for us?

desolve: I think Dima's new approach is acceptable in that it both fixes the doctest and documents the regression. Did anybody ever check if the returned solution with domain: complex is actually correct with some crazy combination of branch cut choices?

### comment:46 in reply to: ↑ 45 Changed 4 years ago by mjo

desolve: I think Dima's new approach is acceptable in that it both fixes the doctest and documents the regression. Did anybody ever check if the returned solution with domain: complex is actually correct with some crazy combination of branch cut choices?

It fixes the doctest in the same way that deleting the doctest fixes the doctest.

All we have to do to fix this ticket is backport simplify_sum.mac:

```sage: maxima.version()
'5.24.0'
sage: n,k = var('n,k')
sage: sum(binomial(n,k)*k^2, k, 2, n)
1/4*(n^2 + n)*2^n - n
```

If we're going to upgrade maxima after that, we should do it properly, i.e. no regressions.

### comment:47 Changed 4 years ago by jpflori

I saw that commit but I must admit I do not understand whether that actually modified the behavior of Maxima and slatec or not.

My point here is that we should get rid of the now useless lines of code and if possible document it properly either in the file itself or in SPKG.txt with a pointer to the corresponding Maxima commit and some explanation.

Deciding whether the current behavior of letting the user know that an error occurred through the last value of the tuple (as maxima does) which is already document, or whether we should actually look at that value in Sage itself and raise an actual Error accordingly can be done later (or not, I'll personally live with it as it is).

### comment:48 Changed 4 years ago by jpflori

I also quite agree with mjo.

If possible, let's just update the 5.25 spkg with an additional patch for the sum (as Dima just did for the false stuff), so that it is quickly reviewed and merged.

And let's upgrade maxima in another ticket (and another one for error handling, and potentially another one for a real solution to the domain stuff, my preferred solution for that is the domain option with try/finally and real by default, for mjo I'm not sure we can call it a regression, did Sage (and Maxima?) actually outputed something different for domain equal to real or complex before? or was it always as it is now with domain= real?).

### comment:49 follow-ups: ↓ 50 ↓ 54 Changed 4 years ago by mjo

It always returned,

```[x - arcsinh(y(x)/x) == c]
```

regardless of the domain.

But, it used to give the correct answer, and now it doesn't.

I'm running a ptestlong now with simplify_sum.mac from 5.26. If that (and the maxima tests) pass, I would certainly give a patched 5.24 a positive review.

### comment:50 in reply to: ↑ 49 ; follow-up: ↓ 51 Changed 4 years ago by dimpase

• Status changed from needs_work to needs_info

It always returned, [x - arcsinh(y(x)/x) == c]

regardless of the domain.

But, it used to give the correct answer, and now it doesn't.

I'm running a ptestlong now with simplify_sum.mac from 5.26. If that (and the maxima tests) pass, I would certainly give a patched 5.24 a positive review.

Can you demonstrate in some way that the solution it gives in the complex domain is wrong? If yes, you have a point. Still, if it were wrong we would have a good Maxima bug report to file. And we have a workaround for this, which is good enough IMHO.

One way or another, I would much prefer have a patch in 5.26 than in 5.24. The latter is an obsolete version, and I see very little value in keeping it (I am sure they fixed a zillion other bugs in 5.26, which were not covered by Sage doctests). Just from the point of view of maintaining good relations with Maxima community it is important to upgrade to 5.26.

### comment:51 in reply to: ↑ 50 Changed 4 years ago by dimpase

(I am sure they fixed a zillion other bugs in 5.26, which were not covered by Sage doctests).

E.g. see my comment above about the quality of numerical approximation in
matrix/matrix2.pyx, line 10369, test, as an example of improvement. Actually, the quality of the approximation returned there by the current maxima spkg is so poor that it can be regarded as a bug!

### comment:52 in reply to: ↑ 43 Changed 4 years ago by dimpase

I would rather follow Nils approach and at least include a try/finally statement in cas an error is raised.

Doing so, we're sure not to let the domain set to real after leaving the function because of an error.

Then, you'd better use "P" already present in the code rather than sage.calculus.calculs.maxima so that you're sure to use the same maxima as the one that will be used for solving the equation. Currently it's equivalent, but if something else changes in the future, such a choice will make everything easier.

I really cannot imagine a scenario when one would run several instances of maxima in a totally uncoordinated way. If such a major change was ever to happen, surely a facility to execute a series of maxima commands in a given maxima instance would be available, and could be used then.

Finally, on my side, I'd really prefer first to use the option approach you already implemented, and then to bakup the previous value of domain and set it back in the finally statement (rather than setting it to complex all the time).*

That's a good point, regardless. As I said, when working on the previous patch I overlooked a way to set domain by another call to maxima.

[...]

A last question is to know what we do when nintegral returns an error code.

Do we let the situation as is? Or do we raise an error. I would agrre that this should be treated in another ticket, but I also don't think there is a better choice to make. Some discussion should happen.

It surely can be handled better, but it surely should happen on another ticket.

### comment:53 Changed 4 years ago by mjo

I would be happy with the following:

1. Add a public method to set the maxima simplification domain. This is useful elsewhere, and is pain-equivalent to adding a domain parameter to the desolve function. Can be a separate ticket.
2. Open another ticket to figure out the domain automatically in desolve where all the variables are known.
3. Update this doctest to use that method, and add a reference to that ticket.

### comment:54 in reply to: ↑ 49 Changed 4 years ago by nbruin

But, it used to give the correct answer, and now it doesn't.

As far as I can check, the newly returned solution is valid:

```sage: y=function('y',x)
sage: sol=(log((2*sqrt(4/x^2)*x^2*sqrt((4*y^4+4*x^2*y^2)/x^2)+8*y^2+4*x^2)/x^2)
....:         -sqrt(4/x^2)*x*asinh(2*y^2/(x*sqrt(4*y^2)))
....:         -sqrt(4/x^2)*x*asinh(2*y/sqrt(4*x^2))+sqrt(4/x^2)*x^2)/(sqrt(4/x^2)*x)
sage: var('DY')
DY
sage: eqn=sol.diff(x).subs_expr(diff(y,x)==DY)
sage: dydx=-eqn.coefficient(DY,0)/eqn.coefficient(DY,1)
sage: assume(x>0)
sage: E=x*dydx-x*sqrt(y^2+x^2)-y
sage: E.factor()
0
```

so I don't think we need to consider maxima 5.26's response a regression in this respect. It is simply returning a different solution than 5.24 did when domain: complex is in effect.

### comment:55 follow-up: ↓ 57 Changed 4 years ago by jpflori

About the multiplie maxima instance, there is indeed only one library instance and I don't think that will change in the near future.

But one can instantiate pexpect interfaces (one or multiple one).

So if I'm a little bizarre and want to use such an interface, I'd like it to behave as much as possible as the library one (e.g. doing something like de_max = maxima(de) and then calling desolve on de_max, in that case we'll get a pexpect interface by the parent() call)..

Of course desole is not really implemented to work in that case and in fact it does not seem to do very well right now (it just returns different kinds of NotImplementedError? telling that Maxima failed).

As a consequence we might also think of getting rid of the _maxima_() and parent() call and explicitely only use sage.calculus.calculus.maxima instead.

In that case I'll be completely happy with the use of sage.calculus.calculus.maxima for setting the domain.

Setting the domain can be domain with a simple call to set('domain','xxx') as mentioned above.

But I just discovered that we (I?) messed up something when building the library interface.

With the pexpect interface (be sure to actually start it with maxima.start()), one can use maxima.domain directly. This does not work with the library interface.

Of course we can also add a new method on top of that as well.

### comment:56 Changed 4 years ago by kcrisman

• Authors changed from Dima Pasechnick to Dima Pasechnik

Disagree with having the other interface instances behaving like the standard one. On the contrary, some people will probably want "vanilla" Maxima in that setting. That said, there should be a nice easy way to say "pexpect interface with the same defaults as the "calculus" copy" or something like that.

### comment:57 in reply to: ↑ 55 ; follow-ups: ↓ 58 ↓ 59 Changed 4 years ago by dimpase

[...]

Setting the domain can be domain with a simple call to set('domain','xxx') as mentioned above.

But I just discovered that we (I?) messed up something when building the library interface.

With the pexpect interface (be sure to actually start it with maxima.start()), one can use maxima.domain directly. This does not work with the library interface.

Of course we can also add a new method on top of that as well.

OK, I am confused now - can we have a new public method to do the simplification domain switching in the library, i.e.
the equivalent of sage.calculus.calculus.maxima('domain:real'), etc., or it would not work due to a bug in the library interface? If the latter, it looks like one can give it the positive review now, right?

I am relieved to hear that there are "real" regressions --- bedankt, Nils!

### comment:58 in reply to: ↑ 57 Changed 4 years ago by dimpase

I am relieved to hear that there are "real" regressions --- bedankt, Nils!

I meant to say "there are no 'real' regressions", of course...

### comment:59 in reply to: ↑ 57 Changed 4 years ago by jpflori

OK, I am confused now - can we have a new public method to do the simplification domain switching in the library, i.e. the equivalent of sage.calculus.calculus.maxima('domain:real'), etc., or it would not work due to a bug in the library interface? If the latter, it looks like one can give it the

You're also confusing me now :)

I really do not like and understand how trac treats newline stuff, whence my maybe confusing message above. I'll try to format that one better.

I'll also says a lot a lot of trivialities, but it's to be sure that we all agree on what we want or do not want.

Fisrt, sage.calculus.calculus.maxima('domain: real') does work as intended indeed.

If you think this is not public enough we can add a new method "set_domain" in MaximaAbstract? class which does the same thing:

```def set_domain(self, domain):
r"""
Sets the working domain of this Maxima instance to `domain`.
"""
self('domain : "+domain)

```

or equivalently

```def set_domain(self, domain):
r"""
Sets the working domain of this Maxima instance to `domain`.
"""
self.set('domain', domain)

```

with some basic doctesting.

In that form I doubt this is really much more useful than calling sage.calculus.calculus.maxima(...).

Maybe you meant a function declared at the top level of sage.calculus.calculus which does the same ?

It's maybe a little more meaningful.

What I meant about the Maxima Library interface bug is that we should be able to access the "domain" variable of Maxima directly as an attribute of the Sage object, so that one could replace the above method by

```def set_domain(self, domain):
...
self.domain = domain

```

or if its a function in sage.calculus.calculus :

```def set_domain(domain):
r"""
Sets the domain of the Maxima instance used by Sage for calculus to `domain`
""""
maxima.domain = domain

```

This works for the pexpect interface but is broken for the library interface.

### comment:60 follow-ups: ↓ 61 ↓ 63 Changed 4 years ago by jpflori

And I do not think that this should be positively reviewed yet: the error handling changes of Maxima in nintegral should be documented and the now useless lines removed.

### comment:61 in reply to: ↑ 60 ; follow-up: ↓ 62 Changed 4 years ago by nbruin

And I do not think that this should be positively reviewed yet: the error handling changes of Maxima in nintegral should be documented and the now useless lines removed.

Indeed, if you read this thread it seems that in maxima 5.16.2 they decided to return the expression in case of error in evaluation and that this decision was reversed, so now they're just reporting the error codes as computed by quadpack. Those error codes are documented in the sage.calculus.calculus.nintegral docstring, so probably the lines

```    #This is just a work around until there is a response to
#http://www.math.utexas.edu/pipermail/maxima/2008/012975.html
raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
```

could just be removed. A workaround is not necessary anymore. If we prefer to raise errors rather than rely on checking error codes, we could do

```    errorcode_dict = {
1 : 'too many sub-intervals were done',
2 : 'excessive roundoff error detected',
3 : 'extremely bad integrand behavior occurs',
4 : 'failed to converge',
5 :  'integral is probably divergent or slowly convergent',
6 : 'invalid input'} #this def can go outside function
cond=int(v[3])
if cond in errorcode_dict:
raise RuntimeError("Maxima (via quadpack) reports an error:'%s'"%error_code_dict(cond)
else:
return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
```

However, once you start changing things like that, we should perhaps change the return value format (or make it a parameter 'raise_errors_rather_than_return_error_codes').

### comment:62 in reply to: ↑ 61 Changed 4 years ago by nbruin

```    #This is just a work around until there is a response to
#http://www.math.utexas.edu/pipermail/maxima/2008/012975.html
raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
```

My apologies! I've looked at the quadpack.lisp source and as you might notice, there's a handler-case ... error construct there, which is a lisp version of try ... except. The error clause constructs a return value which would be an unevaluated quad_qags call in maxima. So the raise above could in principle still be triggered (e.g., give an integrand that triggers an error when evaluated in a way that's not caught otherwise). We need to leave this in. One could leave out the try...except around the maxima call if the error messages produced by maxima are sufficiently informative.

The patch that triggered Mike Hansen's request and Dodier's answer is probably this one. That code is still there as far as I can see. I think we're best off leaving both error checks in there. The change in behaviour between 5.24 to 5.26 is likely caused somewhere else. The 5.26 behaviour is still as documented (provided we consider 1e-14 as 'invalid input').

Whether we want to raise an error rather than return an error code is an independent question.

### comment:63 in reply to: ↑ 60 Changed 4 years ago by dimpase

And I do not think that this should be positively reviewed yet: the error handling changes of Maxima in nintegral should be documented and the now useless lines removed.

I read in Nils' latest comment that there are no useless lines to remove, right? Regarding the documentation, how about adding
the following to the patch:

```diff --git a/sage/calculus/calculus.py b/sage/calculus/calculus.py
--- a/sage/calculus/calculus.py
+++ b/sage/calculus/calculus.py
@@ -632,7 +632,8 @@
- ``5`` - integral is probably divergent or slowly
convergent

-      - ``6`` - the input is invalid
+      - ``6`` - the input is invalid; this includes the case of desired_relative_error
+        being too small to be achieved

ALIAS: nintegrate is the same as nintegral
```

would this be sufficient?

### comment:64 follow-up: ↓ 65 Changed 4 years ago by jpflori

Im currently looking at Maxima's code and am not yet convinced of what changed in the last version.

There' indeed the possibility that an error get raised by Maxima (handler ... error stuff) if some values cannot be converted to floats float-or-lose stuff).

In this case Maxima returns the partially evaluated expression.

But that does not explain why Maxima now returns (0,0,0,6) rather than the old unevaluated expression for integrating x with a crazy expression. I'm currently looking at dqags.lisp but did not find anything meaningful yet.

In particular, it also makes me wonder whether the message of the Sage error (for the "workaround", but also for the earlier try/catch stuff) is meaningful.

Will it ever be related to precision problems ? because Maxima seems not to complain any more.

and as I said earlier, returning Error ourselves when when detect that the last value of the tuple and as Nils proposes to do, seems a good idea to me, but should be done later. I'd just like to actually reflect what Maxima tells at the Sage level even though the end user still has to check this last value for a moment.

### comment:65 in reply to: ↑ 64 Changed 4 years ago by nbruin

OK, this poking around in maxima's source is really quite unproductive, but I think I have located the source of the change in behaviour. With maxima 5.26, observe:

```sage: x.nintegrate(x,0,1,1e-14)
(0.0, 0.0, 0, 6)
sage: from sage.libs.ecl import *
sage: ecl_eval("(setf f2cl-lib::*stop-signals-error-p* T)")
<ECL: T>
sage: x.nintegrate(x,0,1,1e-14)
ValueError: Maxima (via quadpack) cannot compute the integral to that precision
```

whereas with maxima 5.23.2 we have the opposite:

```sage: x.nintegrate(x,0,1,1e-14)
ValueError: Maxima (via quadpack) cannot compute the integral to that precision
sage: from sage.libs.ecl import *
sage: ecl_eval("(setf f2cl-lib::*stop-signals-error-p* T)")
<ECL: T>
sage: x.nintegrate(x,0,1,1e-14)
(0.0, 0.0, 0, 6)
```

I think this commit is responsible]:

``` (defun stop (&optional arg)
(when arg
(format cl::*error-output* "~A~%" arg))
-  (unless *stop-signals-error-p*
+  (when *stop-signals-error-p*
(cerror "Continue anyway" "STOP reached")))
```

In later commits, the following is added:

```+;;; Revision 1.116  2011/02/20 20:51:04  rtoy
+;;; Oops.  STOP should signal an error if *STOP-SIGNALS-ERROR-P* is
+;;; non-NIL.
```

so I think the maxima 5.26 behaviour is as intended. I think you can give Dima's patch a positive review now. Everything is explained and any improvements should be handled on future tickets.

### comment:66 follow-up: ↓ 67 Changed 4 years ago by jpflori

Thanx a lot for locating this! It would have taken me ages to do so. But you did not quite answered one of my questions: we won't get the last error because of precision stuff, so this can get removed? My point is that those last lines were explicitely made to catch that specific error that is not reported as such anymore, so I think it's useless, or if it can then it won't be becase of precision provlems anymore so the stuff about precision should be removed (and the reference to the maxima mailing list should be removed as well).

I also agree it's not really important from a usefulness point of view, but having a good understanding of what we use in Sage is and our code should reflect that.

### comment:67 in reply to: ↑ 66 Changed 4 years ago by nbruin

we won't get the last error because of precision stuff, so this can get removed? My point is that those last lines were explicitely made to catch that specific error that is not reported as such anymore, so I think it's useless, or if it can then it won't be becase of precision provlems anymore so the stuff about precision should be removed

No, you already noted yourself that a failure to convert to float can trigger this:

```sage: var('n')
sage: x.nintegrate(x,0,n)
...
--> 712         raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
...
ValueError: Maxima (via quadpack) cannot compute the integral to that precision
```

that's why I said before that this check is necessary.

### comment:68 follow-up: ↓ 70 Changed 4 years ago by jpflori

Yup but that won't be because of too high precision and that's what our error message says...

So let's just remove the "to that precision" part.

### comment:69 Changed 4 years ago by jpflori

As soon as we agree on this last point, I'll be happy to provide a reviewer patch fixing some formatting problems (spaces at end of lines and stuff like that) and give the ticket a positive review.

To be clear about the domain stuff, I don't think we really need an additional "public" method to set it, fixing the maxima library interface so that sage.calculus.calculus.maxima.domain works as intended (ie can be read and assigned directly) seems enough for me.

And sorry about the typos above, it's late here...

### comment:70 in reply to: ↑ 68 ; follow-up: ↓ 71 Changed 4 years ago by nbruin

Yup but that won't be because of too high precision and that's what our error message says...

So let's just remove the "to that precision" part.

Yes, you're right. Agreed. In fact, I expect that any precision problems will always be reported via the error code, so they will *never* come out as an error or an unevaluated call. This only happened before because fortan's "STOP" raised an error, which it doesn't do now anymore.

That also makes me think that the body should be

```    v = ex._maxima_().quad_qags(x, a, b,
epsrel=desired_relative_error,
limit=maximum_num_subintervals)
#If an error gets raised during evaluation, maxima just returns the unevaluated expression
raise ValueError, "Maxima (via quadpack) cannot compute the integral"
return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
```

i.e., just let quad_qags raise its own errors instead of masking it with a blanket error message. Presently we are not aware of any error that would get raised under normal operation in there anyway, so we would want to know what went wrong.

### comment:71 in reply to: ↑ 70 ; follow-up: ↓ 72 Changed 4 years ago by dimpase

OK, I'll update the patch as follows:

```diff --git a/sage/calculus/calculus.py b/sage/calculus/calculus.py
--- a/sage/calculus/calculus.py
+++ b/sage/calculus/calculus.py
@@ -632,7 +632,8 @@
- ``5`` - integral is probably divergent or slowly
convergent

-      - ``6`` - the input is invalid
+      - ``6`` - the input is invalid; this includes the case of desired_relative_error
+        being too small to be achieved

ALIAS: nintegrate is the same as nintegral

@@ -641,9 +642,6 @@
integration using the GSL C library. It is potentially much faster
and applies to arbitrary user defined functions.

-    Also, there are limits to the precision to which Maxima can compute
-    the integral: here the error code "6" means "the input is invalid".
-
::

sage: f = x
@@ -714,14 +712,14 @@
limit=maximum_num_subintervals)
except TypeError, err:
if "ERROR" in str(err):
-            raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
+            raise ValueError, "Maxima (via quadpack) cannot compute the integral"
else:
raise TypeError, err

#This is just a work around until there is a response to
#http://www.math.utexas.edu/pipermail/maxima/2008/012975.html
-        raise ValueError, "Maxima (via quadpack) cannot compute the integral to that precision"
+        raise ValueError, "Maxima (via quadpack) cannot compute the integral"

return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
```

### comment:72 in reply to: ↑ 71 ; follow-up: ↓ 73 Changed 4 years ago by nbruin

In case anyone is waiting for my response:

• Dima's proposed patch update is fine with me
• I can't find Dima's updated patch on the ticket: Dima, in case you forgot to actually upload the amended patch, please do!
• After that, JP can fix some of the formatting that bothered him and then, for the first time since ages, we can ship an up-to-date maxima with the new sage release!

### comment:73 in reply to: ↑ 72 Changed 4 years ago by dimpase

• Status changed from needs_info to needs_review

In case anyone is waiting for my response:

• Dima's proposed patch update is fine with me
• I can't find Dima's updated patch on the ticket: Dima, in case you forgot to actually upload the amended patch, please do!

I haven't forgotten, I was just waiting for OK...
It's updated now.

• After that, JP can fix some of the formatting that bothered him and then, for the first time since ages, we can ship an up-to-date maxima with the new sage release!

### comment:74 Changed 4 years ago by jpflori

• Status changed from needs_review to needs_info

I'm in the process of fixing some formatting issues (from my point of view) in Sage source tree and in the SPKG itself.

Two remarks:

• I bumped the version number to 5.26.0.p0 because we use a patch on top of 5.26.0
• there is no spk-check although maxima provides a test suite.

I quickly looked for a good reason why there is no such file and could not yet locate one.

Therefore, I'm currently testing "make check" and I got one failure in rtest8 related to some quad_qawf not getting evaluated (this smells like an integration issue like above).

This gives me 1 test failed out of 9,121 in 360 secs (and it does not return an error at the console level...)

Can someone reproduce that failure with SPKG_CHECK=yes and the updated spkg at http://perso.telecom-paristech.fr/~flori/maxima-5.26.0.p0.spkg

### comment:75 follow-up: ↓ 77 Changed 4 years ago by jpflori

Sorry, there is a "sage" missing between flori and maxima...

And this issue is already reported here

although I'm not on OpenBSD/i386 but on a Ubuntu x86_64 system which surely is less important that the underlying lisp.

### comment:76 Changed 4 years ago by jpflori

• Reviewers changed from Jean-Pierre Flori, to Jean-Pierre Flori, Nils Bruin
• Work issues changed from domain issue, error handling in nintegral to spkg-check

Here is the address on the byg tracker:

### comment:77 in reply to: ↑ 75 Changed 4 years ago by kcrisman

although I'm not on OpenBSD/i386 but on a Ubuntu x86_64 system which surely is less important that the underlying lisp.

Correct; the Maxima list is always full of minor failures in their (very large, I believe) test suite. This shouldn't necessarily hold up positive review, because it is quite dependent on the Lisp, if I recall from reading those messages. Also, adding the test suite should in theory be another ticket; no need to hold up having an actual up-to-date Maxima!

### comment:78 follow-up: ↓ 79 Changed 4 years ago by jpflori

• Status changed from needs_info to needs_review
• Work issues spkg-check deleted

Ok, I'll split the addition of the spkg-check.

Nils, could you have a look at my "reviewer" patch (that I'll upload in a minute for the sage tree, you can already look at the mercurial log for the spkg, ignoring the last commit for the spkg-check).

I got rid of a lot of "useless" spaces, jsut want to be sure every one is happy with that.

### Changed 4 years ago by jpflori

Reviewer patch, formatting issues

### Changed 4 years ago by jpflori

spkg diff from author, for review only

### comment:79 in reply to: ↑ 78 Changed 4 years ago by nbruin

• Status changed from needs_review to positive_review

I got rid of a lot of "useless" spaces, jsut want to be sure every one is happy with that.

Appreciated! I happen to believe that while it's nice to fix whitespace, increasing the footprint of patches by fixing it in lines that are otherwise not affected is not worth it, so I'm not too happy about it. But I don't feel sufficiently strongly about it to enter into a discussion about whether you should put the whitespace back. So I won't hold up a positive review on this. You can make up your own mind on whether you want to fix whitespace like this in the future.

### comment:80 follow-up: ↓ 82 Changed 4 years ago by jpflori

• Description modified (diff)

I usually do not do it, especially so many lines at a time, but my reviewer patch was all about removing some white spaces reintroduced by Dima, and split some overlongish lines (that's another tricky choice, I can leave with long lines, but it's much more readable in a terminal without them and according to our coding guidelines we should stick to 79 chars at max...).

So I chose to suppress them this time.

IIRC there was a discussion some time ago on sae-devel about make a one time huge patch for suppressing spaces but some people, especially among sage-combinat, raised the problem that it would break patch queues and that we should rather ensure that no such new lines where added and fix for previous spaces made locally.

I hope my choice won't break anything.

Anyway, good that this ticket is finally positively reviewed.

### comment:81 Changed 4 years ago by jdemeyer

• Description modified (diff)
• Summary changed from sum fails with lower bound != 0 or 1 (upgrade maxima to 5.26) to Upgrade maxima to 5.26

### comment:82 in reply to: ↑ 80 ; follow-up: ↓ 83 Changed 4 years ago by jdemeyer

• Status changed from positive_review to needs_work

IIRC there was a discussion some time ago on sae-devel about make a one time huge patch for suppressing spaces but some people, especially among sage-combinat, raised the problem that it would break patch queues and that we should rather ensure that no such new lines where added and fix for previous spaces made locally.

For precisely this reason I would prefer not to just remove all spaces everywhere. In code which is changed by the patch, it can be done. But not just everywhere like this.

### comment:83 in reply to: ↑ 82 ; follow-up: ↓ 84 Changed 4 years ago by dimpase

• Description modified (diff)
• Status changed from needs_work to needs_review

IIRC there was a discussion some time ago on sae-devel about make a one time huge patch for suppressing spaces but some people, especially among sage-combinat, raised the problem that it would break patch queues and that we should rather ensure that no such new lines where added and fix for previous spaces made locally.

For precisely this reason I would prefer not to just remove all spaces everywhere. In code which is changed by the patch, it can be done. But not just everywhere like this.

would it be OK if we just not apply the reviewer's patch?

### comment:84 in reply to: ↑ 83 ; follow-up: ↓ 85 Changed 4 years ago by jpflori

• Status changed from needs_review to needs_work

IIRC there was a discussion some time ago on sae-devel about make a one time huge patch for suppressing spaces but some people, especially among sage-combinat, raised the problem that it would break patch queues and that we should rather ensure that no such new lines where added and fix for previous spaces made locally.

For precisely this reason I would prefer not to just remove all spaces everywhere. In code which is changed by the patch, it can be done. But not just everywhere like this.

would it be OK if we just not apply the reviewer's patch?

I don't think so. My patch is overkill, but it's no reason to reintroduce additional formatting issues inside Sage code. If we go on like this, its going to be a real mess.

### comment:85 in reply to: ↑ 84 Changed 4 years ago by dimpase

• Status changed from needs_work to needs_info

IIRC there was a discussion some time ago on sae-devel about make a one time huge patch for suppressing spaces but some people, especially among sage-combinat, raised the problem that it would break patch queues and that we should rather ensure that no such new lines where added and fix for previous spaces made locally.

For precisely this reason I would prefer not to just remove all spaces everywhere. In code which is changed by the patch, it can be done. But not just everywhere like this.

would it be OK if we just not apply the reviewer's patch?

I don't think so. My patch is overkill, but it's no reason to reintroduce additional formatting issues inside Sage code. If we go on like this, its going to be a real mess.

I haven't noticed that I have introduced any formatting issues. Could you point me out to this, or even better, of course, create a patch that fixes my faults?
Thanks!

### comment:86 follow-up: ↓ 90 Changed 4 years ago by jpflori

I'll provide a patch asap of course, that's basically trailing whitespaces (e.g. last lines of your patch) in your lines and too long lines (e.g. first line of your patch) and wrong indentation for some doc (e.g. first line of your patch).

It's just that the current patch do not apply on 4.8 and as compiling atlas on my computer takes ages and installing properly atlas on my debian seems broken...

More seriously, you also removed the lines

"Also, there are limits to the precision to which Maxima can compute..."

I agree that we do not raise an error anymore, but Maxima can still not compute up to that precision.

You also left the stuff about waiting for an answer on Maxima's list, but IIRC there has been an answer there, so what is done is not really a workaround, but rather a decision from Sage developpers.

### comment:87 follow-ups: ↓ 89 ↓ 91 Changed 4 years ago by jpflori

Last questions (I hope):

About the lines patched in wester.py, why did you add a semi-colon after "n=var('n')" ?

And what do you mean exactly by "Maxima 5.26 does not do "just" simpify() here." (sic)?

That now Maxima cannot simplify this expression when it is only told that Re(x) >0 and Re(y)>0 and that we have to call simplify_exp() instead ?

I just check that this is indeed the case now and was not with Maxima's 5.24.0.p0 (ok this could be simply deduced from the doctest changes).

Is this a decision from Maxima's developers or a bug introduced recently? if the former holds a pointer to where the change was made would be welcome; if the latter holds, this should be reported upstream.

I guess that this is related to the domain stuff, because setting domain back to real answer right away.

And it's strange that the assumptions is not needed anymore.

In fact looking into the simplify_exp code shows that maxima.eval('domain: real\$') is called before calling the radcan method which actually simplifies the expression.

In particular, this might actually be a Maxima regression.

### comment:88 Changed 4 years ago by jpflori

• Description modified (diff)
• Status changed from needs_info to needs_review

### comment:89 in reply to: ↑ 87 Changed 4 years ago by nbruin

Thanks for updating the patch! It's OK with me. I'll leave the specific questions for Dima to answer.

That now Maxima cannot simplify this expression when it is only told that Re(x) >0 and Re(y)>0 and that we have to call simplify_exp() instead ?
In particular, this might actually be a Maxima regression.

I agree that this is a Maxima regression. First maxima could verify an equality given sufficient information, and now it can't. Assuming everything in sight is real is not the same as assuming some real parts are positive (although who knows what awful bug had as a side-effect this particular equality got verified).

If you look in the history of simp.lisp you'll see there is some recent probably relevant work there, so I think reporting this issue might have a chance of getting it fixed (and will certainly be appreciated).

We are at the mercy of Maxima for these issues. Our alternatives are:

• not upgrade and live with existing bugs
• live with the regression
• fix and patch maxima

I think in this case, Option 2 is the practical one. Option 3 is obviously the best (both for Maxima and for Sage) but doesn't solve anything _now_.

As soon as Dima has commented, back to positive review!

### comment:90 in reply to: ↑ 86 Changed 4 years ago by dimpase

More seriously, you also removed the lines

"Also, there are limits to the precision to which Maxima can compute..."

I agree that we do not raise an error anymore, but Maxima can still not compute up to that precision.

these lines essentially duplicated the description of return code 6 several lines above these.

You also left the stuff about waiting for an answer on Maxima's list, but IIRC there has been an answer there, so what is done is not really a workaround, but rather a decision from Sage developpers.

well, that's the stuff I had no knowledge about, so I left it stand, obviously.

Dima

### comment:91 in reply to: ↑ 87 Changed 4 years ago by dimpase

Last questions (I hope):

About the lines patched in wester.py, why did you add a semi-colon after "n=var('n')" ?

hmm, perhaps a result of a copy/paste from Sage prompt...

And what do you mean exactly by "Maxima 5.26 does not do "just" simpify() here." (sic)?

That now Maxima cannot simplify this expression when it is only told that Re(x) >0 and Re(y)>0 and that we have to call simplify_exp() instead ?

Indeed.

### comment:92 follow-up: ↓ 93 Changed 4 years ago by jpflori

Ok, I understand. Nonetheless I think its better to be more verbose as was the case before, especially that no error is now explicitely returned.

If someone does not pay attention, one might think that Sage returns 0 for the integral of x between 0 an 1.

### comment:93 in reply to: ↑ 92 Changed 4 years ago by dimpase

Ok, I understand. Nonetheless I think its better to be more verbose as was the case before, especially that no error is now explicitely returned.

If someone does not pay attention, one might think that Sage returns 0 for the integral of x between 0 an 1.

OK, I updated the patch as follows, to reflect this and a couple of other minor points raised:

```diff --git a/sage/calculus/calculus.py b/sage/calculus/calculus.py
--- a/sage/calculus/calculus.py
+++ b/sage/calculus/calculus.py
@@ -641,6 +641,8 @@
``numerical_integral`` that implements numerical
integration using the GSL C library. It is potentially much faster
and applies to arbitrary user defined functions.
+    Also, there are limits to the precision to which Maxima can compute
+    the integral to due to limitations in quadpack.

::

diff --git a/sage/calculus/wester.py b/sage/calculus/wester.py
--- a/sage/calculus/wester.py
+++ b/sage/calculus/wester.py
@@ -383,9 +383,10 @@
::

sage: # (YES) Assuming Re(x)>0, Re(y)>0, deduce x^(1/n)*y^(1/n)-(x*y)^(1/n)=0.
-    sage: # Maxima 5.26 does not do "just" simpify() here. Thus simplify_exp() used.
+    sage: # Maxima 5.26 cannot do "just" simpify() here. Thus simplify_exp() is used.
+    sage: # This is a regression from 5.24
sage: # assume(real(x) > 0, real(y) > 0) # (not needed for simplify_exp())
-    sage: n = var('n');
+    sage: n = var('n')
sage: f = x^(1/n)*y^(1/n)-(x*y)^(1/n)
sage: f.simplify_exp()
0
```

### comment:94 follow-up: ↓ 95 Changed 4 years ago by jpflori

In fact, I already corrected these issues in my patch, so if you're happy with my corrections, you'd better reupload the previous version of your patch and let the ticket have positive review.

Otherwise, the reviewer patch will have to be rebased.

### Changed 4 years ago by dimpase

Sage library patch to go with the new spkg

### comment:95 in reply to: ↑ 94 Changed 4 years ago by dimpase

In fact, I already corrected these issues in my patch, so if you're happy with my corrections, you'd better reupload the previous version of your patch and let the ticket have positive review.

Used Apple's Time Machine backup fot the 1st time ever :-)

### comment:96 Changed 4 years ago by kcrisman

```"just" simpify()
```

I assume that, since you are changing that line anyway, it would be trivial to edit "by hand" the text .patch file so that this reads

```"just" simplify()
```

It seems to be the sort of thing we should fix :)

### comment:97 follow-up: ↓ 99 Changed 4 years ago by jpflori

These lines are in fact changed in the reviewer patch.

### comment:98 Changed 4 years ago by kcrisman

Sorry, I just saw what you pasted above. Sorry for the noise.

### comment:99 in reply to: ↑ 97 ; follow-up: ↓ 100 Changed 4 years ago by nbruin

• Status changed from needs_review to needs_work

I'm sorry to report that the new spkg seems to cause problems for the expect interface for maxima. I don't think this problem arose with Dima's original spkg. Of course, I'd be very happy to hear that other people cannot reproduce this problem and that I'm not building ecl/maxima correctly, but currently I'm experiencing:

```sage: maxima(1)
[HANGS]
```

This happens to me on two different builds where I've applied the recipe as outlined in the ticket description. For one, see sage.math.
As it stands, I don't think this ticket will pass doctests (it doesn't for me).

### comment:100 in reply to: ↑ 99 Changed 4 years ago by dimpase

I'm sorry to report that the new spkg seems to cause problems for the expect interface for maxima. I don't think this problem arose with Dima's original spkg. Of course, I'd be very happy to hear that other people cannot reproduce this problem and that I'm not building ecl/maxima correctly, but currently I'm experiencing:

```sage: maxima(1)
[HANGS]
```

This happens to me on two different builds where I've applied the recipe as outlined in the ticket description. For one, see sage.math.
As it stands, I don't think this ticket will pass doctests (it doesn't for me).

I also see weirdness: e.g. some doctests ran out of time with the present spkg/patches pair, while they worked fine with my original spkg/patches pair...

### comment:101 Changed 4 years ago by jpflori

• Status changed from needs_work to needs_info

It's strange that you report that it did not happen with Dima original package because the only changes I've made to the spkg are purely typographical (see trac_10682-spkg.patch Download).
I'll recheck on my computer.

### comment:102 Changed 4 years ago by jpflori

I confirm that something went wrong with the last spkg I uploaded, maybe I downloaded a wrong base spkg, I'll check that now.

### Changed 4 years ago by jpflori

spkg diff from reviewer, for review only

### comment:103 Changed 4 years ago by jpflori

Ok that was an error from my part, some TABS were actually needed in spkg-install.
This is fixed now.
I've also swpped two lines in the sage tree patch.
I'll upload a fixed version in a moment.

Reviewer patch

### comment:104 follow-up: ↓ 105 Changed 4 years ago by jpflori

• Status changed from needs_info to needs_review

make ptest is fine on my computer now (ubuntu/amd64).

### comment:105 in reply to: ↑ 104 Changed 4 years ago by dimpase

make ptest is fine on my computer now (ubuntu/amd64).

make ptestlong went well on MacOSX 10.6.8.
looks like it's good to go!

### comment:106 Changed 4 years ago by nbruin

• Status changed from needs_review to positive_review

no complaints here.

### comment:107 Changed 4 years ago by jdemeyer

• Merged in set to sage-5.0.beta8
• Resolution set to fixed
• Status changed from positive_review to closed

### comment:108 Changed 4 years ago by jdemeyer

• Authors changed from Dima Pasechnik to Dmitrii Pasechnik
Note: See TracTickets for help on using tickets.