Opened 4 years ago
Closed 3 years ago
#21579 closed defect (wontfix)
Errors calculating characteristic polynomials of rational matrices
Reported by:  bober  Owned by:  

Priority:  major  Milestone:  sageduplicate/invalid/wontfix 
Component:  linear algebra  Keywords:  
Cc:  cpernet  Merged in:  
Authors:  Reviewers:  
Report Upstream:  N/A  Work issues:  
Branch:  u/cpernet/errors_calculating_characteristic_polynomials_of_rational_matrices (Commits, GitHub, GitLab)  Commit:  08012ea976866316645c2c63cf834d5317a50301 
Dependencies:  #24214  Stopgaps: 
Description (last modified by )
I'm leaving the description below because that's where I first ran into problems, but this seems like a much more basic issue. I suspect some sort of memory bug, because I'm getting different results based on the order in which I execute functions.
This could just be a problem with polynomials. I haven't quite diagnosed it yet.
A = Matrix([[ 10/7, 3/4, 1, 8/7, 3, 1/7, 1, 1, 7/4, 2, 3, 1/2, 8/3, 3/8, 10/7, 3, 1/2, 4/7, 9/8, 4], [ 3/2, 3/5, 3/10, 4/3, 2/7, 3/2, 1, 3, 2, 0, 4/3, 2, 4/9, 1/2, 1/5, 1, 7/8, 1, 5/2, 10/9], [ 9/10, 7, 3, 0, 3/2, 1/3, 5, 8/3, 6/7, 0, 7/4, 4/3, 1/2, 1/6, 1/4, 3/7, 6/5, 1/9, 1/3, 8/5], [ 8, 1/2, 8/7, 1, 3/8, 7, 10/7, 1, 4/3, 0, 0, 4, 3, 1/4, 5, 1, 1, 0, 2/7, 10/3], [ 8/3, 5, 2/3, 9/4, 9/8, 1/3, 2/3, 1, 8, 4/3, 3/4, 5/4, 2/3, 1/2, 5, 1/3, 10, 1/2, 5/3, 9/2], [10/9, 3, 7, 1/2, 5/9, 0, 3, 1, 9/4, 3/8, 10, 0, 6, 1/6, 2, 9/8, 2/3, 6, 0, 0], [ 2/5, 6/5, 0, 0, 6/7, 9/10, 2/7, 1/10, 1/3, 7/2, 1/3, 3/2, 7/2, 3/2, 1/3, 0, 7/2, 2/9, 2/3, 1/8], [ 1/3, 3, 10/3, 1, 1/2, 3/8, 5, 3/2, 1/3, 0, 1/2, 4/9, 7/5, 0, 3/2, 9/10, 1, 1/2, 1, 9/10], [ 1, 3/8, 1, 2/5, 1/3, 1, 3/2, 4, 1, 4, 1, 1, 9/10, 0, 1, 1, 7/6, 0, 10/9, 0], [ 6, 0, 8/5, 9/7, 7, 1/3, 8/9, 8, 7/9, 1/2, 1, 1, 1, 2/7, 8/3, 1/3, 4/5, 2/3, 5, 1/7], [ 5/3, 1/3, 9/7, 2, 4/5, 5/7, 1/4, 5/8, 1, 5/7, 3/2, 3/2, 2/5, 0, 1/2, 1/4, 0, 3, 1, 0], [ 7/8, 1, 6, 9/7, 3/5, 1/2, 2, 7/2, 1, 6, 7/2, 1/9, 1, 7/5, 5, 5, 2/7, 5/4, 3, 3], [ 7/3, 1/5, 2, 8, 1/2, 2, 0, 2/7, 1/4, 0, 0, 4/7, 1/3, 1, 5, 3/2, 10, 3/10, 1/3, 3/4], [ 8/5, 0, 4/5, 0, 1/4, 8/5, 4/9, 5, 7/9, 1/9, 7/8, 5/8, 1/4, 10, 1, 8/9, 1/2, 1/6, 8/7, 1/7], [ 2/5, 5/4, 2, 1/2, 3/5, 5/4, 10/9, 1/2, 3/5, 9/7, 5/6, 3/7, 7/6, 4/5, 9/5, 7/3, 1/3, 1, 5, 2/5], [ 1/2, 2, 0, 3, 5/6, 0, 3, 1/3, 1/2, 1, 2/9, 5/7, 1, 8, 1/2, 10, 5, 1/2, 1, 5], [ 8/3, 2, 2, 4/3, 9/2, 4/5, 10/7, 1/2, 9, 1/3, 10/7, 4, 7, 1/10, 10/9, 7/6, 2/7, 2/7, 0, 8], [ 2/3, 1/2, 4, 3/5, 2/3, 2/7, 9/10, 1/4, 5/3, 6, 1/3, 1, 5/4, 4/3, 6, 8/7, 2, 1, 5/6, 1/3], [10/3, 5, 0, 2/5, 4/3, 1/2, 1/2, 9/4, 6, 1/5, 1/9, 2, 6/5, 2, 4, 3/4, 1, 1/4, 3/8, 1/2], [ 7/6, 7/10, 0, 1/5, 3/4, 3/5, 2, 3/2, 1, 6/7, 5/7, 0, 9/5, 1/2, 1, 2, 7, 0, 4/5, 9/2]]) k = int(sys.argv[1]) if k == 0: sage.matrix.matrix_integer_dense.USE_LINBOX_POLY = True # True is the default A._clear_cache() f = A.charpoly() # should be using linbox sage.matrix.matrix_integer_dense.USE_LINBOX_POLY = False # Setting this False should force # integer charpoly to use flint when # algorithm=='linbox'. (Rational # matrices don't have a flint option.) A._clear_cache() g = A.charpoly() # should use flint, I think A._clear_cache() h = A.charpoly(algorithm='generic') elif k == 1: sage.matrix.matrix_integer_dense.USE_LINBOX_POLY = False A._clear_cache() g = A.charpoly() # should use flint, I think sage.matrix.matrix_integer_dense.USE_LINBOX_POLY = True A._clear_cache() f = A.charpoly() # should be using linbox A._clear_cache() h = A.charpoly(algorithm='generic') if k == 0 or k == 1: print 'flint == linbox:', f == g print 'flint == generic:', g == h print 'linbox == generic:', f == h print 'linbox (maybe) correct? ', f(A) == 0 print 'flint (maybe) correct?', g(A) == 0 print 'generic (maybe) correct?', h(A) == 0 sage.matrix.matrix_integer_dense.USE_LINBOX_POLY = True if k in [0,1,2]: # let's try computing the charpoly the way that # B.charpoly() does it to see what goes wrong. B, denom = A._clear_denom() B._clear_cache() f1 = B.charpoly(algorithm='linbox') B._clear_cache() g1 = B.charpoly(algorithm='flint') print 'flint == linbox', f1 == g1 print 'linbox (maybe) correct?', f1(B) == 0 print 'flint (maybe) correct?', g1(B) == 0 x = f1.parent().gen() F = f1(x * denom) * (1 / (denom**f1.degree())) x = g1.parent().gen() G = g1(x * denom) * (1 / (denom**g1.degree())) print 'flint == linbox', F == G if x == 0 or x == 1: print 'linbox == linbox', f == F print 'flint == flint?', g == G print '2nd linbox (maybe) correct?', F(A) == 0 print '2nd flint (maybe) correct?', G(A) == 0 if k == 3: # let's try computing the charpoly the way that # B.charpoly() does it to see what goes wrong. B, denom = A._clear_denom() f1 = B.charpoly(algorithm='linbox') print 'linbox (maybe) correct?', f1(B) == 0 x = f1.parent().gen() F = f1(x * denom) * (1 / (denom**f1.degree())) print '2nd linbox (maybe) correct?', F(A) == 0 if k == 4: # let's try computing the charpoly the way that # B.charpoly() does it to see what goes wrong. B, denom = A._clear_denom() g1 = B.charpoly(algorithm='flint') print 'flint (maybe) correct?', g1(B) == 0 x = g1.parent().gen() G = g1(x * denom) * (1 / (denom**g1.degree())) print '2nd flint (maybe) correct?', G(A) == 0
On the most recent version of the develop branch, here's output that I get consistently:
jb12407@lmfdb5:~/sagebug$ sage bug.sage 0 flint == linbox: False flint == generic: False linbox == generic: True linbox (maybe) correct? True flint (maybe) correct? False generic (maybe) correct? True flint == linbox True linbox (maybe) correct? True flint (maybe) correct? True flint == linbox True 2nd linbox (maybe) correct? False 2nd flint (maybe) correct? False jb12407@lmfdb5:~/sagebug$ sage bug.sage 1 flint == linbox: True flint == generic: False linbox == generic: False linbox (maybe) correct? False flint (maybe) correct? False generic (maybe) correct? True flint == linbox True linbox (maybe) correct? True flint (maybe) correct? True flint == linbox True 2nd linbox (maybe) correct? False 2nd flint (maybe) correct? False jb12407@lmfdb5:~/sagebug$ sage bug.sage 2 flint == linbox True linbox (maybe) correct? True flint (maybe) correct? True flint == linbox True 2nd linbox (maybe) correct? False 2nd flint (maybe) correct? False jb12407@lmfdb5:~/sagebug$ sage bug.sage 3 linbox (maybe) correct? True 2nd linbox (maybe) correct? True jb12407@lmfdb5:~/sagebug$ sage bug.sage 4 flint (maybe) correct? True 2nd flint (maybe) correct? False
All those False
s ought to be True
. Option 3 does get the right answer, but I was getting (less directly repeatable) bugs just by calling charpoly() repeatedly on the same (rational) matrix, which should have been using only linbox.
In another (nonrepeated, interactive) test, I wound up in a situation that was something like
sage: f = A.charpoly(algorithm='flint') sage: g = A.charpoly(algorithm='linbox') sage: f == g True sage: f(A) == 0 True sage: g(A) == 0 False sage: print f ... sage: print g ... # yes, the certainly looked the same!
(I might have the True and False mixed up, though.)

Old Description:
I ran the following code for a bit
@parallel def test(n): M = ModularSymbols(3633, 2, 1) S = M.cuspidal_subspace().new_subspace() f = S.hecke_matrix(2).characteristic_polynomial() return S.dimension(), f for result in test(range(5000)): print result
and line 530 of the output is (((15,), {}), (171, 1))
.
line 621 is (((22,), {}), (171, 1))
line 663 is (((156,), {}), (171, 1))
Line 757 starts (((281,), {}), (171, 2167647436431476168725215705905073292546268...
and has a bit over 14 million characters. I haven't checked the degree of the polynomial yet.
(Most lines look something like (((28,), {}), (171, x^171  x^170 ...
and have a bit less than 6800 characters.)
That's a random but pretty high failure rate, and also suggests that when whatever might go wrong does go wrong, the computation takes a lot longer to run than normal.
I may follow up with more tests. I have no idea what is going wrong or if there is something special about the level. I only tried level 3633 because in a separate computation I saw an error there which I was trying to reproduce. In that case I was actually getting a degree 170 polynomial somehow (and I did not separately record the reported dimension of the space), so I may leave this running for a while to see if that ever happens.
Attachments (2)
Change History (28)
comment:1 Changed 4 years ago by
comment:2 Changed 4 years ago by
 Component changed from modular forms to linear algebra
 Description modified (diff)
 Summary changed from Modular symbols computations are inconsistent and occasionally wrong to Errors calculating characteristic polynomials of rational matrices
comment:3 Changed 4 years ago by
It would be good if you could reduce the issue further :)
comment:4 Changed 4 years ago by
sigh...
I think that while trying to diagnose one bug I found another, which is now #21596.
comment:5 Changed 4 years ago by
 Cc cpernet added
comment:6 Changed 4 years ago by
This might be related to #15535 (and related tickets therein).
What happens is that charpoly over ZZ uses chinese remaindering of charpoly's over Z/p_iZ. The termination condition of the chinese remaindering is that the reconstructed charpoly stabilizes when 3 more moduli are added. The infinite loop presumably occurs when one charpoly mod p_i is erroneous, thus preventing the stabilization to happen. I've been torture testing the charpoly mod p without success for the moment. I'll try to use you test case matrix to trace what happens.
About the variance in the timings: this is a Las Vegas probabilistic algorithm: it may fail but then raises failure, hence restart a computation, hence the result is supposed to be always correct bu the computation time may vary. Still this large variance is really strange.
comment:7 Changed 4 years ago by
(I tried to send this via email this morning, but that functionality is either broken or temporarily down or, perhaps, a figment of my imagination.)
I did see the "running out of primes" message the first time I ran into problems. That was on an "ancient" version of Sage, though. I'm not seeing that message now, unless it has been redirected to stderr and lost.
My current test has still not produced an error running under the current Sage beta, but it does seem as if some processes might run forever. If I don't eventually see an error, perhaps I'm going to have to look carefully at Sage's charpoly for rational matrices to see if it is doing anything questionable, since I did see an error there. (I doubt that some sort of conversion bug would only occur 1 out of 5000 times, though.) It does fail fairly quickly in Sage 7.3, preLinbox upgrade.
I took a quick look at your paper http://arxiv.org/abs/cs/0501074, which I guess is roughly what Linbox does. About the algorithm used: if the charpoly mod p is failing more often than you expect, a simple hack to avoid the infinite loop over ZZ would be to use the Hadamard bound to detect problems with the early (never) abort.
Also, is the charpoly() in Linbox really not proven correct? If you do what the paper describes, it seems it isn't, but I'd be surprised to find a genuine example where it gives the wrong answer that's not a bug. That sort of thing is generally fine by me if it really works as advertized, but I think it is frowned upon by Sage, which wants everything to be a "Theorem".
(I've been doing the same sort of thing recently, building polynomials over ZZ by computing them mod p and waiting for them to stabilize, and I've been wondering whether I'll ever actually be able to go back and prove them correct.)
If it helps, you can post or send me some pure C++ Linbox code, and I'll run it to test. (I've never used Linbox directly before.)
comment:8 Changed 4 years ago by
I think I finally found the bug(s). The one causing the neverending computation is in Givaro's random iterator choice for a seed. It used the 6 digits of the microseconds of current time. Hence once in a blue moon (with proba 10^{6), the seed would be 0 and the congruential random generator would always output 0, which would cause an infinite loop in the minpoly, searching for a nonzero vector to iterate the Krylov basis. }
In this process, I also found what could also be a bug in the LUKrylov charpoly. I attach the 2 patches to this ticket: one to be applied on givaro, the other one on fflasffpack. Could you let me know if it fixes all your bugs/infinite loop problem.
comment:9 Changed 4 years ago by
 Branch set to u/cpernet/errors_calculating_characteristic_polynomials_of_rational_matrices
comment:10 Changed 4 years ago by
 Commit set to 0e8bc88f63a0988fc0047c1bc3c84b588b05290a
The 2 patches are now in the ticket's branch.
New commits:
0e8bc88  add givaro and fflasffpack patches attached to ticket #21579

comment:11 followup: ↓ 15 Changed 4 years ago by
I've been testing for a while now and see no noticeable difference. There are 12 processes that have been running for more than an hour now.
I hope I'm testing correctly. I checked out the branch and then forced a rebuild of the two packages. Then I realized that might not be enough, so I updated the packageversion.txt for each of them (which your patch should probably do) and ran make build.
The random seed bug is a funny one. I haven't figured out what the charpoly patch is doing.
I wish I had saved the output the one time I got a wrong answer on the current development branch. If I get a wrong answer again, it should be saved. I wonder if the result will be a polynomial which has the correct reduction mod p for lots of p but is wrong somewhere and for which the CRT will miraculously stabilize.
comment:12 Changed 4 years ago by
Actually, there is a difference. It does seem to be always finishing now, but it still sometimes takes ridiculously long. Actually, ignoring the ones that ran forever before, the distribution of timings is about the same.
comment:13 Changed 4 years ago by
After letting the test run to completion 100000 times, I finally got a wrong answer again. The correct answer and the wrong answer differ only in the constant term. With the correct answer being f, and the wrong answer being g, I get
sage: f  g 151012103848425355244200954703540426043866904165204055830261178949029704179162112427542664943806786298858476093532538835074924828446540823256367823339218258609681625753661508234861960817139615529399832981296111900679567024247219254268104611895097155970778236782887697656762154816061591814647082042340089443747291763928372681213079578424656507828771748777750852982258479888495869657460409348964325583114546973737498080511951329032911155629250820799309534264678106380858222598586826132396800534297679417898365363687982919973098676243412382505418978680668442099909762634433088449097140756248589241278295016885267708001337470599645081698363942965176093920140998988309374034976444892009811758973426469656640517612333957869075629697472441972189753967972750212471081376890624479658069564505036282529179135314797769829482780785188541453131105257999316397948389490215366051042230181219888105765051958866258221137838619488085694866455742489857058716310280616273213285476379508677332371296020658867691461563709440790778840690519721186167443771058589072654152844625174186571503521067198235222815898281363463492162550437482326665450066304688121300443342676890827316634527457180578708847236138890147832449069212579845664125867630345819952085579949872263779625303575099423861420547519107530023922093527691981952824468598382154585250245375315585543414389978966017515442661974659972069518242643099845111128042081400000970203042255427235699986741046359370297735867675150765488038014558644477415356490752764824701203628534645224670777179182750201270498922201705310852108254575237091719841041565378392871836672668185130693637807064973171425532512837006227167559335323068665370581813047716584004816839780780960796454333944491643551958879552213915396578724553718876517985064637899902165323417711016024737765666077679740021880911972907069038019105632520554466212908430287361402774270457049773592966623777465415630965326144463370205562537302258219907811599149769018975454937337470118499694397417996996793255845351935822462305292835717260066984966995095284742142675760478095595955303563397501653503811337582774010336310591161025026071916964766548895503395690162923675437405050957840442976595709886399994537045799288783349672622616060981316691085665035970016969244413193376889078634084996319260895896065892088462526538163831611354025202314551952799228805662030953907364050606529627591008411535902635924153255682592597222741558810265423802000799368481733033178008937529779901311327951965057427776889891335528073171132536428566058277469502207610991305800271349804773157529312296863578237698631687562830954703653493867705233537044684301823920336293430938536558711784539126939841241784406694213376143587334812926722610378920618688919469457196321821565456484443531968149943745115938277102261047615104875351799786733744292494807293742761354118284382898289012583171259164384012677917047674083644685911717643199987025562791925507393639022448335226983468938886302274678797912661136810819248611930125512663475201464873166136850888931341829013631445011342465235306106492081497150421961902438922399779253466938390253041808421848908966144097213141920889571560337205119907190343697616554437492788409961571209346081107495663816341735296635376065934630408876652626371492352582326904728077848256028260564646026799006001417139296979385513878199935667261985468913676293901703803768939425359580444756822981507606020306753977496086878692933752436987168178347148832881996614083930419165658242804898326875894619835
comment:14 Changed 4 years ago by
Here is the factorization of f  g:
1 * 5 * 2102873 * 2105993 * 2107607 * 2108699 * 2111471 * 2119483 * 2121121 * 2126767 * 2126923 * 2129749 * 2130853 * 2135141 * 2145937 * 2151241 * 2156303 * 2156897 * 2157511 * 2161927 * 2164177 * 2172613 * 2174197 * 2174617 * 2176501 * 2183453 * 2189599 * 2190521 * 2191949 * 2192327 * 2194811 * 2195801 * 2197333 * 2199143 * 2206387 * 2207437 * 2214713 * 2219177 * 2221507 * 2226131 * 2227109 * 2229767 * 2231011 * 2239319 * 2249537 * 2251279 * 2251517 * 2261807 * 2262017 * 2268977 * 2273651 * 2275027 * 2277809 * 2282141 * 2282173 * 2282699 * 2286883 * 2287381 * 2292181 * 2303663 * 2307821 * 2319613 * 2324233 * 2325871 * 2329597 * 2342437 * 2344939 * 2345209 * 2347129 * 2347559 * 2357231 * 2366057 * 2371073 * 2373431 * 2375917 * 2376079 * 2380607 * 2382337 * 2382937 * 2387039 * 2389643 * 2390389 * 2396309 * 2403887 * 2406227 * 2414053 * 2414371 * 2414717 * 2415839 * 2425043 * 2431453 * 2431523 * 2431691 * 2434207 * 2437147 * 2438587 * 2455889 * 2462689 * 2469919 * 2481697 * 2482177 * 2485579 * 2487883 * 2490973 * 2494433 * 2511541 * 2513519 * 2516351 * 2517677 * 2518823 * 2519747 * 2526361 * 2531099 * 2534993 * 2536477 * 2546749 * 2547029 * 2548391 * 2550133 * 2557099 * 2558909 * 2565191 * 2568457 * 2568983 * 2571017 * 2572589 * 2575219 * 2579051 * 2583661 * 2591513 * 2600267 * 2601089 * 2602331 * 2607359 * 2611157 * 2612143 * 2615177 * 2622013 * 2624059 * 2626049 * 2628169 * 2628191 * 2628817 * 2637647 * 2642771 * 2650931 * 2657119 * 2657939 * 2661419 * 2662811 * 2662981 * 2666633 * 2667131 * 2668793 * 2668993 * 2676637 * 2679277 * 2681191 * 2683477 * 2686451 * 2689711 * 2692367 * 2693707 * 2698387 * 2703137 * 2705447 * 2706059 * 2708711 * 2711327 * 2714783 * 2716699 * 2718139 * 2719291 * 2722043 * 2725829 * 2732473 * 2742317 * 2743381 * 2745419 * 2749441 * 2750431 * 2751247 * 2753111 * 2756609 * 2761091 * 2765953 * 2766979 * 2771687 * 2777473 * 2778901 * 2780753 * 2782147 * 2793407 * 2793731 * 2794241 * 2794879 * 2795693 * 2795809 * 2806439 * 2812609 * 2818561 * 2826379 * 2829049 * 2839049 * 2840681 * 2841199 * 2845789 * 2848031 * 2860457 * 2874559 * 2877221 * 2878847 * 2882441 * 2887513 * 2897689 * 2905667 * 2907301 * 2915069 * 2915587 * 2916203 * 2922833 * 2946589 * 2957321 * 2958821 * 2960681 * 2960801 * 2960917 * 2967647 * 2972309 * 2977651 * 2977679 * 2983397 * 2987651 * 2989573 * 2990497 * 2999863 * 3005771 * 3009967 * 3011167 * 3011887 * 3012413 * 3014699 * 3015527 * 3015787 * 3016777 * 3017921 * 3025619 * 3028997 * 3031421 * 3036797 * 3048901 * 3050381 * 3061087 * 3061379 * 3063799 * 3065507 * 3065893 * 3068939 * 3069841 * 3070541 * 3080513 * 3084743 * 3088439 * 3095681 * 3099359 * 3105007 * 3105749 * 3112757 * 3115559 * 3120217 * 3120497 * 3139837 * 3149221 * 3151399 * 3152557 * 3153173 * 3155923 * 3162829 * 3165739 * 3166253 * 3170807 * 3171107 * 3179027 * 3180167 * 3180553 * 3185543 * 3200909 * 3205837 * 3206509 * 3214867 * 3217399 * 3217843 * 3218027 * 3219107 * 3219757 * 3219973 * 3220381 * 3221353 * 3224647 * 3224677 * 3225647 * 3226393 * 3228469 * 3231643 * 3233221 * 3238453 * 3241933 * 3246527 * 3246857 * 3247781 * 3253417 * 3258389 * 3261641 * 3262499 * 3266957 * 3272779 * 3275803 * 3276719 * 3278461 * 3286991 * 3295471 * 3296863 * 3298577 * 3299677 * 3302987 * 3304079 * 3320969 * 3323633 * 3325097 * 3325159 * 3329033 * 3329251 * 3335561 * 3341501 * 3345961 * 3350621 * 3352919 * 3366109 * 3367421 * 3367811 * 3372623 * 3373681 * 3378103 * 3383959 * 3389161 * 3391813 * 3393317 * 3396793 * 3401899 * 3408869 * 3409853 * 3429389 * 3429581 * 3434177 * 3434399 * 3434903 * 3436709 * 3437743 * 3439861 * 3453943 * 3461593 * 3475931 * 3476029 * 3477037 * 3481691 * 3485291 * 3485743 * 3486019 * 3486997 * 3494269 * 3497827 * 3501413 * 3507059 * 3507421 * 3507877 * 3510317 * 3524197 * 3530851 * 3533207 * 3535369 * 3537337 * 3558371 * 3559519 * 3567149 * 3569369 * 3570923 * 3571067 * 3575057 * 3579647 * 3581419 * 3583309 * 3596893 * 3603619 * 3605087 * 3608543 * 3621577 * 3622097 * 3623171 * 3624977 * 3625621 * 3639107 * 3639941 * 3647851 * 3649417 * 3653579 * 3664487 * 3665939 * 3666197 * 3669727 * 3676549 * 3682801 * 3685369 * 3688631 * 3688757 * 3694877 * 3715559 * 3717533 * 3717751 * 3726259 * 3739793 * 3741197 * 3741497 * 3744179 * 3747671 * 3751567 * 3755333 * 3764041 * 3764659 * 3768983 * 3770009 * 3770917 * 3772319 * 3776153 * 3785021 * 3785371 * 3787997 * 3788903 * 3795059 * 3797603 * 3801241 * 3805849 * 3815711 * 3824413 * 3830711 * 3831043 * 3833629 * 3833867 * 3835123 * 3847351 * 3852419 * 3853909 * 3856109 * 3858523 * 3858929 * 3859879 * 3860723 * 3864521 * 3868591 * 3871271 * 3877597 * 3878263 * 3881849 * 3884011 * 3898747 * 3902243 * 3908027 * 3916799 * 3924719 * 3929161 * 3930461 * 3935089 * 3936623 * 3942751 * 3943591 * 3948247 * 3951707 * 3953249 * 3953393 * 3954617 * 3959651 * 3961057 * 3970567 * 3972877 * 3973877 * 3974911 * 3977663 * 3979561 * 3980819 * 3985889 * 3993083 * 3994747 * 3996743 * 4000589 * 4001071 * 4003453 * 4003841 * 4004591 * 4005613 * 4007723 * 4007789 * 4011011 * 4011883 * 4014271 * 4015877 * 4020733 * 4026053 * 4029049 * 4044959 * 4050433 * 4056847 * 4057723 * 4072417 * 4075831 * 4077247 * 4078367 * 4078561 * 4087801 * 4087873 * 4092757 * 4108387 * 4113719 * 4120621 * 4122067 * 4126693 * 4127467 * 4131293 * 4143187 * 4146203 * 4146811 * 4147357 * 4152653 * 4162783 * 4166527 * 4168097 * 4178483 * 4190663
comment:15 in reply to: ↑ 11 Changed 4 years ago by
Replying to bober:
I wonder if the result will be a polynomial which has the correct reduction mod p for lots of p but is wrong somewhere and for which the CRT will miraculously stabilize.
It seems that this is indeed the case here.
comment:16 Changed 4 years ago by
Nice catch. This seems like we hit the unfortunate case of the randomized algorithm. Still I don't understand why:
 only the constant coefficient is erroneous. Normally (i.e. with random matrices), the bitsize of the coefficients grows a bit with the degree, before decreasing, hence, there should be other coefficients for which the theoretical bound was not met when the early termination was decided. Maybe this is due to the special shape of your instance. Could you check whether the bitsize of the coeffs of degree>0 is lower than that of the determinant ?
 the erroneous factor (5) is so small.
comment:17 Changed 4 years ago by
The constant term is the largest, and perhaps it is significantly larger than the others compared to the size of the primes used.
sage: [(n, RR(abs(f[n])).log(2)) for n in range(f.degree() + 1)] [(0, 11507.7121078501), (1, 11445.6397830061), (2, 11382.5167712473), (3, 11320.7650151477), (4, 11255.3983428913), (5, 11194.5632626618), (6, 11126.7122981548), (7, 11067.4853919579), (8, 11001.6807388619), (9, 10939.7232563914), (10, 10874.4419357066), (11, 10811.3885084419), (12, 10746.3717429455), (13, 10682.5576543468), (14, 10617.7114947932), (15, 10553.2870290630), (16, 10488.5711336986), (17, 10423.6198028301), (18, 10359.0179669343), (19, 10293.5899620934), (20, 10229.0995750147), (21, 10163.2246427098), (22, 10098.8524964071), (23, 10032.5453615314), (24, 9968.30626839587), (25, 9901.56837048656), (26, 9837.48555321369), (27, 9770.30407146403), (28, 9706.41136311503), (29, 9638.75492457264), (30, 9575.10182905012), (31, 9506.90990249518), (32, 9443.57272569099), (33, 9374.72836752184), (34, 9311.83785818533), (35, 9242.07747017647), (36, 9179.90936377286), (37, 9108.21101765289), (38, 9047.79795588630), (39, 8976.57271469779), (40, 8915.51312600324), (41, 8845.59047596339), (42, 8783.06331251847), (43, 8713.73701731272), (44, 8650.45604290155), (45, 8581.50166814549), (46, 8517.69805378384), (47, 8449.00548925298), (48, 8384.79539263804), (49, 8316.29994296649), (50, 8251.75350404040), (51, 8183.41307414749), (52, 8118.57730299435), (53, 8050.36262646348), (54, 7985.27123738308), (55, 7917.16099364580), (56, 7851.83934127697), (57, 7783.81744797804), (58, 7718.28528053566), (59, 7650.33927435814), (60, 7584.61239190519), (61, 7516.73240307098), (62, 7450.82371660982), (63, 7383.00178932404), (64, 7316.92202927124), (65, 7249.15165461196), (66, 7182.90986284850), (67, 7115.18564781036), (68, 7048.78953017677), (69, 6981.10695712348), (70, 6914.56314258685), (71, 6846.91839056158), (72, 6780.23262600801), (73, 6712.62243546612), (74, 6645.79973489026), (75, 6578.22130359821), (76, 6511.26606422705), (77, 6443.71696597046), (78, 6376.63305991297), (79, 6309.11118018752), (80, 6241.90202763238), (81, 6174.40551217495), (82, 6107.07414044169), (83, 6039.60135360460), (84, 5972.15044518027), (85, 5904.69993594567), (86, 5837.13186782041), (87, 5769.70234181290), (88, 5702.01921784611), (89, 5634.60951410369), (90, 5566.81319173168), (91, 5499.42226328811), (92, 5431.51437557435), (93, 5364.14127312167), (94, 5296.12324691991), (95, 5228.76710498120), (96, 5160.64017580597), (97, 5093.30020096972), (98, 5025.06542503318), (99, 4957.74088589370), (100, 4889.39914966128), (101, 4822.08936818055), (102, 4753.64139571211), (103, 4686.34573977418), (104, 4617.79209804691), (105, 4550.50997501900), (106, 4481.85107736832), (107, 4414.58192851731), (108, 4345.81803627930), (109, 4278.56133191904), (110, 4209.69255430945), (111, 4142.44778957637), (112, 4073.47408179433), (113, 4006.24077296599), (114, 3937.16193246343), (115, 3869.93961374891), (116, 3800.75527455574), (117, 3733.54349529817), (118, 3664.25312023730), (119, 3597.05144247813), (120, 3527.65431303874), (121, 3460.46230940104), (122, 3390.95751296019), (123, 3323.77476481471), (124, 3254.16117880050), (125, 3186.98727468302), (126, 3117.26354715052), (127, 3050.09808140388), (128, 2980.26260733757), (129, 2913.10517895485), (130, 2843.15607140574), (131, 2776.00628305433), (132, 2705.94133794560), (133, 2638.79879515434), (134, 2568.61544821863), (135, 2501.47975871214), (136, 2431.17503251431), (137, 2364.04580568054), (138, 2293.61624396705), (139, 2226.49309044543), (140, 2155.93467604803), (141, 2088.81720742694), (142, 2018.12525847659), (143, 1951.01308708969), (144, 1880.18212411403), (145, 1813.07486292622), (146, 1742.09843608671), (147, 1674.99569866089), (148, 1603.86615920601), (149, 1536.76755974258), (150, 1465.47575140842), (151, 1398.38090485021), (152, 1326.91573700683), (153, 1259.82425920215), (154, 1188.17209928773), (155, 1121.08360720524), (156, 1049.22738562251), (157, 982.141497611510), (158, 910.059332134092), (159, 842.975668224361), (160, 770.638635072389), (161, 703.556817308040), (162, 630.925083892562), (163, 563.844736694591), (164, 490.860204624151), (165, 423.780955176645), (166, 350.351296274330), (167, 283.272774929579), (168, 209.228831342903), (169, 142.150672040888), (170, 67.0781593020146), (171, 0.000000000000000)]
comment:18 Changed 4 years ago by
What happens with the termination condition if you randomly pick some primes that you have already used for the CRT? I don't find it so easy to read the linbox source code, but judging from the factorization of the difference, you are randomly picking primes between 2^{21} and 2^{22}, or something close to that. That is a relatively small number of possible primes (140336) and in this case the coefficients are really big, so you need a lot of them. And, after you've picked something like 532 primes (so all but one of the coefficients have stabilized) there is about a 1 in 70000 chance that the next two primes you pick will be primes you have already used, so of course nothing will change.
Maybe the 5 is a red herring. The difference has to be divisible by all the primes you used, but perhaps it also has a 1 in p chance of being divisible by any given small prime p.
(I'm just making guesses here, of course, but 1 in 70000 does seem to be about the probability of failure.)
comment:19 Changed 4 years ago by
 Commit changed from 0e8bc88f63a0988fc0047c1bc3c84b588b05290a to 3faea7eeba2d9c83924770b72d8a43871d960f30
Branch pushed to git repo; I updated commit sha1. New commits:
56be9c8  Merge branch 'develop' into t/21579/errors_calculating_characteristic_polynomials_of_rational_matrices

da952e1  replace early terminated Chinese remaindering by the deterministic version for charpoly

3faea7e  remove useless patch application

comment:20 Changed 4 years ago by
Sorry for this long period of inactivity. My latest commit is to
 rebase the ticket to upstream develop
 remove obsolete patch applier in spkginstall
 get rid of the early termination Chinese Remaindering for charpoly. It should fix the last remaining problem, where we got sufficiently unlucky to have picked 3 primes in a row at random which let the algorithm believe that the termination was reached.
As for bober's last question: linbox makes sure to pick distinct primes for the Chinese remaindering. This is not the cause of the bug.
comment:21 Changed 3 years ago by
 Commit changed from 3faea7eeba2d9c83924770b72d8a43871d960f30 to 08012ea976866316645c2c63cf834d5317a50301
Branch pushed to git repo; I updated commit sha1. New commits:
08012ea  Adding a patch fixing a rare bug in charpoly (fflasffpack PR #130)

comment:22 Changed 3 years ago by
After a torture test of fflasffpack charpoly, I found a rare bug, appearing only when a LasVegas probabilistic algorithm failed once but not on the second attempt. I do not know whether it corresponds to the cause of failure here, but it can't hurt. I added the patch to this branch (and it is merged in upstream fflasffpack).
@Bober, do you still manage to reproduce the bug?
comment:24 in reply to: ↑ 23 Changed 3 years ago by
 Milestone changed from sage7.4 to sageduplicate/invalid/wontfix
 Status changed from new to needs_review
comment:25 Changed 3 years ago by
 Status changed from needs_review to positive_review
Ok then, I think it's safe to set it as positive review as it was fixed in #24214.
comment:26 Changed 3 years ago by
 Resolution set to wontfix
 Status changed from positive_review to closed
closing positively reviewed duplicates
Just a comment: PARI supports modular symbols now: http://pari.math.ubordeaux.fr/dochtml/html/cont_Modular_symbols.html (but don't ask me anything more, this is not my mathematics).