Opened 3 years ago

Closed 15 months ago

# Errors calculating characteristic polynomials of rational matrices

Reported by: Owned by: bober major sage-duplicate/invalid/wontfix linear algebra cpernet N/A u/cpernet/errors_calculating_characteristic_polynomials_of_rational_matrices (Commits) 08012ea976866316645c2c63cf834d5317a50301 #24214

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:~/sage-bug\$ 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:~/sage-bug\$ 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:~/sage-bug\$ 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:~/sage-bug\$ sage bug.sage 3
linbox (maybe) correct? True
2nd linbox (maybe) correct? True

jb12407@lmfdb5:~/sage-bug\$ 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.

### comment:1 Changed 3 years ago by jdemeyer

Just a comment: PARI supports modular symbols now: http://pari.math.u-bordeaux.fr/dochtml/html/cont_Modular_symbols.html (but don't ask me anything more, this is not my mathematics).

### comment:2 Changed 3 years ago by bober

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

It would be good if you could reduce the issue further :-)

### comment:4 Changed 3 years ago by bober

sigh...

I think that while trying to diagnose one bug I found another, which is now #21596.

### comment:6 Changed 3 years ago by cpernet

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 3 years ago by bober

(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, pre-Linbox 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 3 years ago by cpernet

I think I finally found the bug(s). The one causing the never-ending computation is in Givaro's random iterator choice for a seed. It used the 6 digits of the micro-seconds 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 non-zero 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 fflas-ffpack. Could you let me know if it fixes all your bugs/infinite loop problem.

### Changed 3 years ago by cpernet

Fix to Givaro randiter

### Changed 3 years ago by cpernet

Fix to fflas-ffpack LUKrylov charpoly

### comment:9 Changed 3 years ago by cpernet

• Branch set to u/cpernet/errors_calculating_characteristic_polynomials_of_rational_matrices

### comment:10 Changed 3 years ago by cpernet

• Commit set to 0e8bc88f63a0988fc0047c1bc3c84b588b05290a

The 2 patches are now in the ticket's branch.

New commits:

 ​0e8bc88 `add givaro and fflas-ffpack patches attached to ticket #21579`

### comment:11 follow-up: ↓ 15 Changed 3 years ago by bober

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 package-version.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 3 years ago by bober

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 3 years ago by bober

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 3 years ago by bober

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

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 3 years ago by cpernet

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.
Last edited 3 years ago by cpernet (previous) (diff)

### comment:17 Changed 3 years ago by bober

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 3 years ago by bober

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 221 and 222, 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 2 years ago by git

• 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 2 years ago by cpernet

Sorry for this long period of inactivity. My latest commit is to

• rebase the ticket to upstream develop
• remove obsolete patch applier in spkg-install
• 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 2 years ago by git

• 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 (fflas-ffpack PR #130)`

### comment:22 Changed 2 years ago by cpernet

After a torture test of fflas-ffpack charpoly, I found a rare bug, appearing only when a Las-Vegas 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 fflas-ffpack).

@Bober, do you still manage to reproduce the bug?

### comment:23 follow-up: ↓ 24 Changed 21 months ago by cpernet

• Dependencies set to #24214

#24214 should fix it.

### comment:24 in reply to: ↑ 23 Changed 17 months ago by klee

• Milestone changed from sage-7.4 to sage-duplicate/invalid/wontfix
• Status changed from new to needs_review

#24214 should fix it.

Sure. bober's all tests now pass. So we close this?

### comment:25 Changed 17 months ago by cpernet

• 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 15 months ago by vdelecroix

• Resolution set to wontfix
• Status changed from positive_review to closed

closing positively reviewed duplicates

Note: See TracTickets for help on using tickets.