Opened 10 months ago

Closed 8 months ago

#31137 closed defect (fixed)

Fix subs failure involving integer-valued rational exponents

Reported by: charpent Owned by:
Priority: minor Milestone: sage-9.3
Component: symbolics Keywords: factor
Cc: slelievre Merged in:
Authors: Dave Morris Reviewers: Dima Pasechnik
Report Upstream: N/A Work issues:
Branch: c2b45c0 (Commits, GitHub, GitLab) Commit: c2b45c0b3181ff8803e116e4c92ecf4c4bd748d9
Dependencies: #30786, #30446 Stopgaps:

Status badges

Description (last modified by slelievre)

From Ask Sage question 54986: subs fails with a memory error (or doesn't return) on some expressions.

The following works fine:

sage: _ = var('A, L, G, R, f, k, n, q, u, beta, gamma', domain="positive")
sage: a = I*R^2*f^3*k*q*A*u
sage: b = 2*pi*L*R^2*G*f^4*k^2*q - 2*pi*L*R^2*G*f^4*q - 2*pi*L*R^2*beta^2*G*q
sage: c = (2*I*pi*L*R^2*beta*gamma*q + 2*I*pi*L*R*(beta + q))*G*f^3
sage: d = 2*(pi*(beta^2 + 1)*L*R^2*q + pi*L*R*beta*gamma*q + pi*L*beta)*G*f^2
sage: e = (-2*I*pi*L*R^2*beta*gamma*q - 2*I*pi*(beta^2*q + beta)*L*R)*G*f
sage: expr = a / ((b + c + d + e)*n)
sage: R1 = (expr.real()^2 + expr.imag()^2).factor()
sage: R2 = (sqrt(expr.real()^2 + expr.imag()^2)^2).factor()
sage: R3 = ((sqrt(expr.real()^2 + expr.imag()^2).factor())^2).factor()
sage: all((R1 == R2, R1 == R3, R2 == R3))

These expressions have the same string representation:

sage: [len(str(ex)) for ex in (R1, R2, R3)]
[734, 734, 734]
sage: all((str(R1) == str(R2), str(R1) == str(R3), str(R2) == str(R3)))

Substituting f = 2*beta in R1 R2 R3 works:

sage: R1s = R1.subs(f = 2*beta)
sage: R2s = R2.subs(f = 2*beta)
sage: R3s = R3.subs(f = 2*beta)

However, the following never return:

sage: bool(R3s == R1s)  # hangs
sage: len(str(R1s))
sage: len(str(R2s))
sage: len(str(R3s))  # hangs

Attempting to interrupt the last call can result in a Sage crash; the original poster reports a "Memory Error" exception.

As a comparison, substituting in R3 via SymPy works:

sage: import sympy
sage: bool(sympy.sympify(R3).subs(f, 2*beta)._sage_() == R1.subs(f = 2*beta))

but seems to fail via Mathematica somehow:

sage: mma = mathematica
sage: bool(mma.Replace(R3, mma.Rule(f, 2*beta)).sage() == R1.subs(f = 2*beta))

Priority set to critical, because it affects a very basic feature of Sage.

RESOLUTION: solved by the patch to pynac in #30446. There was a bug in evaluating an expression (perhaps sqrt(expr.real())^2) that has a rational exponent (such as 2/2) that simplifies to an integer. This ticket just adds a doctest to ensure that the problem is not allowed to return.

Change History (13)

comment:1 Changed 10 months ago by slelievre

  • Cc slelievre added
  • Description modified (diff)

Seems to have a lot to do with factor.

I also got a crash when interrupting the computation:

Exception ignored in: 'sage.libs.pynac.pynac.py_repr'
Traceback (most recent call last):
  File "sage/rings/integer.pyx", line 1013, in sage.rings.integer.Integer.__repr__ (build/cythonized/sage/rings/integer.
  File "sage/rings/integer.pyx", line 1134, in sage.rings.integer.Integer.str (build/cythonized/sage/rings/integer.c:901
Unhandled SIGSEGV: A segmentation fault occurred.
This probably occurred because a *compiled* module has a bug
in it and is not properly wrapped with sig_on(), sig_off().
Python will now terminate.
Segmentation fault

comment:2 Changed 10 months ago by slelievre

  • Description modified (diff)
  • Keywords factor added

Seems expressions resulting from factor behave strangely. Not sure what role subs plays here.

comment:3 Changed 10 months ago by slelievre

What seemed to be hanging was just taking a long time. Still, the result is surprising.

sage: len(str(R3s))  # long time -- 1,292,914,505

comment:4 Changed 10 months ago by gh-DaveWitteMorris

  • Branch set to public/31137

comment:5 Changed 10 months ago by gh-DaveWitteMorris

  • Authors set to Dave Morris
  • Commit set to 5da750a7fca49b55bd131585f0c620490611834b
  • Dependencies set to #30786
  • Priority changed from critical to minor
  • Status changed from new to needs_review

The problem in this ticket seems to be solved by the patch to pynac in #30446. The PR just adds a doctest. It is based on #30786 in order to avoid a merge conflict.

New commits:

dc073dffixes for trac #30446
4a5d053doctests for trac 28620, 30304, 30786
5da750adoctest for trac 31137

comment:6 Changed 9 months ago by dimpase

  • Dependencies changed from #30786 to #30786, #30446
  • Status changed from needs_review to needs_work

only the patch for expression.pyx is needed here, with the rest in #30446

comment:7 Changed 9 months ago by git

  • Commit changed from 5da750a7fca49b55bd131585f0c620490611834b to c2b45c0b3181ff8803e116e4c92ecf4c4bd748d9

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

f74f66csage --package update-latest: Accept package classes :standard:, :optional: etc., restrict to normal Python packages
182b3d2sage -package fix-checksum: Handle package classes, ignore non-normal packages
9a57cf6pynac update
563940edoctest from #30446
51def3csane version name
4834dc8remove upstreamed patch
214712dtarball update
6ee3113doctests for trac 28620, 30304, 30786
6161215fixes for trac #30446
c2b45c0doctest for trac 31137

comment:8 Changed 9 months ago by dimpase

  • Reviewers set to Dima Pasechnik
  • Status changed from needs_work to positive_review


comment:9 Changed 9 months ago by slelievre

Can we improve the ticket summary and ticket description?

comment:10 Changed 9 months ago by dimpase

it just adds an extra doctest, which used to fail before #30446

comment:11 Changed 9 months ago by gh-DaveWitteMorris

  • Description modified (diff)

comment:12 Changed 9 months ago by slelievre

  • Description modified (diff)
  • Summary changed from subs fails in obscure circumstances to Fix subs failure involving integer-valued rational exponents

Thanks for the explanation.

comment:13 Changed 8 months ago by vbraun

  • Branch changed from public/31137 to c2b45c0b3181ff8803e116e4c92ecf4c4bd748d9
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.