Crash when expanding square root of -x

[EmmanuelCharpentier]

Before this ticket, this crashes Sage:

sage: assume(x < 0)
sage: sqrt(-x).expand()

See comment:12 for an explanation of the cause.

---

Original title and description:

Some algebraic substitutions may crash Sage

History : see this sage-release post and the following.

TL;DR: Sage crashes reproducibly on some substitutions:

f = function("f")
var("k")
Eq = f(x).diff(x, 2)/f(x) == k
with assuming(k > 0):
    Sp = desolve(Eq, f(x), ivar=x)
with assuming(k == 0):
    Sz = desolve(Eq, f(x), ivar=x)
with assuming(k < 0):
    Sn = desolve(Eq, f(x), ivar=x)
var("_K1, _K2")
with assuming(k < 0):
    Sn1 = Sn.subs(sqrt(k) == I*sqrt(abs(k)))
# Crashes Sage

This may result from an incorrect backtranslation from Maxima:

$ sage
... SageMath version 9.2.beta14, Release Date: 2020-09-30
... Using Python 3.8.6. Type "help()" for help.
...
sage: var("k")
k
sage: assume(k < 0)
sage: f = function("f")
sage: Eq = f(x).diff(x, 2)/f(x) == k
sage: S = desolve(Eq, f(x), ivar=x); S
_K2*cosh(sqrt(k)*x) + I*_K1*sinh(sqrt(k)*x)
sage: SS = maxima.subst([sqrt(k) == I*sqrt(-k)], S); SS
_SAGE_VAR__K2*cos(sqrt(-_SAGE_VAR_k)*_SAGE_VAR_x)-_SAGE_VAR__K1*sin(sqrt(-_SAGE_VAR_k)*_SAGE_VAR_x)
sage: SS.sage()
# CRASH:
------------------------------------------------------------------------
/usr/local/sage-9/local/lib/python3.8/site-packages/cysignals/signals.cpython-38-x86_64-linux-gnu.so(+0x842b)[0x7fac8352242b]
# Etc. ad nauseam...

I mark this as blocker since it crashes Sage on (relatively) trivial algebraic operations.

Interestingly, I can't reproduce the problem in expression coming from another source:

$ sage
... SageMath version 9.2.beta14, Release Date: 2020-09-30
... Using Python 3.8.6. Type "help()" for help.
...
sage: k = var("k")
sage: assume(k < 0)
sage: f = function("f")
sage: Eq = f(x).diff(x, 2)/f(x) == k
sage: from sympy import sympify, dsolve
sage: S = dsolve(*map(sympify, (Eq, f(x)))); S
Eq(f(x), C1*exp(-sqrt(k)*x) + C2*exp(sqrt(k)*x))
sage: S._sage_()
f(x) == C2*e^(sqrt(k)*x) + C1*e^(-sqrt(k)*x)
sage: S.subs({sympify(sqrt(k)): sympify(I*sqrt(-k))})
Eq(f(x), C1*exp(-I*x*sqrt(-k)) + C2*exp(I*x*sqrt(-k)))
sage: S.subs({sympify(sqrt(k)): sympify(I*sqrt(-k))})._sage_()
f(x) == C2*e^(I*sqrt(-k)*x) + C1*e^(-I*sqrt(-k)*x)
sage: S._sage_().subs(sqrt(k) == I*sqrt(-k))
f(x) == C2*e^(I*sqrt(-k)*x) + C1*e^(-I*sqrt(-k)*x)

Since this problem appears in compiled code, I don't (yet) know how to tackle it. Advice welcome...

Depends on #30446

CC: @kcrisman @slel

Component: symbolics

Author: Dave Morris

Branch: 20067d3

Reviewer: Dima Pasechnik

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

[EmmanuelCharpentier]

comment:1

Not the backtranslation from Maxima (or not only) ; I have the same problem backtranslating from Sympy :

charpent@zen-book-flip:~$ sage
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 9.2.beta14, Release Date: 2020-09-30              │
│ Using Python 3.8.6. Type "help()" for help.                        │
└────────────────────────────────────────────────────────────────────┘
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Warning: this is a prerelease version, and it may be unstable.     ┃
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
sage: k=var("k")                                                                
sage: assume(k<0)                                                               
sage: f=function("f")                                                           
sage: Eq=f(x).diff(x,2)/f(x)==k                                                 
sage: SM=desolve(Eq, f(x), ivar=x)                                              
sage: SM                                                                        
_K2*cosh(sqrt(k)*x) + I*_K1*sinh(sqrt(k)*x)
sage: from sympy import sympify                                                 
sage: sympify(SM).subs({sympify(sqrt(k)):sympify(I*sqrt(-k))})                  
-_K1*sin(x*sqrt(-k)) + _K2*cos(x*sqrt(-k))
sage: sympify(SM).subs({sympify(sqrt(k)):sympify(I*sqrt(-k))})._sage_()         
# CRASH : 
------------------------------------------------------------------------
/usr/local/sage-9/local/lib/python3.8/site-packages/cysignals/signals.cpython-38-x86_64-linux-gnu.so(+0x842b)[0x7f5dcdf4e42b]

This seems to point to pynac...

[EmmanuelCharpentier]

comment:2

#30378 seems relevant.

[EmmanuelCharpentier]

comment:3

The error is trigerred by the existence of the assumption about k, which might trigger further transformations (which turn out to be toxic in this case).

This works:

f = function("f")
var("k")
Eq = f(x).diff(x, 2)/f(x) == k
with assuming(k > 0):
    Sp = desolve(Eq, f(x), ivar=x)
with assuming(k == 0):
    Sz = desolve(Eq, f(x), ivar=x)
with assuming(k < 0):
    Sn = desolve(Eq, f(x), ivar=x)
Sn1 = Sn.subs(sqrt(k) == I*sqrt(abs(k)))

This crashes Sage:

assume(k < 0)
Sn2 = Sn.subs(sqrt(k) == I*sqrt(abs(k)))

Ditto when using with assuming....

[fchapoton]

comment:4

Given that our release is long overdue, and nobody works on the present ticket, I turn it to critical. ANybody proposing a branch in the next few days can of course reconsider this.

[EmmanuelCharpentier]

comment:5

Replying to @fchapoton:

Given that our release is long overdue, and nobody works on the present ticket, I turn it to critical. ANybody proposing a branch in the next few days can of course reconsider this.

I'd be delighted to work on this, but the problem resides in Cythonized library code, the care and feeding of which I haven't the slightest idea...

Do you have pointers to an introduction to working on Sage library code ?

[EmmanuelCharpentier]

comment:6

A very quick check with Cocalc's available kernels (back to 8.2) shows that the following:

from __future__ import print_function
print(sage.version.version)
f = function("f")
var("k")
Eq = f(x).diff(x, 2)/f(x) == k
with assuming(k < 0):
    Sn = desolve(Eq, f(x), ivar=x)
var("K1,_K2")
with assuming (k < 0):
    print(Sn.subs(sqrt(k) == I*sqrt(abs(k))))

crashed even in 8.2 days...
removing the with assuming (x < 0): allows it to run.

[EmmanuelCharpentier]

comment:7

Mark as blocker, but for 9.3...

[mjungmath]

comment:8

Might be related to #30446.

[dimpase]

comment:9

Replying to @mjungmath:

Might be related to #30446.

still the example from comment:6 with #30446 installed leads to a crash.

It certainly should not be a blocker, though. There are many other ways to crash Sage, after all.

[DaveWitteMorris]

comment:10

Here is a shorter session that demonstrates the crash:

sage: assume(x < 0)
sage: sin(sqrt(-x))
------------------------------------------------------------------------
(no backtrace available)
------------------------------------------------------------------------
Unhandled SIGSEGV: A segmentation fault occurred.
    ...

On CoCalc (v9.2), this hangs, instead of crashing, and the memory meter quickly gets over 1GB. Pynac does seem to be the problem (as suggested in comment:1), because I got this trace when I interrupted the process:

Kernel last terminated by signal SIGINT.
...
093e2c1f06]
/ext/sage/sage-9.2/local/lib/libpynac.so.18(_ZNK5GiNaC5power6expandEj+0x3d8)[0x7f093e38fcb8]
/ext/sage/sage-9.2/local/lib/libpynac.so.18(_ZNK5GiNaC2ex6expandEj+0x36)[0x7f093e2c1f06]
/ext/sage/sage-9.2/local/lib/libpynac.so.18(_ZNK5GiNaC5power6expandEj+0x3d8)[0x7f093e38fcb8]
/ext/sage/sage-9.2/local/lib/libpynac.so.18(_ZNK5GiNaC2ex6expandEj+0x36)[0x7f093e2c1f06]
    <these two lines repeated many more times, then>
------------------------------------------------------------------------
Attaching gdb to process id 7387.
Traceback (most recent call last):
  File "/ext/sage/sage-9.2/local/bin/cysignals-CSI", line 237, in <module>
    main(args)
  File "/ext/sage/sage-9.2/local/bin/cysignals-CSI", line 186, in main
    trace = run_gdb(args.pid, not args.nocolor)
  File "/ext/sage/sage-9.2/local/bin/cysignals-CSI", line 110, in run_gdb
    stdout, stderr = cmd.communicate(gdb_commands(pid, color))
  File "/usr/lib/python3.8/subprocess.py", line 1024, in communicate
    stdout, stderr = self._communicate(input, endtime, timeout)
  File "/usr/lib/python3.8/subprocess.py", line 1866, in _communicate
    ready = selector.select(timeout)
  File "/usr/lib/python3.8/selectors.py", line 415, in select
    fd_event_list = self._selector.poll(timeout)
KeyboardInterrupt

[DaveWitteMorris]

comment:11

I think the problem is even worse than it appeared.

sage: assume(x < 0)                                                                    
sage: ((-x)^(1/2)).expand()                                                            
---------------------------------------------------------------------------
SignalError                               Traceback (most recent call last)
<ipython-input-2-1d555cbb3a07> in <module>
----> 1 ((-x)**(Integer(2)/Integer(3))).expand()

~/development/sage-9.2/local/lib/python3.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.expand (build/cythonized/sage/symbolic/expression.cpp:30019)()
   4820 
   4821         cdef GEx x
-> 4822         sig_on()
   4823         try:
   4824             x = self._gobj.expand(0)

SignalError: Segmentation fault

I think you can replace 1/2 with 2/3 or 3/5 or pretty much any other non-integer fraction, and still get a segfault.

[DaveWitteMorris]

comment:12

I think I found the bug. The definition of the power::expand method in ginac's power.cpp file includes the following snippet:

// search for positive/negative factors
                for (const auto & elem : m.seq) {
	const ex& e = m.recombine_pair_to_ex(elem);
	if (e.is_positive())
		prodseq.push_back(pow(e, exponent).expand(options));
	else if (e.info(info_flags::negative)) {
		prodseq.push_back(pow(-e, exponent).expand(options));
		possign = !possign;
	} else
		powseq.push_back(elem);
}

This is part of the code that calculates the expansion of a power of a product: (a1 * a2 * ... * an)^c, when c is not an integer. The else if clause says that if some factor ai is negative, then we first need to calculate the expansion of (-ai)^c. However, suppose n = 2 and a1 = -1. Then this says that before we calculate (-a2)^c, we need to calculate (-a2)^c. This is an infinite loop.

I think it will be easy to fix. I can't do it right now, but I should have time to write a patch in the next couple of days.

[dimpase]

comment:13

Thanks for working on this. As you know, pynac is based on an old fork of ginac - so the question is also whether this bug was fixed upstream.

[dimpase]

comment:14

This is what one has in Ginac upstream: https://ginac.de/reference/power_8cpp_source.html

        // search for positive/negative factors
         for (auto & cit : m.seq) {
             ex e=m.recombine_pair_to_ex(cit);
             if (e.info(info_flags::positive))
                 prodseq.push_back(pow(e, exponent).expand(options));
             else if (e.info(info_flags::negative)) {
                 prodseq.push_back(pow(-e, exponent).expand(options));
                 possign = !possign;
             } else
                 powseq.push_back(cit);
         }

Looks like it's still basically the same, but I'm not 100% sure.

[DaveWitteMorris]

comment:15

I don't know c++ or ginac, but I don't see any meaningful difference, so I don't think the ginac upstream fixed this.

Patching pynac to comment out the else if clause does eliminate the errors, so I am confident that this is indeed the bug we were looking for. However, it might be better to fix this else if clause, instead of just removing it. For example, with the clause removed, we have

sage: var("a b");
sage: assume(a < 0)
sage: assume(b < 0)
sage: sqrt(a*b).expand()
sqrt(a*b)

If we fix the clause, then I think the output would instead be sqrt(-a) * sqrt(-b). I don't know whether that's really an improvement, but I think it's what the developers of ginac and pynac intended.

I will think about how to do that, but it's great if someone who knows something about ginac (and/or c++) wants to chime in.

[DaveWitteMorris]

Branch: public/30688

[DaveWitteMorris]

comment:17

Setting to "needs review" to ask the patchbots whether they think that deleting the else if clause is a viable option (and to encourage comments from human reviewers).

dimpase already knows this, but for others who are following the thread: As pointed out by mkoeppe in #30446 comment:17, we will actually need a new pynac tarball, not just a patch, when we have decided how to fix this. But I think a patch is good enough for now.

---

New commits:

050b28e

trac 30688 patch infinite loop in pynac

[DaveWitteMorris]

Commit: 050b28e

[dimpase]

comment:20

Opened pynac/pynac#368 - hopefully people knowing GiNaC code can comment.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

6ea4bc6

doctests and minor edits

[sagetrac-git]

Changed commit from 050b28e to 6ea4bc6

[DaveWitteMorris]

comment:22

I suggest we close this critical ticket by applying my patch that comments out the else if clause. A follow up ticket can discuss a better way to handle factors that are negative, but that ticket would have minor priority.

I don't know anything about packaging and making tarballs. @dimpase, would that be easy for you to do? (Or perhaps add it to the tarball in #30446?) I could probably figure out how to post a PR to https://github.com/pynac/pynac, if that would be helpful.

I downloaded vanilla ginac and did some testing in C++. I was unable to get an infinite loop. It may be there, but it is definitely harder to access there than in sage. For example, vanilla ginac does not provide a straightforward way to declare a variable to be negative. (Of course, my lack of expertise is also an issue.)

What I don't understand is why only negative terms are a problem.

The expression e (or what I called ai) is only one term, so it is definitely shorter than the entire product, so the recursive call pow(e, exponent).expand(options) is safe. That is why positive terms do not give an infinite loop. The problem with negative terms is that they call pow(-e, exponent).expand(options). But -e (in other words, -1 * e) might be the entire product, in which case this is an infinite loop.

[dimpase]

comment:23

Replying to @DaveWitteMorris:

I suggest we close this critical ticket by applying my patch that comments out the else if clause. A follow up ticket can discuss a better way to handle factors that are negative, but that ticket would have minor priority.

Commenting out the whole else if codebranch would result in some expressions being left
unexpanded, no?

I don't know anything about packaging and making tarballs. @dimpase, would that be easy for you to do? (Or perhaps add it to the tarball in #30446?) I could probably figure out how to post a PR to https://github.com/pynac/pynac, if that would be helpful.

we can just do a package patch, no need to make a new tarball.
And this patch is just the output of git diff in the right repo,
stored as build/pkgs/pynac/patches/<name>.patch

I downloaded vanilla ginac and did some testing in C++. I was unable to get an infinite loop. It may be there, but it is definitely harder to access there than in sage. For example, vanilla ginac does not provide a straightforward way to declare a variable to be negative. (Of course, my lack of expertise is also an issue.)

Could you put your try somewhere on the github?
It'd save me time to try to get a better test case.

What I don't understand is why only negative terms are a problem.

The expression e (or what I called ai) is only one term, so it is definitely shorter than the entire product, so the recursive call pow(e, exponent).expand(options) is safe. That is why positive terms do not give an infinite loop. The problem with negative terms is that they call pow(-e, exponent).expand(options). But -e (in other words, -1 * e) might be the entire product, in which case this is an infinite loop.

yes, this looks to be the case. A proper way to fix this might be to get the
expression "length" in this case and if it is 2, tag it as "expanded".
(Doing this in C++ is a bit of work - mostly to understand the API)

[DaveWitteMorris]

comment:24

I think that almost every calculation that uses the else if clause gives a crash now. Since there haven't been any complaints until recently, it must be very uncommon. (Further evidence is that the patch doesn't give any doctest failures on my computer with sage -t -l -p 2 --optional=sage --all.) This is really a corner case that has little or no impact, so I think that deleting the clause is the right thing to do for this ticket. It actually would not mean that the expressions would not be expanded, just that they would not be expanded in this particular way: I think sqrt(A*B) would sometimes be expanded as sqrt((A*B).expand()), instead of sqrt((-A).expand()) * sqrt((-B).expand()). However, as I said, I think this "improved" expansion has rarely been applied to a negative factor.

Do I need to do something to get the patchbots to test this PR? I edited build/pkgs/pynac/package-version.txt in this latest PR, hoping that it would help.

[DaveWitteMorris]

comment:25

I tried lots of things, but here is the best one that I could come up with by using my almost nonexistent understanding of ginac and c++:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

class negsymbol : public symbol { 
    public:
    bool info(unsigned inf) const
    {
        if (inf == info_flags::negative) {
            return 1;}
        return inherited::info(inf);
    }
};

int main()
{
    realsymbol x("x");
    negsymbol y;
    numeric r(3,5);    
    cout << pow(x * y, r).expand() << endl;
    cout << y.info(info_flags::negative) << endl;
    cout << tree << x * y ;
}

After a warning message from -Wundefined-var-template, the output is

(x*symbol3)^(3/5)
1
mul @0x7fa050c091e0, hash=0x7ece0009, flags=0x3, nops=2
    x (symbol) @0x7fa050c09190, serial=2, hash=0x4d1fe6f1, flags=0xf, domain=1
    1 (numeric) @0x7fa050c06410, hash=0x160af41d, flags=0xf
    -----
     (symbol) @0x7fa050c090b0, serial=3, hash=0x7ecebd84, flags=0xf, domain=0
    1 (numeric) @0x7fa050c06410, hash=0x160af41d, flags=0xf
    =====

I think the else if clause is not triggered, because y is stored as y^1 in the product, and ginac does not know that y^1 is negative (even though y is negative). I do not know how to get ginac to realize that y^1 is negative.

If I change negsymbol y to possymbol y, then the output is

symbol3^(3/5)*x^(3/5)
0
mul @0x7fd3add07800, hash=0x693c001c, flags=0x3, nops=2
     (symbol) @0x7fd3add076d0, serial=3, hash=0x693c85bf, flags=0xf, domain=2
    1 (numeric) @0x7fd3add04a30, hash=0x160af41d, flags=0xf
    -----
    x (symbol) @0x7fd3add077b0, serial=2, hash=0x76bbb6f1, flags=0xf, domain=1
    1 (numeric) @0x7fd3add04a30, hash=0x160af41d, flags=0xf
    =====

This means that the if clause was triggered, as it should be.

[dimpase]

Dependencies: #30446

[dimpase]

Reviewer: Dima Pasechnik

[dimpase]

comment:26

OK, I'm at loss here.

Anyway,let's get your patch in, but please rebase it over
the branch of #30446
(i.e. it should be for pynac 0.7.27, not the current one)

[DaveWitteMorris]

Changed branch from public/30688 to public/30688r1

[DaveWitteMorris]

comment:28

OK, I always having trouble doing a rebase, but I think I eventually got it correct.

I changed the package version to 0.7.27.p1 -- is that correct?

---

Last 10 new commits:

182b3d2	sage -package fix-checksum: Handle package classes, ignore non-normal packages
9a57cf6	pynac update
563940e	doctest from #30446
51def3c	sane version name
4834dc8	remove upstreamed patch
214712d	tarball update
ccdf77c	update SPKG.rst
fab5559	trac 30688 patch infinite loop in pynac
23875bd	doctests and minor edits
20067d3	update package version

[DaveWitteMorris]

Author: Dave Morris

[DaveWitteMorris]

Changed commit from 6ea4bc6 to 20067d3

[DaveWitteMorris]

comment:29

Followup ticket: I opened #31337 for reimplementing the else if clause.

[dimpase]

comment:31

OK,good.

[vbraun]

Changed branch from public/30688r1 to 20067d3

[slel]

Changed commit from 20067d3 to none