Opened 7 years ago

Closed 7 years ago

Last modified 7 years ago

#13314 closed defect (fixed)

segmentation fault with roots over ComplexField

Reported by: zimmerma Owned by: burcin
Priority: major Milestone: sage-5.4
Component: calculus Keywords:
Cc: jdemeyer, 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:

Description (last modified by jdemeyer)

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 jdemeyer 7 years ago.
13314_polroots_test.patch (5.2 KB) - added by jdemeyer 7 years ago.

Download all attachments as: .zip

Change History (31)

comment:1 Changed 7 years ago by leif

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 7 years ago by leif

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 7 years ago by zimmerma

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 7 years ago by zimmerma

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 7 years ago by zimmerma

investigation is underway with the help of the Pari developers.

Paul

comment:6 Changed 7 years ago by zimmerma

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 7 years ago by zimmerma

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 7 years ago by zimmerma

  • Report Upstream changed from N/A to Reported 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 7 years ago by zimmerma

  • Report Upstream changed from Reported upstream. Developers acknowledge bug. to 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 follow-up: Changed 7 years ago by zimmerma

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 7 years ago by burcin

  • Cc jdemeyer 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 follow-up: Changed 7 years ago by zimmerma

  • Authors set to Paul Zimmermann
  • Status changed from new to needs_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 ; follow-up: Changed 7 years ago by 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.

comment:14 Changed 7 years ago by jdemeyer

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

comment:15 Changed 7 years ago by jdemeyer

  • Description modified (diff)

Changed 7 years ago by jdemeyer

comment:16 Changed 7 years ago by jdemeyer

  • Description modified (diff)

comment:17 in reply to: ↑ 13 Changed 7 years ago by zimmerma

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 7 years ago by jdemeyer

  • Authors changed from Paul Zimmermann to Paul Zimmermann, Jeroen Demeyer
  • Description modified (diff)

Added a regression test ;-) Needs review.

Changed 7 years ago by jdemeyer

comment:19 Changed 7 years ago by jdemeyer

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

Still needs review.

comment:20 follow-up: Changed 7 years ago by zimmerma

  • Cc cremona added

John, please could you review that easy one?

Paul

comment:21 in reply to: ↑ 20 Changed 7 years ago by 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 7 years ago by jdemeyer

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 7 years ago by 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 follow-up: Changed 7 years ago by 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 7 years ago by jdemeyer

  • Dependencies set to #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 follow-up: Changed 7 years ago by cremona

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

comment:27 in reply to: ↑ 26 Changed 7 years ago by cremona

  • Reviewers set to John Cremona
  • Status changed from needs_review to positive_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 7 years ago by jdemeyer

  • Merged in set to sage-5.4.beta0
  • Milestone changed from sage-5.3 to sage-5.4
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:29 Changed 7 years ago by zimmerma

thank you very much John for the efficient review!

Paul

Note: See TracTickets for help on using tickets.