jdemeyer
commented
6 years ago comment:1
Open mezzarobba opened 8 years ago
jdemeyer
commented
6 years ago 
mezzarobba
commented
6 years ago /home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.real_roots (build/cythonized/sage/rings/polynomial/real_roots.c:47075)()
4053
4054 oc = ocean(ctx, bernstein_polynomial_factory_ratlist(b), linear_map(left, right))
-> 4055 oc.find_roots()
4056 oceans.append((oc, factor, exp))
4057
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.ocean.find_roots (build/cythonized/sage/rings/polynomial/real_roots.c:35103)()
3069 """
3070 while not self.all_done():
-> 3071 self.refine_all()
3072 self.increase_precision()
3073
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.ocean.refine_all (build/cythonized/sage/rings/polynomial/real_roots.c:35430)()
3114 while isle is not self.endpoint:
3115 if not isle.done(self.ctx):
-> 3116 isle.refine(self.ctx)
3117 isle = isle.rgap.risland
3118
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine (build/cythonized/sage/rings/polynomial/real_roots.c:37966)()
3367 self.shrink_bp(ctx)
3368 try:
-> 3369 self.refine_recurse(ctx, self.bp, self.ancestors, [], True)
3370 except PrecisionError:
3371 pass
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine_recurse (build/cythonized/sage/rings/polynomial/real_roots.c:39863)()
3515 return
3516 else:
-> 3517 self.refine_recurse(ctx, p1, ancestors, history, False)
3518 assert(self.lgap.upper == p2.lower)
3519 bp = p2
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine_recurse (build/cythonized/sage/rings/polynomial/real_roots.c:39705)()
3510 if not lisland.done(ctx):
3511 try:
-> 3512 lisland.refine_recurse(ctx, p1, ancestors, history, True)
3513 except PrecisionError:
3514 pass
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine_recurse (build/cythonized/sage/rings/polynomial/real_roots.c:39863)()
3515 return
3516 else:
-> 3517 self.refine_recurse(ctx, p1, ancestors, history, False)
3518 assert(self.lgap.upper == p2.lower)
3519 bp = p2
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine_recurse (build/cythonized/sage/rings/polynomial/real_roots.c:39705)()
3510 if not lisland.done(ctx):
3511 try:
-> 3512 lisland.refine_recurse(ctx, p1, ancestors, history, True)
3513 except PrecisionError:
3514 pass
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.refine_recurse (build/cythonized/sage/rings/polynomial/real_roots.c:38939)()
3471 rv = bp.try_rand_split(ctx, [])
3472 if rv is None:
-> 3473 (ancestors, bp) = self.more_bits(ctx, ancestors, bp, rightmost)
3474 if rv is None:
3475 rv = bp.try_split(ctx, [])
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.island.more_bits (build/cythonized/sage/rings/polynomial/real_roots.c:41622)()
3617 assert(rel_bounds[1] == 1)
3618
-> 3619 ancestor_val = split_for_targets(ctx, ancestor_val, [(self.lgap, maybe_rgap, target_lsb_h)])[0]
3620 # if rel_lbounds[1] > 0:
3621 # left_split = -exact_rational(simple_wordsize_float(-rel_lbounds[1], -rel_lbounds[0]))
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.split_for_targets (build/cythonized/sage/rings/polynomial/real_roots.c:32908)()
2880 split = wordsize_rational(split_targets[best_index][0], split_targets[best_index][1], ctx.wordsize)
2881
-> 2882 (p1_, p2_, ok) = bp.de_casteljau(ctx, split, msign=split_targets[best_index][2])
2883 assert(ok)
2884
/home/marc/co/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx in sage.rings.polynomial.real_roots.interval_bernstein_polynomial_integer.de_casteljau (build/cythonized/sage/rings/polynomial/real_roots.c:9092)()
728 msign = sign
729 elif sign != 0:
--> 730 assert(msign == sign)
731
732 cdef Rational absolute_mid = self.lower + mid * (self.upper - self.lower
mkoeppe
commented
4 years ago bug is still present in 9.2.beta2
orlitzky
commented
4 years ago FWIW you can change this to a recursion error by increasing the precision:
diff --git a/src/sage/rings/polynomial/real_roots.pyx b/src/sage/rings/polynomi$
index ba57c1ba4f..e0226fc53a 100644
--- a/src/sage/rings/polynomial/real_roots.pyx
+++ b/src/sage/rings/polynomial/real_roots.pyx
@@ -2864,6 +2864,7 @@ def split_for_targets(context ctx, interval_bernstein_pol$
goal = Integer(65535)/65536
else:
goal = Integer(255)/256
+ goal = Integer(4294967296)/Integer(4294967295)
if ntargs == 1 and not(out_of_bounds) and split_targets[1][0] - split_targ$
return [bp]
Maybe that helps debug it. Or maybe this polynomial is just too crazy.
mezzarobba
commented
4 years ago Description changed:
---
+++
@@ -8,3 +8,5 @@
AssertionError:+ +See also #26955.
Moving to 9.4, as 9.3 has been released.
mezzarobba
commented
1 year ago Here are the corresponding logs (extracted by patching the code so that it returns the context object on error):
sage: ctx.get_be_log()
[(2.0742416381835938e-05, 594, 22, True, 5, -1952891, 1, 41, 30, (3, 39))]
sage: ctx.get_dc_log()
[('halff', 736, 53, True, []),
('halff', 736, 10, False, []),
('pulling',
736,
0,
615,
615,
0.500000000000000,
682,
746,
736,
698,
682,
650),
('split_for_targets', 1, 1/2, 615, 121, 78, 121),
('down_degree', False, ('>', 691)),
('halff', 693, 53, False, []),
('divf', 693, 53, 393/1024, False, []),
('halfi', 654, 39, False, []),
('divi', 654, 39, 5/16, False, []),
('splitting', (0, 0), (1/2, 1), 615),
('split_for_targets', 1, 1/2, 495, 241, 198, 241),
('down_degree', False, ('>', 691)),
('halfi', 594, 99, True, []),
('halff', 652, 53, True, []),
('halff', 652, 16, False, []),
('pulling',
652,
0,
594,
594,
0.500000000000000,
594,
668,
652,
620,
604,
572),
('split_for_targets', 1, 1/2, 594, 58, 21, 58),
('down_degree', True, ('<', 615)),
('down_degree', False, ('$', 632)),
('halff', 616, 53, True, []),
('halff', 616, 52, True, []),
('halff', 616, 52, False, []),
('divf', 616, 52, 517/1024, False, []),
('pulling',
693,
0,
495,
495,
0.0625000000000000,
495,
668,
616,
578,
548,
488)]Does anyone understand the code, or at least the algorithm, well enough to stand a chance of understanding the issue? (@zimmermann6 maybe??) The code is quite well documented, but it is a bit much for me to digest.
zimmermann6
commented
1 year ago who is the author of that code ? Where is it located ? Is there a simpler example (I cannot copy&paste the long input) ?
mezzarobba
commented
1 year ago It is located in src/sage/rings/polynomial/real_roots.pyx and was written by Carl Witty in 2007. I'm asking because @videlec suggested you already had a look at it in the past.
mezzarobba
commented
1 year ago No simpler example, unfortunately, but here is this one as a .sage file
mezzarobba
commented
1 year ago Actually some other examples are listed at #26957, but none of them fails on my system.
However, @fchapoton mentioned a few years ago that they did fail on his. This could point at an issue with the use of hardware floating-point numbers...
zimmermann6
commented
1 year ago I cloned Sage from github (revision 84eebf0). I see two assert(msign == sign) lines in real_roots.pyx: one at line 736 and one at line 1606. Is that the one at line 736 that fails? (I cannot build that version from source.)
zimmermann6
commented
1 year ago @mezzarobba please can you create a Sage file with the inputs to the de_casteljau routine that exhibit the issue?
mezzarobba
commented
1 year ago Yes, the failing assertion is the one on line 736. Here is a script that reproduces just the failing call to de_casteljau.
de_casteljau.sage.txt
zimmermann6
commented
1 year ago thanks. Assuming the polynomial corresponds to the list [-1, -1, -1, ...] with low-order coefficients first, the given value of msign (-1) which should be its sign at mid (1/16) seems to be wrong:
sage: p
214703574267988798662235357982294318934218074726925561209802*x^41 - 31927758327618367206289302020296005463331863499287179368579*x^40 + 2464456455393036467931853180993144281825912858646329643946*x^39 - 57282753586007415684000058798450052338020227846080090858*x^38 + 929876275164479656695807002583912108582639045389126219*x^37 - 12485235866954306130263093645484782044646194469136114*x^36 + 148281870295947612175805931658004559611524961794303*x^35 - 1612012770123305002535201993110614970838127144475*x^34 + 16366955764162863616656242426151528955066499910*x^33 - 157215788846403459357317335681167746154153432*x^32 + 1441419942233426111247575987038886388280618*x^31 - 12693969499550675470426730641925419720964*x^30 + 107882329734489301124300685555932947101*x^29 - 887961286084820649106528456808614966*x^28 + 7097839572742785961854456784079973*x^27 - 55219175251568110333400180827958*x^26 + 418829784219772636034142555868*x^25 - 3101527576670298970223292214*x^24 + 22449047472657242788847144*x^23 - 158968162795708635878326*x^22 + 1102172203562783475570*x^21 - 7486857251770649620*x^20 + 49854376824325898*x^19 - 325591388291103*x^18 + 2086404903275*x^17 - 13123677469*x^16 + 81060571*x^15 - 491839*x^14 + 2932*x^13 - 18*x^12 - x^11 - x^10 - x^9 - x^8 - x^7 - x^6 - x^5 - x^4 - x^3 - x^2 - x - 1
sage: p(1/16)
75047351740104879556819916558221890244724875220211010179789/11692013098647223345629478661730264157247460343808so the question is: which routine did compute that wrong msign value? It seems to me that the routine de_casteljau_intvec uses only integer/rational operations, thus the value sign which comes from it should be correct. What happens if one replaces assert(msign == sign) at line 736 by msign = sign?
mezzarobba
commented
1 year ago Thanks a lot for taking a look.
Assuming the polynomial corresponds to the list [-1, -1, -1, ...] with low-order coefficients first,
I don't think this is the case, but maybe I don't understand what you are saying. According to the docstring of interval_bernstein_polynomial_integer, [-1, -1, -1, ...] is a list of mantissas of coefficients in the Bernstein basis (over [0,1/2] in this case, I think), with a shared exponent (here 495). So, if I understand correctly, our polynomial is
sage: n = len(l) - 1
sage: p = sum(binomial(n, k)*x^k*(x-1/2)^(n-k)*2^k*c for k, c in enumerate(l))
sage: p
-11716129527125443095007855851527849572802892536021721259349096585967715*x^41 - 327788757504994699118480733259240434826210766607779868419405528115794615/2*x^40 + 169745677322830168727478208204350303533738043048201991338758432271900465*x^39 + 30518169743483218957876977480082746009756008610997750753460971801551745/2*x^38 + 5025171278930782510856480878530732019369088819640724830109264565062095/8*x^37 + 258933372047300549313593976285491438692158966750290650906536164390043/16*x^36 + 4726850796975746698971118054711938473220099879026915640063413930853/16*x^35 + 130816635848700882701978507628044109625209105936458857674799044525/32*x^34 + 11453425753446299922185045731467250907975932724862281780063364225/256*x^33 + 204132812557731234761308125916013158614992144749917830083722765/512*x^32 + 189100359293917214410366580540198688222765234772918230027229/64*x^31 + 2368095326449555871775767356451902884208637603999130483869/128*x^30 + 50729112844859655318824337161741558415609435556377355685/512*x^29 + 469226351788792872024291590501295510036914157665358745/1024*x^28 + 1888400776771596521643802209482867502871459384860365/1024*x^27 + 13307830442454493177247521040645993242814206766469/2048*x^26 + 660093059428341569870349018272882673957097426965/32768*x^25 + 3614898607392167732691699770203397140919868225/65536*x^24 + 4384628493144803218603771156545379030794375/32768*x^23 + 18890506975954197957791923766031587599575/65536*x^22 + 144787224840721134283558378567412905485/262144*x^21 + 494099670608773184823693672101217105/524288*x^20 + 751184946082778433029394265214475/524288*x^19 + 2035244741197169983115524200675/1048576*x^18 + 19645639416973503476396336925/8388608*x^17 + 42185885510208558346520409/16777216*x^16 + 10061133997315944088449/4194304*x^15 + 17032035756835107105/8388608*x^14 + 47414737873413365/33554432*x^13 + 593458526091145/67108864*x^12 - 34980575577831/67108864*x^11 + 4137487433937/134217728*x^10 - 6895812389895/4294967296*x^9 + 626892035445/8589934592*x^8 - 12292000695/4294967296*x^7 + 819466713/8589934592*x^6 - 91051857/34359738368*x^5 + 4101435/68719476736*x^4 - 71955/68719476736*x^3 + 1845/137438953472*x^2 - 123/1099511627776*x + 1/2199023255552
sage: p(1/16)
167734730327691663936749081528367790836259256068994501684866046481370335973/23384026197294446691258957323460528314494920687616so the question is: which routine did compute that wrong msign value?
split_for_target, but based on the sign fields of the rr_gap objects it was passed, and that's where it gets complicated! Ultimately, it could be one of the constructors of the Bernstein polynomial classes, or shrink_pb, or (maybe most likely given how this bug seems to occur only in extreme cases) more_bits.
I was hoping that you understood the algorithm well enough to make a plausible guess... :-)
It seems to me that the routine
de_casteljau_intvecuses only integer/rational operations, thus the valuesignwhich comes from it should be correct. What happens if one replacesassert(msign == sign)at line 736 bymsign = sign?
Then I get
sage: ibp.de_casteljau(ctx, 1/16)
(<IBP: ((-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, -2, -2, -1, -2, -1, -2, -1, -2, -1, -2, -1, -2, -2, -3, -2, -6, 19, -108, 460, -1561, 2083, 22624, -245343, 1313184, -2002406, -42379996, 569661788, -4354890636, 23877560352) + [0 .. 4)) * 2^495 over [0 .. 1/32]; lsign -1>,
<IBP: ((23877560352, 447364325174, 8331176857701, 154266207126873, 2841025539290806, 52049654442092246, 948807530809418581, 17211465959776074626, 310730026553772917470, 5583489482455840182304, 99862602115934902277802, 1777783203337883935294765, 31500980456632009102118537, 555540389709548050733964486, 9750260899743195116964849484, 170283865194501927358041835497, 2958820910074760843578864154706, 51140491814494121311022921130344, 879032346066839601043638138571069, 15021310154293731935170089565190618, 255103097093067960135980438051905274, 4303662807719911273781493233899063074, 72085195531703401341742558772819474090, 1198015785831024629866599863258633509621, 19740220642975151903183233720579633724597, 322184445471886005108584545983332110168257, 5202534317300063066059410863923155869868156, 82995191471425326555364656946313503465801362, 1305632571746143983709096408111859213422207807, 20206530015097254084832421548391959936353293933, 306699164482870665709538267526188846822461523780, 4546253874722966609050264977317618301943981603866, 65426461468145659494659210542267982117661506417207, 906282427445111119788053476486919539443728242095895, 11922390007792935437590310274647322489089554794903150, 145619389914953066764086839837869848193603051076670907, 1580953356712087172113097181357376120571400507801578596, 13730043066567805006089835614502160881546198268020372111, 60579410049484097241141959422504281415028595670106126133, -736822160535660204069440378527845240448857176001037695156, -16513300040392919339506510770134110188484992360148883082431, 214703574267988798662235357982294318934218074726925561209802) + [0 .. 4)) * 2^495 over [1/32 .. 1/2]>,
True)And the call to real_roots returns 12 isolating intervals
[((0, 0), 11),
((37465753076616074685581795154893275042071217153399687052041634117786117404748180825691/62165404551223330269422781018352605012557018849668464680057997111644937126566671941632,
38364931150454860478035762194091702820461747831488769767720688614395449715913347605715559/63657374260452690195888927762793067532858387302060507832379389042324415617604272068231168),
1),
((268435103/268435456, 134217965/134217728), 11),
((109396387/67108864, 328189161/134217728), 1),
((5789/2048, 12405/4096), 11),
((817424974259450517108044311062521326211325218106685357371502090135192085255/226156424291633194186662080095093570025917938800079226639565593765455331328,
51089060891215657320094234704185588586169705961560086962110752375649875785/14134776518227074636666380005943348126619871175004951664972849610340958208),
1),
((27606985387162255149739023427370689427422170705051945374080728870453965/6901746346790563787434755862277025452451108972170386555162524223799296,
110427941548649020598956093896341333562562541192944547040497565473380985/27606985387162255149739023449108101809804435888681546220650096895197184),
11),
((4135/256, 9097/512), 1),
((2481/128, 827/32), 1),
((356437/4096, 827/8), 1),
((4135/16, 9097/32), 1),
((2481/8, 827/2), 1)]while there are only 11 real roots. The interval is (356437/4096, 827/8) does not contain any root; the others seem to be correct.
zimmermann6
commented
1 year ago I was hoping that you understood the algorithm sorry I don't know de Casteljau's algorithm. What I know well is Uspensky's algorithm, which I have implemented for example in CADO-NFS (directory utils, file usp.c). This is a 600-line file (excluding the main() part), compared to 4650 lines for
real_roots.pyx.
See also #26955.
CC: @videlec @dimpase @orlitzky @sagetrac-cwitty
Component: algebra
Issue created by migration from https://trac.sagemath.org/ticket/20390