Infinite loop in gcd() via pynac-0.2.1

[sagetrac-bgoodri]

This bug was fixed in GiNaC 1.5.x but pynac forked off the 1.4.x branch. See

http://www.ginac.de/ginac.git?p=ginac.git;a=commit;h=edf1ae46a926d0a718063c149b78c1b9a7ec2043

Translated to sage syntax from the GiNaC 1.5.x unit tests:

u = var('u')
v = var('v')
w = var('w')
x = var('x')
y = var('y')
z = var('z')

e = 792*z^8*w^4*x^3*y^4*u^7 + 24*z^4*w^4*x^2*y^3*u^4 + \
264*z^8*w^3*x^2*y^7*u^5 + 198*z^4*w^5*x^5*y*u^6  + 110*z^2*w^3*x^5*y^4*u^6 \
- 120*z^8*w*x^4*u^6 - 480*z^5*w*x^4*y^6*u^8 - 720*z^7*x^3*y^3*u^7 + \
165*z^4*w^2*x^4*y*u^5 + 450*z^8*w^6*x^2*y*u^8 + 40*z^2*w^3*x^3*y^3*u^6 - \
288*z^7*w^2*x^3*y^6*u^6  + 250*z^6*w^4*x^2*y^4*u^8 + \
576*z^7*w^7*x^2*y^4*u^8  - 80*z^6*w^2*x^5*y^3*u^7 - 144*z^8*w^4*x^5*u^7 + \
120*z^4*w*x^2*y^6*u^6 + 320*z^5*w^5*x^2*y^7*u^8 + 192*z^7*w^6*x*y^7*u^6 - \
12*z^4*w^3*x^3*y^5*u^6  - 36*z^4*w^4*x^4*y^2*u^8 + 72*z^4*w^5*x^3*u^6  - \
20*z^2*w^2*x^4*y^5*u^8 + 660*z^8*w*x^2*y^4*u^6 + 66*z^4*w^4*x^4*y^4*u^4 + \
440*z^6*w^2*x^3*y^7*u^7  - 30*z^4*w*x^3*y^2*u^7 - 48*z^8*w^3*x^4*y^3*u^5 + \
72*z^6*w^2*x*y^6*u^4 - 864*z^7*w^3*x^4*y^3*u^8 + 480*z^7*w^4*x*y^4*u^7 + \
60*z^4*w^2*x^2*u^5 + 375*z^8*w^3*x*y*u^7 + 150*z^8*w^5*x*y^4*u^6 + \
180*z^6*x*y^3*u^5 + 216*z^6*w^3*x^2*y^3*u^6;

E = e.polynomial(QQ)
G = gcd(E, e.diff(x).polynomial(QQ)) 
G # Singular gets it correct: u^4*z^2

g = gcd(e, e.diff(x)) # pynac falls into an infinite loop 
# also there is potential for a wrong answer

So, it appears as if we need to use Singular or sync to at least GiNaC 1.5.8 .

Depends on #21963

Upstream: Fixed upstream, in a later stable release.

Component: symbolics

Keywords: gcd

Reviewer: Frédéric Chapoton

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

[burcin]

Replying to @sagetrac-bgoodri:

This bug was fixed in GiNaC 1.5.x but pynac forked off the 1.4.x branch. See

http://www.ginac.de/ginac.git?p=ginac.git;a=commit;h=edf1ae46a926d0a718063c149b78c1b9a7ec2043

> So, it appears as if we need to use Singular or sync to at least GiNaC 1.5.8 .

If there isn't a straightforward way to adopt the patch that went into GiNaC 1.5.x, I don't think it's worth the time to sync the gcd code in pynac with the latest version in GiNaC.

Multivariate gcd is a hard problem. IMHO, we should have a single, well-tested, fast implementation and use that everywhere in Sage, rather than rely on the efforts of each upstream developer team to implement their own with a unique set of bugs.

Using Factory as a library doesn't look too hard:

http://www.singular.uni-kl.de/svn/trunk/factory/examples/gcd.cc

The README file is helpful, although I'm not sure if there is any other documentation available:

http://www.singular.uni-kl.de/svn/trunk/factory/README

I'm quite sure we would get good feedback from the current developers if the issues are discussed at the libsingular-devel list:

http://groups.google.com/group/libsingular-devel

It seems that once the necessary changes are made in pynac to use Factory for gcd's, it would be straightforward to plug in another library, perhaps giac, later.

[sagetrac-bgoodri]

comment:2

Replying to @burcin:

All valid points. Maybe we should move the discussion to sage-devel to get more input?

If you want a quick hack fix for this particular bug, I think it would suffice to add some logic to py_gcd to coerce rational expressions with a .polynomial(QQ) and call gcd() from sage/rings/polynomial/multi_polynomial_libsingular instead of gcd() from sage/rings/arith . I tried that a couple of days ago and it was insufficient to fix #10268, but it would presumably fix this bug. Not that #10268 is a killer feature but GiNaC does seem to have a significant speed advantage over Maxima at simplifying rationals.

As for the more general GiNaC / pynac vs. Singular / Factory comparison (not to mention giac), I haven't investigated enough to have a super strong opinion about it. I agree it would be easier to maintain if there were only one interface. On the other hand, it is good to make it possible to use whatever might be fastest and have the most features for a particular problem. Also, in both cases, we are lagging behind upstream, but they have a lot of the same algorithms now. Singular 3.1-2 fixed / improved a lot of multivariate polynomial stuff but sage currently ships with 3.0-something. Similarly, GiNaC 1.5.8 fixed this bug and 1.5.0 added a lot of functionality in the polynomials/ directory, but pynac is based on 1.4-something.

What I like about (the latest) GiNaC is that it is closer to the sage SR concept, whereas Singular is more literal about its rings. For example, you can do this

goodrich@Y560:/tmp$ ginsh
ginsh - GiNaC Interactive Shell (ginac V1.5.8)
  __,  _______  Copyright (C) 1999-2010 Johannes Gutenberg University Mainz,
 (__) *       | Germany.  This is free software with ABSOLUTELY NO WARRANTY.
  ._) i N a C | You are welcome to redistribute it under certain conditions.
<-------------' For details type `warranty;'.

Type ?? for a list of help topics.
> f = Pi^2 - 1;
-1+Pi^2
> g = Pi + 1;
1+Pi
> normal(f/g);
-1+Pi
> quit;

and similarly with other transcendental numbers and trig and so forth. AFAIK, Singular would throw an error about QQness instead of determining that Pi - 1 was the greatest common factor.

There is a lot to like about Singular too, which is probably why William wrote those things in TODO.txt and the docstrings. Some of that seems to have changed in the last couple of years though.