Opened 10 years ago

Closed 10 years ago

Last modified 10 years ago

#13314 closed defect (fixed)

segmentation fault with roots over ComplexField

Reported by: Paul Zimmermann Owned by: Burcin Erocal
Priority: major Milestone: sage-5.4
Component: calculus Keywords:
Cc: Jeroen Demeyer, John Cremona Merged in: sage-5.4.beta0
Authors: Paul Zimmermann, Jeroen Demeyer Reviewers: John Cremona
Report Upstream: Fixed upstream, but not in a stable release. Work issues:
Branch: Commit:
Dependencies: #13320, #12346 Stopgaps:

Status badges

Description (last modified by Jeroen Demeyer)

this was reported to me by Benjamin Dadoun:

sage: CC=ComplexField(128)
sage: z=polygen(CC, 'z')
sage: p=340282366920938463463374607431768211457*z^84 - 446662175392119146634315708974074167296*z^80 + 279189840125316599983095153558720348160*z^76 - 110522928999596165272398266482588385280*z^72 + 31087466142071580858231009754176552960*z^68 - 6606701181888932299640001500276588544*z^64 + 1101219301431275611150985571298443264*z^60 - 147498447216213728441335387395194880*z^56 + 16134143217737051761322600788131840*z^52 - 1456690062677843407408776850964480*z^48 + 109261914698248556477361418141696*z^44 - 6829504668478678587441535778816*z^40 + 355736441058167833075415777280*z^36 - 15393873200469425036030115840*z^32 + 549832294277433090197422080*z^28 - 16038265915694612515651584*z^24 + 375931794881970935169024*z^20 - 6911153168471496130560*z^16 + 95997158355731169280*z^12 - 947428403886210560*z^8 + 5921977682886976*z^4 - 17626570956801
sage: p.roots(ring=CC)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)

/users/caramel/bdadoun/Desktop/ceval/<ipython console> in <module>()

/usr/local/sage-5.1-linux-64bit-ubuntu_12.04_lts-x86_64-Linux/local/lib/python2\
.7/site-packages/sage/rings/polynomial/polynomial_element.so in sage.rings.poly\
nomial.polynomial_element.Polynomial.roots (sage/rings/polynomial/polynomial_el\
ement.c:31093)()

/usr/local/sage-5.1-linux-64bit-ubuntu_12.04_lts-x86_64-Linux/local/lib/python2\
.7/site-packages/sage/libs/pari/gen.so in sage.libs.pari.gen.gen.polroots (sage\
/libs/pari/gen.c:35318)()

RuntimeError: Segmentation fault

spkg: http://boxen.math.washington.edu/home/jdemeyer/spkg/pari-2.5.2.p1.spkg (diff: pari-2.5.2.p1.diff)

apply 13314_polroots_test.patch

Attachments (2)

pari-2.5.2.p1.diff (2.1 KB) - added by Jeroen Demeyer 10 years ago.
13314_polroots_test.patch (5.2 KB) - added by Jeroen Demeyer 10 years ago.

Download all attachments as: .zip

Change History (31)

comment:1 Changed 10 years ago by Leif Leonhardy

I get the same with Sage 5.2 [rc0] built from source on Ubuntu 10.04.4 LTS x86_64 (GCC 4.7.0).

comment:2 Changed 10 years ago by Leif Leonhardy

FWIW, the segfault occurs in the PARI library, and apparently depends on the precision of the complex field (works for me with up to 64 bits).

comment:3 Changed 10 years ago by Paul Zimmermann

yes the Seg. fault occurs for a precision of 65 bits or more. However I was unable to reproduce it using the GP interpreter (with version 2.5.0). It would be easier to report upstream if we can reproduce it using GP than the Pari library.

Paul

comment:4 Changed 10 years ago by Paul Zimmermann

I tried to trace the error down to the Pari library. In file rootpol.c, function roots_aux, the segmentation fault happens in the call to roots_com with l=4 and bit_accuracy(l)=128 where p is (as printed by the Pari output function):

(1.00000000000000000000000000000000000000+0.E-38*I)*z^84+(0.E-38+0.E-38*I)*z^83+(0.E-38+0.E-38*I)*z^82+(0.E-38+0.E-38*I)*z^81+(-1.31262215974857454625635000411421060562+0.E-38*I)*z^80+(0.E-38+0.E-38*I)*z^79+(0.E-38+0.E-38*I)*z^78+(0.E-38+0.E-38*I)*z^77+(0.82046519968571818282043750514276325703+0.E-38*I)*z^76+(0.E-38+0.E-38*I)*z^75+(0.E-38+0.E-38*I)*z^74+(0.E-38+0.E-38*I)*z^73+(-0.324797696688395171049634768678515683860+0.E-38*I)*z^72+(0.E-38+0.E-38*I)*z^71+(0.E-38+0.E-38*I)*z^70+(0.E-38+0.E-38*I)*z^69+(0.091357852078460689143613038254443381448+0.E-38*I)*z^68+(0.E-38+0.E-38*I)*z^67+(0.E-38+0.E-38*I)*z^66+(0.E-38+0.E-38*I)*z^65+(-0.0194153497922034252412722132863365231970+0.E-38*I)*z^64+(0.E-38+0.E-38*I)*z^63+(0.E-38+0.E-38*I)*z^62+(0.E-38+0.E-38*I)*z^61+(0.00323619266962232567325444265726730463939+0.E-38*I)*z^60+(0.E-38+0.E-38*I)*z^59+(0.E-38+0.E-38*I)*z^58+(0.E-38+0.E-38*I)*z^57+(-0.000433458978644296491914756665198016349904+0.E-38*I)*z^56+(0.E-38+0.E-38*I)*z^55+(0.E-38+0.E-38*I)*z^54+(0.E-38+0.E-38*I)*z^53+(4.74139855195190713770398956497520437380E-5+0.E-38*I)*z^52+(0.E-38+0.E-38*I)*z^51+(0.E-38+0.E-38*I)*z^50+(0.E-38+0.E-38*I)*z^49+(-4.2808273489419220927898795386466919422E-6+0.E-38*I)*z^48+(0.E-38+0.E-38*I)*z^47+(0.E-38+0.E-38*I)*z^46+(0.E-38+0.E-38*I)*z^45+(3.21091908719544601469774304001807661846E-7+0.E-38*I)*z^44+(0.E-38+0.E-38*I)*z^43+(0.E-38+0.E-38*I)*z^42+(0.E-38+0.E-38*I)*z^41+(-2.00701103917778153737692276628195892519E-8+0.E-38*I)*z^40+(0.E-38+0.E-38*I)*z^39+(0.E-38+0.E-38*I)*z^38+(0.E-38+0.E-38*I)*z^37+(1.04541544211375485185820631903782289784E-9+0.E-38*I)*z^36+(0.E-38+0.E-38*I)*z^35+(0.E-38+0.E-38*I)*z^34+(0.E-38+0.E-38*I)*z^33+(-4.52385274610660402740066744025679312462E-11+0.E-38*I)*z^32+(0.E-38+0.E-38*I)*z^31+(0.E-38+0.E-38*I)*z^30+(0.E-38+0.E-38*I)*z^29+(1.61581188955695039273058373628568909992E-12+0.E-38*I)*z^28+(0.E-38+0.E-38*I)*z^27+(0.E-38+0.E-38*I)*z^26+(0.E-38+0.E-38*I)*z^25+(-4.71322274522116476238604321139902691134E-14+0.E-38*I)*z^24+(0.E-38+0.E-38*I)*z^23+(0.E-38+0.E-38*I)*z^22+(0.E-38+0.E-38*I)*z^21+(1.10476425294559942483778527936674858155E-15+0.E-38*I)*z^20+(0.E-38+0.E-38*I)*z^19+(0.E-38+0.E-38*I)*z^18+(0.E-38+0.E-38*I)*z^17+(-2.03100537680144918840320161959898939477E-17+0.E-38*I)*z^16+(0.E-38+0.E-38*I)*z^15+(0.E-38+0.E-38*I)*z^14+(0.E-38+0.E-38*I)*z^13+(2.82110293355386359309869812795381913797E-19+0.E-38*I)*z^12+(0.E-38+0.E-38*I)*z^11+(0.E-38+0.E-38*I)*z^10+(0.E-38+0.E-38*I)*z^9+(-2.78424184144204274350181426076550186938E-21+0.E-38*I)*z^8+(0.E-38+0.E-38*I)*z^7+(0.E-38+0.E-38*I)*z^6+(0.E-38+0.E-38*I)*z^5+(1.74031282798232506192330105295904480920E-23+0.E-38*I)*z^4+(0.E-38+0.E-38*I)*z^3+(0.E-38+0.E-38*I)*z^2+(0.E-38+0.E-38*I)*z+(-5.1799836460219446701130005742251516370E-26+0.E-38*I)

Paul

comment:5 Changed 10 years ago by Paul Zimmermann

investigation is underway with the help of the Pari developers.

Paul

comment:6 Changed 10 years ago by Paul Zimmermann

for the record, here are the instructions from the Pari developers to (try to) reproduce the problem with Pari: first remove the static declaration for the function roots_com, then add gpwritebin("/tmp/gp.tmp", p); before the roots_com call in rootpol.c, then under GP:

? install("roots_com", GL)
? pol = read("/tmp/gp.tmp");
> ? roots_com(pol, 4)
  ***   bug in PARI/GP (Segmentation Fault), please report
  ***   Break loop: type 'break' to go back to GP

Paul

comment:7 Changed 10 years ago by Paul Zimmermann

we are unable to reproduce the Segmentation fault from GP:

? install("roots_com", GL)
? pol = read("/tmp/gp.tmp");
? var = roots_com(p, 128)
%2 = [0.E-38]

Paul

comment:8 Changed 10 years ago by Paul Zimmermann

Report Upstream: N/AReported upstream. Developers acknowledge bug.

Karim Belabas tells me he can now reproduce the Segmentation fault with Pari 2.5.2, but not with 2.6.*.

Paul

comment:9 Changed 10 years ago by Paul Zimmermann

Report Upstream: Reported upstream. Developers acknowledge bug.Fixed upstream, but not in a stable release.

the bug is fixed with the following Pari patch (from Karim Belabas):

    fix typo introduced in 7b985ff8a726e89191d2733b9052b2e2c228bfde
    "minor improvements root_error()"
    prod was not re-initialized to 1.

diff --git a/src/basemath/rootpol.c b/src/basemath/rootpol.c
index 3099bb2..0be400f 100644
--- a/src/basemath/rootpol.c
+++ b/src/basemath/rootpol.c
@@ -1611,8 +1611,8 @@ split_0(GEN p, long bit, GEN *F, GEN *G)
 static GEN
 root_error(long n, long k, GEN roots_pol, long pari_err, GEN shatzle)
 {
-  GEN rho, d, eps, epsbis, eps2, prod, aux, rap = NULL;
-  long i, j, m;
+  GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
+  long i, j;

   d = cgetg(n+1,t_VEC);
   for (i=1; i<=n; i++)
@@ -1627,11 +1627,11 @@ root_error(long n, long k, GEN roots_pol, long pari_err\
, GEN shatzle)
   if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
   eps = mulrr(rho, shatzle);
   aux = shiftr(powru(rho,n), pari_err);
-  prod = NULL; /* 1. */

   for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
   {
-    m = n;
+    GEN prod = NULL; /* 1. */
+    long m = n;
     epsbis = mulrr(eps, dbltor(1.25));
     for (i=1; i<=n; i++)
     {

Note this bug is fixed in the "master" branch upstream (2.6.*, no official release yet).

Paul

PS: I leave for 3 weeks of vacations, feel free to submit a proper patch in the meantime.

comment:10 Changed 10 years ago by Paul Zimmermann

since nobody did a patch, I will try. How to do a patch to a standard spkg? Should I build a new one?

Paul

comment:11 in reply to:  10 Changed 10 years ago by Burcin Erocal

Cc: Jeroen Demeyer added

Replying to zimmerma:

since nobody did a patch, I will try.

Jeroen usually merges such fixes to the pari package immediately. Perhaps he hasn't seen this yet.

How to do a patch to a standard spkg? Should I build a new one?

The usual way is to create a new package and bump the revision number, from pari-2.5.2.p0 to pari-2.5.2.p1 for example.

Thanks for following up on this.

comment:12 Changed 10 years ago by Paul Zimmermann

Authors: Paul Zimmermann
Status: newneeds_review

I've put a new spkg at http://www.loria.fr/~zimmerma/pari-2.5.1.p4.spkg. Please be cool with me since it is the first time I update a spkg.

Paul

comment:13 in reply to:  12 ; Changed 10 years ago by Jeroen Demeyer

Could you please put a link to the upstream bug tracker for futute reference?

For the spkg, this should be based on the most recent version, which is pari-2.5.2. I can take care of that.

comment:14 Changed 10 years ago by Jeroen Demeyer

Except for not committing the changes the to hg repository, it seems you made the spkg correctly.

comment:15 Changed 10 years ago by Jeroen Demeyer

Description: modified (diff)

Changed 10 years ago by Jeroen Demeyer

Attachment: pari-2.5.2.p1.diff added

comment:16 Changed 10 years ago by Jeroen Demeyer

Description: modified (diff)

comment:17 in reply to:  13 Changed 10 years ago by Paul Zimmermann

Jeroen,

Replying to jdemeyer:

Could you please put a link to the upstream bug tracker for futute reference?

For the spkg, this should be based on the most recent version, which is pari-2.5.2. I can take care of that.

yes please do!

Paul

comment:18 Changed 10 years ago by Jeroen Demeyer

Authors: Paul ZimmermannPaul Zimmermann, Jeroen Demeyer
Description: modified (diff)

Added a regression test ;-) Needs review.

Changed 10 years ago by Jeroen Demeyer

Attachment: 13314_polroots_test.patch added

comment:19 Changed 10 years ago by Jeroen Demeyer

There was one doctest failure on 32-bit systems, fixed now.

Still needs review.

comment:20 Changed 10 years ago by Paul Zimmermann

Cc: John Cremona added

John, please could you review that easy one?

Paul

comment:21 in reply to:  20 Changed 10 years ago by John Cremona

Replying to zimmerma:

John, please could you review that easy one?

Paul

OK -- should I install the new spkg using sage -f and then apply the patch and then test? I think other spkgs depend on the pari one (e.g. eclib does) so I am not sure how that will work.

John

comment:22 Changed 10 years ago by Jeroen Demeyer

First of all: reviewing a ticket is mainly about looking at the patch and checking that it makes sense.

If eclib uses PARI to do root-finding over the complex numbers, then yes, it could affect eclib. However, assuming eclib does dynamic linking, then it's not required to actually recompile eclib. The upgraded PARI would automatically be used by eclib.

So, for a ticket like this, I would suggest to run all doctests in Sage (root-finding affects many things) but I don't think a full rebuild is needed.

In any case, any failures would be caught by the Buildbot. So the worst that can happen is that you waste the time of the release manager and possibly anger him slightly if the failure is particularly obvious (e.g. if somebody edits file X and there are doctests failures in file X, then the author/reviewer did not do their due diligence).

comment:23 Changed 10 years ago by John Cremona

Thanks. I did what I said (built the new spkg with sage -f, then sage -b to rebuild which did trigger a lot of compilations, then applied the patch). The touched files test OK, and I'm now testing everything (i.e. sage -t -long).

The last thing I would want to do is anger the release manager (even if only slightly)!

The Sage patch looks fine, and also the spkg patching (assuming that Karim has fixed this properly, which seems likely).

comment:24 Changed 10 years ago by John Cremona

All good except:

sage -t -long devel/sage-main/sage/matrix/matrix_integer_dense.pyx
**********************************************************************
File "/home/jec/sage-5.3.beta1/devel/sage-main/sage/matrix/matrix_integer_dense.pyx", line 4967:
    sage: pari('mathnf(Mat([0,1]), 4)')
Expected:
    [Mat([0, 1]), [1, 0; 0, 1]]
Got:
    [Mat(1), [1, 0; 0, 1]]
**********************************************************************
File "/home/jec/sage-5.3.beta1/devel/sage-main/sage/matrix/matrix_integer_dense.pyx", line 5040:
    sage: pari('mathnf(Mat([0,1]), 4)')
Expected:
    [Mat([0, 1]), [1, 0; 0, 1]]
Got:
    [Mat(1), [1, 0; 0, 1]]
**********************************************************************

which look strange to me (and also not relevant to the changes made)!

comment:25 in reply to:  24 Changed 10 years ago by Jeroen Demeyer

Dependencies: #13320, #12346

Replying to cremona:

All good except:

sage -t -long devel/sage-main/sage/matrix/matrix_integer_dense.pyx
**********************************************************************
File "/home/jec/sage-5.3.beta1/devel/sage-main/sage/matrix/matrix_integer_dense.pyx", line 4967:
    sage: pari('mathnf(Mat([0,1]), 4)')
Expected:
    [Mat([0, 1]), [1, 0; 0, 1]]
Got:
    [Mat(1), [1, 0; 0, 1]]
**********************************************************************
File "/home/jec/sage-5.3.beta1/devel/sage-main/sage/matrix/matrix_integer_dense.pyx", line 5040:
    sage: pari('mathnf(Mat([0,1]), 4)')
Expected:
    [Mat([0, 1]), [1, 0; 0, 1]]
Got:
    [Mat(1), [1, 0; 0, 1]]
**********************************************************************

Apply the patch from #12346 (merged in sage-5.3.beta2, I assume you use an older version).

comment:26 Changed 10 years ago by John Cremona

Sorry, I was using beta1 and did not know about the dependencies. Watch this space...

comment:27 in reply to:  26 Changed 10 years ago by John Cremona

Reviewers: John Cremona
Status: needs_reviewpositive_review

Replying to cremona:

Sorry, I was using beta1 and did not know about the dependencies. Watch this space...

OK now with that patch!

comment:28 Changed 10 years ago by Jeroen Demeyer

Merged in: sage-5.4.beta0
Milestone: sage-5.3sage-5.4
Resolution: fixed
Status: positive_reviewclosed

comment:29 Changed 10 years ago by Paul Zimmermann

thank you very much John for the efficient review!

Paul

Note: See TracTickets for help on using tickets.