AssertionError in real_roots.real_roots()

[mezzarobba]

sage: Pol.<n> = QQ[]
sage: pol = (-28524712624242039055315336000413967778733418480863118977869865132075390079059109742198859964232127005/1451626239969468099340993140755597642170368*n^85 + 106129053170798145292157267201895276884016106806128019074207858099654340989739989819086948789647680411711/6532318079862606447034469133400189389766656*n^84 - 48187461003922327740551674859848326512999070262843108330739019140511010180683826810454677328853749567659529/9798477119793909670551703700100284084649984*n^83 + 42762289309042999830570418433539358762819280522472412204583485774468713937906840315002959365007403789113983333/58790862718763458023310222200601704507899904*n^82 - 2618262455746996034793717082427769260012747160045137879009436993298900640475242462655902215693284389473549525853/39193908479175638682206814800401136338599936*n^81 + 21184278324043255960865454202595855100443658574318565625522527050972139625986684343158454636555343817236600979557/4899238559896954835275851850050142042324992*n^80 - 43366900101688156340338924312027109375611629318571702313036703102942950280704023013723734060949072250338097966263/204134939995706451469827160418755918430208*n^79 + 13572135639646239474763103960005344090575994043190066312356381948195830305343420437915768559698229005465629179977771/1633079519965651611758617283350047347441664*n^78 - 2621950604238651742641700807396667564130622602387996727209846221545805994994062896408231034527495711169978850180116385/9798477119793909670551703700100284084649984*n^77 + 213541759148758791831119082131136221510572110310516167585008032697272377714912315344199105604174026356568352724364112571/29395431359381729011655111100300852253949952*n^76 - 51857932671475150757351386835473140539389211185900765104711839807127640289329023366920218966169252835344633871733079851/306202409993559677204740740628133877645312*n^75 + 33702068710711968826608440662977108099998066632506421948624544601354866434420790701794657988262260923208242812844178431373/9798477119793909670551703700100284084649984*n^74 - 1809429934990413763624860745414371427044930547717750702590658010864143409159073381528696539124438620399757738683135469876707/29395431359381729011655111100300852253949952*n^73 + 2400324019270149165515899070382263844082746063647750106754802717767023339894255870442874721441739829329917440319024965982791/2449619279948477417637925925025071021162496*n^72 - 8563390247528340528803408428350155033736370396166974037074003432558960979000532892313091186161660022197220844642101865296931/612404819987119354409481481256267755290624*n^71 + 441074701733100682032818138203764777734702675472074278029850710261302021681370301127998740635090960193435126788416296708873997/2449619279948477417637925925025071021162496*n^70 - 41228517545242808616083101310512767978764159118535358927545416990880852361293634185679618332142179596636676160119205219464875241/19596954239587819341103407400200568169299968*n^69 + 219597858716984702598026480348555525565465597400033851995577594297145515979345111884385202040155124918058177443774789295204315343/9798477119793909670551703700100284084649984*n^68 - 3212472730211623163653269657051594893068965842200301918670128210554692406592337278919297930353479378182025697704072123206071433569/14697715679690864505827555550150426126974976*n^67 + 19190664269507997348121536288459319782965183142745351113470155438684896178634719723368680838643642647036315341062757807537906154017/9798477119793909670551703700100284084649984*n^66 - 316993213555679391726210288889440824656684302928783359105489910754134411008475416023517691831185162336039875428364650552399081287643/19596954239587819341103407400200568169299968*n^65 + 453700925885694764994609854627854639453901855081341311952679430525187241742746023627902014882378414521161140252336619806026588066259/3674428919922716126456888887537606531743744*n^64 - 44565821448374992726481694737317462894037084667361374427083302044226433288822386946993072592760771028569874689570492426160471274853/51033734998926612867456790104688979607552*n^63 + 3511791416828609821132523477752726791348888447614622129379516952226997340972812157324216161592305714494138617534853515688835690265891/612404819987119354409481481256267755290624*n^62 - 343243069912338715119597661756552935584869354986864776842628621594257110669008625530358180165633959774582416557434506894159576800813631/9798477119793909670551703700100284084649984*n^61 + 61058231989649998211427458286138470940800166799620647173200436577482489844896721411902702381559157615799572644306120571256347297715033/306202409993559677204740740628133877645312*n^60 - 2594857585618410103811665572559679461480204576934813372672198100081875578449829536676827679192156087507735621503592686015318319569054463/2449619279948477417637925925025071021162496*n^59 + 38643205696871931088560343392723663793721450517311357581197348085972006850651379397908469890727322014156754805531372384917934258884082233/7348857839845432252913777775075213063487488*n^58 - 26587101559599614179573549328248440296070700391106639612726146191189142418008441481898231843143150160848122488375760942488533270381965789/1088719679977101074505744855566698231627776*n^57 + 260172203220611089845780395010476873024675206853025482743895281428542524137915961075729990235009825605734086355583701542828515755744654977/2449619279948477417637925925025071021162496*n^56 - 795465537263071015534002209903430699836574801760918154756784997483043711478399066381468461072138847229375269678063461818223616667147525357/1837214459961358063228444443768803265871872*n^55 + 4055881045582182715054862248335323340571359119514227264959526835412120843088363219481502169719358570215952452846674295573109383381966568399/2449619279948477417637925925025071021162496*n^54 - 232938877013346540679103993878875270665268548639082935555531923233007720357876752038759810387582440553127791742725418153751269272823934452071/39193908479175638682206814800401136338599936*n^53 + 130881972954338736275105980232446932551270675717252216944343334423633686100353250981832325881704829585502714817879812312170071547665564628535/6532318079862606447034469133400189389766656*n^52 - 207295247430367670137381172254411163734326309667243882423819320158131653043316826143196215768336414753194112028756007010334377063964668452811/3266159039931303223517234566700094694883328*n^51 + 3703293710299778086428298075091197811183582634673428689744302209457330235368908199547161021358591598816601332430134512885389938520810937253071/19596954239587819341103407400200568169299968*n^50 - 62203375686471580123784123343410176154057852085289755393064632850731108640967287392595505024286394080915958966085089070168189979119216390880575/117581725437526916046620444401203409015799808*n^49 + 6823163443776981282501002797761805215178988579986714299579334365231130497515469507434830327694924450864403662544894491405592164260113267066031/4899238559896954835275851850050142042324992*n^48 - 2111701208091545604778113051968487213334613488901973060099644573102137636009663456749637159895668136311584812046403255236901911618844144835701/612404819987119354409481481256267755290624*n^47 + 118016656801039867033437791270830297232471341070477873874305984917489400400671377828386021027222101116836987985477338640752126685723502004381085/14697715679690864505827555550150426126974976*n^46 - 86150610675202429982375394372258682929437537147248110966706766731353988695097803881031062025306332157122668012800535492031845329073886020908841/4899238559896954835275851850050142042324992*n^45 + 354813086964037595886750116765126589484373671214105887912883366693443435758600432415651151159926861481543862242245874345000131098589138093684603/9798477119793909670551703700100284084649984*n^44 - 171722909472605390672928310434419590310342927502145772447297828137268659373758006392455417810005950974519882469235840215205074324446403331952599/2449619279948477417637925925025071021162496*n^43 + 1249740677122157841340510720173938982496818680305902066804258835378301132201858801777222065667985704579900970800395319700749127035347303069599779/9798477119793909670551703700100284084649984*n^42 - 133506603371010833132077589933388900916093368576226847867889485617428438674242642766423717090778912114772694554322256935224588309192071225725789/612404819987119354409481481256267755290624*n^41 + 642816318522377843400706703271031804131615785949614263831252339816454051934802079935946547915154891076175998222986212756569771409923774940691045/1837214459961358063228444443768803265871872*n^40 - 20170491501078288581134446640983080002349678089800195837555131264631574033178569040849437674526748426112354165792856744137636796355979355131683/38275301249194959650592592578516734705664*n^39 + 227923784857334381888441976923957413217948530851541376848317379652083235826411690844774576956250569916054003213483419513131311478796934373199691/306202409993559677204740740628133877645312*n^38 - 7240265005176416379053889690882759632762746706289971842422529713015257851057336027687211334756912305077855278328722677019778552164065532729223357/7348857839845432252913777775075213063487488*n^37 + 996882558724619246002003227904557978410704119122718995495379097845567797133963152619588565361410391334636694693278109825629483337968717512808997/816539759982825805879308641675023673720832*n^36 - 1733004085692285900260265585477331156516157016949479940952683307572277609930404252405999890118847400808302285329191868347409967982896415409799385/1224809639974238708818962962512535510581248*n^35 + 117265588340919225361893400152616559392585635699286767863088277529094611306041370243184355655344655945946210297134862660975180909525232420210427/76550602498389919301185185157033469411328*n^34 - 29611263800068275566230203638186900257373612279283812642972792584103598272398962493951716081283975603231524571102070432046643344132637383684229/19137650624597479825296296289258367352832*n^33 + 6965531536301560074645713272094785795253228012325824333420443411433782816880232631896315019052272449901841991391455796806046125376187806894307/4784412656149369956324074072314591838208*n^32 - 9142970171987438597615997005152146581775186485881951321344462268857027158568665487333696909542767828150272682197911086792122026821513583910947/7176618984224054934486111108471887757312*n^31 + 12893524121483164761537310026969146503292685575182290476950350242914657414391293053309633978507935082501960432131112445708436256817316489239/12459407958722317594593942896652582912*n^30 - 6467354526540620489642420756923455173462978731025907853169965733676787089535857410729802244748839753795163548108696074172235639671537100759/8306271972481545063062628597768388608*n^29 + 2248198215893424811476813300390874102843931372799615421331899101546483463602279740457249425339371795735262078632116839134988360939715223723/4153135986240772531531314298884194304*n^28 - 80020564367455378146684948091148376360562062241125870664822725933825169415158183384353136876972744429856426569768858099840647877705690007/230729777013376251751739683271344128*n^27 + 490585126108898559374431784920654270112603802910220279530374628644096434864728916969337852643644926171227475356337479142258234339070031/2403435177222669289080621700743168*n^26 - 66073504438022085528467639050728145775784269135333710822304366556644940625846993430955633327788851038115280922413767720008185805414403/600858794305667322270155425185792*n^25 + 901348438337341035274724698441799840705173584801487754149095521681387839822659827057552501530610206533328065459726466445471298098837/16690522064046314507504317366272*n^24 - 66929813811935451825454410566167178573240768796194075400923295970402593210981550939143630462172584295448457604630405070600696853719/2781753677341052417917386227712*n^23 + 1681295438300968627769053952625967798309111365872391720500275328597486352140241954859373243298123154434989365676432582334283611893/173859604833815776119836639232*n^22 - 5607051673965509789953274338155832185101485475533415602600000257942892945334509734473787755910703785087486186306170609745870133/1609811155868664593702191104*n^21 + 3590637386720442005402477768988546933545272496775843175260358165964885850288875947727556104978246161793242426123216758011053447/3219622311737329187404382208*n^20 - 506058176618441495038176401204534369182994150081466863944671455666741210467793304506009389747049444131711956835710759233397941/1609811155868664593702191104*n^19 + 861538898207043961047669466632570506687338219508418340094737780497618128612605571852860272753799185426692758249207969703345/11179244137976837456265216*n^18 - 15072584106783073774904665954367274557501302637483181456155309205953606894402191641741301955971958200170822616793794577585/931603678164736454688768*n^17 + 663804358691226265591073685749685611746140406890360067849368534141655460689054294361547637680734933893196467848359580143/232900919541184113672192*n^16 - 5301240044439525252134982074833108095455875778905014197655556283070369216994066731156090716107200104890178264317038977/12938939974510228537344*n^15 + 6220084452552043967635060245197312388501290170430731456055726736507132343414103218664244762853792885616026198198505/134780624734481547264*n^14 - 514810176078266640412018821895444140385603064721798916602413778330571818858202473017015450291035298094664430873725/134780624734481547264*n^13 + 1546468651746133868702294176342894895267439358913165081888164869199810853512873675122324919432312199003142551027/7487812485248974848*n^12 - 758920559048027343748732545186719290162745430942201696641161166096623883661324931066893164360154890187720241/138663194171277312*n^11)
sage: from sage.rings.polynomial.real_roots import real_roots
sage: real_roots(pol, (QQ(0), QQ(827)))
...

AssertionError:

See also #26955.

CC: @videlec @dimpase @orlitzky @sagetrac-cwitty

Component: algebra

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

[jdemeyer]

comment:1

Traceback please

[mezzarobba]

comment:2

/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]

comment:3

bug is still present in 9.2.beta2

[orlitzky]

comment:5

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.

[mkoeppe]

comment:8

Moving to 9.4, as 9.3 has been released.

[mezzarobba]

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)]

[mezzarobba]

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]

who is the author of that code ? Where is it located ? Is there a simpler example (I cannot copy&paste the long input) ?

[mezzarobba]

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]

No simpler example, unfortunately, but here is this one as a .sage file

[mezzarobba]

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]

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]

@mezzarobba please can you create a Sage file with the inputs to the de_casteljau routine that exhibit the issue?

[mezzarobba]

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]

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/11692013098647223345629478661730264157247460343808

[zimmermann6]

so 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]

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/23384026197294446691258957323460528314494920687616

so 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_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?

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]

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.