Opened 2 years ago

Last modified 5 months ago

#31411 new defect

Expand isn't absorbing (Iterated expand may fail).

Reported by: Emmanuel Charpentier Owned by:
Priority: major Milestone: sage-9.8
Component: symbolics Keywords: symbolics expand
Cc: Karl-Dieter Crisman Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description

Inspired by this question on ask.sagemath.org :

    sage: var("q, jt")
    (q, jt)
    sage: A=q^(2/3)+q**(2/5)
    sage: B=product(1 -q**jt, jt, 1 , 31)*q**(1 /24 )
    sage: bool((A*B).expand()==(A*B.expand()).expand())
    False
    sage: (((A*B).expand())/(A*B.expand()).expand()).factor()
    2*q^(1/24)/(q^(17/24) + q^(53/120))

However :

    sage: bool(A*B==A*B.expand()) # Damn slow...
    True
    sage: ((A*B)/(A*B.expand())).factor()
    1

The iteration of expand should not change the value of the expanded quantity.

Severity set to blocker because this can give mathematically incorrect results on basic computations met in everyday situations.

Change History (12)

comment:1 Changed 2 years ago by Emmanuel Charpentier

This seems not to be a Maxima bug :

;;; Loading #P"/usr/local/sage-9/local/lib/ecl/sb-bsd-sockets.fas"
;;; Loading #P"/usr/local/sage-9/local/lib/ecl/sockets.fas"
Maxima 5.44.0 http://maxima.sourceforge.net
using Lisp ECL 20.4.24
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) A:q^(2/3)+q**(2/5);

(%o2) q^(2/3)+q^(2/5)
(%i3) B:product(1 -q**jt, jt, 1 , 31)*q**(1 /24 );

(%o3) (1-q)*q^(1/24)*(1-q^2)*(1-q^3)*(1-q^4)*(1-q^5)*(1-q^6)*(1-q^7)*(1-q^8)
           *(1-q^9)*(1-q^10)*(1-q^11)*(1-q^12)*(1-q^13)*(1-q^14)*(1-q^15)
           *(1-q^16)*(1-q^17)*(1-q^18)*(1-q^19)*(1-q^20)*(1-q^21)*(1-q^22)
           *(1-q^23)*(1-q^24)*(1-q^25)*(1-q^26)*(1-q^27)*(1-q^28)*(1-q^29)
           *(1-q^30)*(1-q^31)
(%i4) is(equal(expand(A*B),expand(A*expand(B))));

(%o4) true

Pynac ?

comment:2 Changed 2 years ago by Emmanuel Charpentier

Sympy seems exempt of the problem :

>>> from sympy import symbols, Rational
>>> from math import prod
>>> A = q**Rational(2,3)+q**Rational(2,5)
>>> B = prod([1-q**(u+1) for u in range(31)])*q**Rational(1,24)
>>> ((A*B).expand()==(A*B.expand()).expand())
True

comment:3 Changed 2 years ago by Martin Rubey

I wouldn't be surprised of a connection with #31077 (which should be a blocker, too, in my opinion).

comment:4 Changed 2 years ago by Dima Pasechnik

Priority: blockercritical

this is an oldish pynac bug, same happens with almost 3 years old https://github.com/pynac/pynac/releases/tag/pynac-0.7.19

thus, not a blocker, just crucial, I think.

comment:5 Changed 2 years ago by Karl-Dieter Crisman

Cc: Karl-Dieter Crisman added

comment:6 Changed 2 years ago by Dave Morris

The patch at #31077 seems to solve this problem.

comment:7 Changed 2 years ago by Dave Morris

This is not the same issue as #31077 (because there are no negative exponents here). The fact that the patch there fixes this ticket should mean that the problem is related to singular (or the interface between pynac and singular).

comment:8 Changed 23 months ago by Dave Morris

Milestone: sage-9.3sage-9.4
Priority: criticalmajor

There are two problems.

The first problem is that if a variable appears in a sum with two (or more) different fractional exponents, such as x^(1/2) + x^(1/3), then pynac's poly_mul_expand function finds the gcd of the exponents (call it d), and introduces a new variable (call it z) to represent x^d. Then the original terms can be represented as integer powers of the new variable. In our example they are z^3 and z^2. The first problem is that pynac calculates the gcd incorrectly if one of the exponents is 1. I have opened ticket #31477 to fix this.


The other problem shows up when a variable arises with a fractional exponent, but also arises with no exponent at all, such as in x^(1/2) + x^(1/3) + x. In this case, poly_mul_expand erroneously substitutes x^d for the term that has no exponent (where d is the gcd that was calculated above). This is is clearly wrong, because x will almost never be equal to x^d.

The code should be rewritten so that a term with no exponent is treated as if it has 1 as an exponent. This will not be difficult, but it will require some refactoring to avoid excessive code duplication.

I have opened metaticket #31478 to track this and other problems with poly_mul_expand. Given that the code has these two bugs, plus the bug in #31077, I think it needs a thorough examination, instead of just patches for the bugs that we have noticed. Ticket #31479 disables the use of poly_mul_expand, so this is not an urgent problem.

comment:9 Changed 19 months ago by Matthias Köppe

Milestone: sage-9.4sage-9.5

comment:10 Changed 13 months ago by Matthias Köppe

Milestone: sage-9.5sage-9.6

comment:11 Changed 10 months ago by Matthias Köppe

Milestone: sage-9.6sage-9.7

comment:12 Changed 5 months ago by Matthias Köppe

Milestone: sage-9.7sage-9.8
Note: See TracTickets for help on using tickets.