Fix sign of symbolic product product(1 - q^k, k, 1, N)

[sagetrac-tmonteil]

Reported at Ask Sage question 53345.

The output for this product
has an incorrect leading minus sign:

sage: N, q, k = SR.var('N, q, k')
sage: product(1 - q^k, k, 1, N)
-(-1)^N*product(q^k - 1, k, 1, N)

Expected:

sage: N, q, k = SR.var('N, q, k')
sage: product(1 - q^k, k, 1, N)
(-1)^N*product(q^k - 1, k, 1, N)

Fixed in Maxima by this commit:

• https://sourceforge.net/p/maxima/code/ci/762fd858d2a25fa1b8b6f5423909adbc365b3c5a/

CC: @DaveWitteMorris @macrakis @robert-dodier @kcrisman @nbruin @slel

Component: symbolics

Keywords: maxima

Author: Dave Morris

Branch/Commit: 69c0235

Reviewer: Samuel Lelièvre

Issue created by migration from https://trac.sagemath.org/ticket/30520

[EmmanuelCharpentier]

comment:1

FWIW, Sympy doesn't seem to be able to pull this one:

sage: from sympy import Product                                                 
sage: from sympy.abc import j, p, q                                             
sage: foo=Product(1-q^j, [j,1,p]); foo                                          
Product(1 - q**j, (j, 1, p))
sage: type(foo)                                                                 
<class 'sympy.concrete.products.Product'>
sage: foo.doit()                                                                
Product(1 - q**j, (j, 1, p))

Mathematica isn't more helpful:

sage: bar=mathematica("Product[1-q^j, {j, 1, p}]"); bar                         
QPochhammer[q, q, p]

Giac's answer can't be believed :

age: var("j,q,p")                                                              
(j, q, p)
sage: gargs=list(map(giac, [1-q^j,j,1,p,1]))                                    
sage: gargs                                                                     
[-q^j+1, j, 1, p, 1]
sage: libgiac.product(*gargs)                                                   
-q+1

HTH, but doubting it...

[heluani]

comment:2

for what it's worth, it's either in the interface to maxima in maxima_eval or in maxima itself, what ultimately gets called is:

sage: var('n,q,i')
(n, q, i)
sage: expr = 1 - q**i
sage: from sage.interfaces import maxima_lib
sage: sage.interfaces.maxima_lib.maxima.sr_prod(expr, i, 1, n)
-(-1)^n*product(q^i - 1, i, 1, n)

I couldn't reproduce this on pure maxima, but I have never used it before, so it may well be there.

EDIT: it may also be in the final conversion max_to_sr

[fchapoton]

comment:3

A simple case

sage: from sage.calculus.calculus import symbolic_product                       
sage: var('x,n')                                                                
(x, n)
sage:  symbolic_product(-x**2,x,1,n)                                            
-(-1)^n*factorial(n)^2

[fchapoton]

comment:4

This seems to happen in max_simplify_prod. Before:

ipdb> p maxima_eval(([max_prod],[sr_to_max(SR(a)) for a in args]))              
<ECL: ((%PRODUCT SIMP) ((MTIMES SIMP) -1 ((MEXPT SIMP) |$_SAGE_VAR_x| 2))
 |$_SAGE_VAR_x| 1 |$_SAGE_VAR_n|)>

and after:

ipdb> p maxima_eval([max_simplify_prod, ([max_prod],[sr_to_max(SR(a)) for a in a
rgs])])                                                                         
<ECL: ((MTIMES SIMP) ((MEXPT SIMP) -1 ((MPLUS SIMP) -1 |$_SAGE_VAR_n|))
 ((MEXPT SIMP) ((MFACTORIAL SIMP) |$_SAGE_VAR_n|) 2))>

[slel]

comment:6

Any insights from Maxima or Sage symbolics experts?

[macrakis]

comment:7

Hi, thanks for the Maxima problem report.

I was unable to reproduce this in Maxima 5.44.0/SBCL 2.0.0/Windows 10.

Could you please report the details necessary to reproduce the problem:

• The exact version of Maxima/SBCL/OS you're using.

• The exact commands that were passed to Maxima.

• Any non-default global settings that Sage uses when calling Maxima.

• What does max_simplify_product do?

• What the var(...) declaration/statement does in Maxima terms. Does it declare anything special about those variables?

• What does "before this ticket" and "after this ticket" mean? Was it run on different versions of Maxima? If so, which ones?
  Thanks,

        -s
  

[slel]

Changed keywords from none to maxima

[slel]

comment:8

Apologies for cc-ing you. Rather than in Maxima,
the bug might well be in Sage's interface to it.

In Sage, the bug can lead to wrong results as in

sage: from sage.calculus.calculus import symbolic_product                       
sage: x, n = SR.var('x,n')
sage: symbolic_product(-x**2, x, 1, n)
-(-1)^n*factorial(n)^2

or cause an error to be raised as in

sage: symbolic_product(-x, x, 1, n)
Traceback (most recent call last)
...
RuntimeError: ECL says: Error executing code in Maxima: factorial: factorial of negative integer -1 not defined.

Not sure what exactly gets sent to Maxima in those cases.

I have Maxima 5.44.0, ECL 21.2.1, macOS 10.14.6.

sage: pv = installed_packages()  # package versions
sage: print(', '.join(f'{p} {pv[p]}' for p in ['maxima', 'ecl']))
maxima 5.44.0, ecl 21.2.1

Sage calls Maxima with the following global settings:

besselexpand : true ;
display2d : false ;
domain : complex ;
keepfloat : true ;
load(to_poly_solve) ;
load(simplify_sum) ;
load(diag) ;
nolabels : true ;

"Before this ticket" is the observed bug.
"After this ticket" is the goal (so far not reached).

[kcrisman]

comment:9

• What does max_simplify_product do?

Here:

max_simplify_prod=EclObject("$SIMPLIFY_PRODUCT")

which is presumably in the contributed simplify_sum package, more precisely here. I had trouble finding references in the official documentation, though.

• What the var(...) declaration/statement does in Maxima terms. Does it declare anything special about those variables?

No.

[robert-dodier]

comment:10

"Before this ticket" is the observed bug. "After this ticket" is the goal (so far not reached). I know this is just a convention, but, um, this is pretty confusing; I don't see how anyone who was not already among the cognoscenti would know what it means. "Observed" and "expected behavior" seems adequate.

I can't be sure about what's going on here, but I speculate the error is coming out of simplify_product in the simplify_sum code.

simplify_product ('product(1 - f(k), k, 1, N));
 => (-1)^(N-1)*'product(f(k)-1,k,1,N)

The leading term is off by 1 (right?) -- I think this is the origin of the stray -1 which was observed.

[slel]

comment:11

Thanks for the rephrasing suggestion and your thoughts on the problem.

[DaveWitteMorris]

comment:12

@robert-dodier: Yes, I think you are right. In line 1112 of the file suggested by kcrisman, the code for simplify_product says:

(-1)^(hi-lo)*simplify_product(apply(nounify('product), [-term, %n, lo, hi]))

If I understand correctly, then this is an off-by-one error: the number of terms in the product is hi - lo + 1 because hi and lo are both included.

[slel]

comment:13

I located the copy of that file in my Sage installation:

$ cd $(sage -c 'print(SAGE_ROOT)')
$ find . -name "solve_rec.mac"
./local/share/maxima/5.44.0/share/solve_rec/solve_rec.mac

and changed line 1112 as suggested by Dave.

Now I get the correct product from the ticket description:

sage: N, q, k = SR.var('N, q, k')
sage: product(1 - q^k, k, 1, N)
(-1)^N*product(q^k - 1, k, 1, N)

and for the sum of -x^2:

sage: from sage.calculus.calculus import symbolic_product
sage: x, n = SR.var('x, n')
sage: symbolic_product(-x**2, x, 1, n)
(-1)^n*factorial(n)^2

but this still gives trouble:

sage: from sage.calculus.calculus import symbolic_product
sage: x, n = SR.var('x, n')
sage: symbolic_product(-x, x, 1, n)
#0: simplify_product(product='product(-_SAGE_VAR_x,_SAGE_VAR_x,1,_SAGE_VAR_n))
...
Traceback (most recent call last)
...
RuntimeError: ECL says: Error executing code in Maxima: factorial: factorial of negative integer -1 not defined.

It would be nice to fix that too.

[robert-dodier]

comment:14

I've come to the same conclusion. I think this patch fixes it:

$ git diff
diff --git a/share/solve_rec/solve_rec.mac b/share/solve_rec/solve_rec.mac
index 2727db808..8979e3f1d 100644
--- a/share/solve_rec/solve_rec.mac
+++ b/share/solve_rec/solve_rec.mac
@@ -1109,7 +1109,7 @@ simplify_product(prod) := block(
         prod
     )
     else if part(term, 0)="-" and length(args(term))=1 then
-      (-1)^(hi-lo)*simplify_product(apply(nounify('product), [-term, %n, lo, hi]))
+      (-1)^(hi-lo+1)*simplify_product(apply(nounify('product), [-term, %n, lo, hi]))
     /* Give up! */
     else
       p%*apply(nounify('product), [factor(term), %n, lo, hi])

Although the existing, incorrect code is exercised by the unit tests (rtest_simplify_sum and rtest_solve_rec), fixing it does not change any results! Not sure what's going on there; I guess it is a happy accident. Or not.

Anyway I will add a few unit tests and commit the patch.

[robert-dodier]

comment:15

Replying to @slel:

but this still gives trouble:

sage: from sage.calculus.calculus import symbolic_product
sage: x, n = SR.var('x, n')
sage: symbolic_product(-x, x, 1, n)
#0: simplify_product(product='product(-_SAGE_VAR_x,_SAGE_VAR_x,1,_SAGE_VAR_n))
...
Traceback (most recent call last)
...
RuntimeError: ECL says: Error executing code in Maxima: factorial: factorial of negative integer -1 not defined.

It would be nice to fix that too.

Well, that's a different bug, although it's in the same file; looks like it is in the stanza which begins if freeof(%n,aa) and freeof(%n,bb) around line 1096. In the interest of clarity, I would suggest opening a separate bug report about it.

[slel]

comment:16

Replying to @robert-dodier:

Replying to @slel:

but this still gives trouble:

sage: from sage.calculus.calculus import symbolic_product
sage: x, n = SR.var('x, n')
sage: symbolic_product(-x, x, 1, n)
#0: simplify_product(product='product(-_SAGE_VAR_x,_SAGE_VAR_x,1,_SAGE_VAR_n))
...
Traceback (most recent call last)
...
RuntimeError: ECL says: Error executing code in Maxima: factorial: factorial of negative integer -1 not defined.

It would be nice to fix that too.

Well, that's a different bug, although it's in the same file; looks like it is in the stanza which begins if freeof(%n,aa) and freeof(%n,bb) around line 1096. In the interest of clarity, I would suggest opening a separate bug report about it.

You are right. This is now tracked at #31557.

[robert-dodier]

comment:17

Fixed (to the best of my knowledge; no doubt you'll want to verify for Sage) by commit 762fd85 on maxima-code/master which applies the patch for simplify_product shown above.

[DaveWitteMorris]

comment:18

Thanks!

[slel]

comment:19

Should we apply the upstream commit fixing this
as a patch or wait for the next Maxima release?

[robert-dodier]

comment:20

Replying to @slel:

Should we apply the upstream commit fixing this
as a patch or wait for the next Maxima release?

For what it's worth, a new Maxima version is on the horizon; it is being discussed on the mailing list, although no work on it has been announced. So my guess is that a new version will appear in, let's say, one to three months. I don't know whether that implies you should wait or not.

[mkoeppe]

comment:21

Moving to 9.4, as 9.3 has been released.

[DaveWitteMorris]

Branch: public/30520

[DaveWitteMorris]

Commit: 69c0235

[DaveWitteMorris]

Author: Dave Morris

[DaveWitteMorris]

comment:23

This was fixed in maxima 5.45. Ticket #31876 did the upgrade in sage. The PR just adds a doctest.

[slel]

Reviewer: Samuel Lelièvre

[slel]

comment:24

Lovely.

[vbraun]

Changed branch from public/30520 to 69c0235