Test failure with `multiple_of_order_using_frobp`

[user202729]

Steps To Reproduce

src/bin/sage -t --warn-long 5.0 --random-seed=210700063500461381029255797937294367397 src/sage/modular/abvar/torsion_subgroup.py
**********************************************************************
Error: Failed example:: Got: 1
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
    G.multiple_of_order_using_frobp()
Expected:
    97
Got:
    1
**********************************************************************
Error: Failed example:: Got: [92645296242160800, 7275, 291, 1]
                                                                                                                                                                                                                                                                                                                                                                                                                                                          
    [G.multiple_of_order_using_frobp(p) for p in prime_range(3,13)]
Expected:
    [92645296242160800, 7275, 291, 97]
Got:
    [92645296242160800, 7275, 291, 1]
**********************************************************************
Error: Failed example:: Got: [92645296242160800, 7275, 291, 1, 1, 1]
                                                                                                                                                                                                                                                                                                                                                                                                                                                    
    [G.multiple_of_order_using_frobp(p) for p in prime_range(3,19)]
Expected:
    [92645296242160800, 7275, 291, 97, 97, 97]
Got:
    [92645296242160800, 7275, 291, 1, 1, 1]
**********************************************************************
1 item had failures:
   3 of  25 in sage.modular.abvar.torsion_subgroup.RationalTorsionSubgroup.multiple_of_order_using_frobp

Expected Behavior

Actual Behavior

Additional Information

No response

Environment

• OS:
• Sage Version:

Checklist

• I have searched the existing issues for a bug report that matches the one I want to file, without success.
• I have read the documentation and troubleshoot guide

[user202729]

might be linbox issue again? cf #40856

[user202729]

or maybe not, although for now I can't rule out that is linbox is not used yet.

For reference, the test looks like this

J = J0(389)
G = J.rational_torsion_subgroup(); G
G.multiple_of_order_using_frobp()

[user202729]

note that you can't increase the probability of reproducing it by repeatedly running this a lot of times because of the caching. As you can observe the first run is slow (200ms) but subsequent times not so much.

here is %prun of first time run

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        6    0.078    0.013    0.084    0.014 ambient.py:2908(_compute_hecke_matrix_prime)
        6    0.070    0.012    0.164    0.027 hecke_operator.py:270(charpoly)
        8    0.028    0.004    0.028    0.004 {built-in method builtins.compile}
      279    0.018    0.000    0.019    0.000 dynamic_class.py:339(dynamic_class_internal)
        1    0.017    0.017    0.038    0.038 relation_matrix.py:258(gens_to_basis_matrix)
       30    0.008    0.000    0.010    0.000 matrix_space.py:1051(_element_constructor_)
       17    0.005    0.000    0.005    0.000 {built-in method _imp.create_dynamic}
  231/163    0.004    0.000    0.025    0.000 category.py:2332(join)
        1    0.004    0.004    0.020    0.020 matrix_morphism.py:864(kernel)
     4615    0.004    0.000    0.004    0.000 manin_symbol_list.py:317(manin_symbol)
     1274    0.003    0.000    0.006    0.000 verbose.py:116(verbose)
        6    0.003    0.000    0.089    0.015 submodule.py:224(_compute_hecke_matrix)
    90/70    0.003    0.000    0.010    0.000 category.py:839(_all_super_categories)
     17/7    0.002    0.000    0.006    0.001 {built-in method _imp.exec_dynamic}
       65    0.002    0.000    0.002    0.000 element.py:307(<listcomp>)

indeed characteristic polynomial is called.

    def charpoly(self, var='x'):
        """
        Return the characteristic polynomial of this Hecke operator.

        INPUT:

        - ``var`` -- string (default: ``'x'``)

        OUTPUT: a monic polynomial in the given variable

        EXAMPLES::

            sage: M = ModularSymbols(Gamma1(6),4)
            sage: M.hecke_operator(2).charpoly('x')
            x^6 - 14*x^5 + 29*x^4 + 172*x^3 - 124*x^2 - 320*x + 256
        """
        return self.matrix().charpoly(var)

where the self could be e.g.

Hecke operator T_3 on Modular Symbols subspace of dimension 64 of Modular Symbols space of dimension
 65 for Gamma_0(389) of weight 2 with sign 0 over Rational Field

and the matrix is e.g.

64 x 64 dense matrix over Rational Field

charpoly of rational_matrix has

        poly = self.fetch('charpoly')
        if poly is not None:
            return poly.change_variable_name(var)

        if algorithm is None:
            algorithm = 'flint' if self._nrows <= 40 else 'linbox'

        if algorithm == 'flint' or algorithm == 'linbox':
            A, denom = self._clear_denom()
            f = A.charpoly(var, algorithm=algorithm)
            x = f.parent().gen()
            g = f(x * denom) / denom ** f.degree()
        elif algorithm == 'generic':
            g = Matrix_dense.charpoly(self, var)
        else:
            raise ValueError("no algorithm '%s'" % algorithm)

        self.cache('charpoly', g)
        return g

so indeed linbox would use this. This in turn calls the charpoly method for integers…

to get the matrix, one need to do either

G.abelian_variety().hecke_operator(3).matrix()
G.abelian_variety().modular_symbols().hecke_operator(3).matrix()

I've no idea why one returns a matrix over ℤ but the other returns a matrix over ℚ, but their values are equal. the version passing through modular_symbols is much faster (144ms versus 2.75s)

[user202729]

apparently on my machine, there are only 10^9 distinct values of seed (the nanosecond part of time). not that it's feasible to brute force through all of them but…?

On the other hand each charpoly() needs two seeds.

[user202729]

there's a benchmark in #22970 (comment) . If flint has gotten faster than linbox then we can just default to flint

[user202729]

okay, I find a failing seed value.

to reproduce the issue

J = J0(389)
G = J.rational_torsion_subgroup(); G
m = G.abelian_variety().modular_symbols().hecke_operator(11).matrix()
m._clear_cache()
m.charpoly()

this returns 1 when seed is 818466634. (to modify seed, I patch /include/linbox/randiter/random-prime.h like this:

		PrimeIterator(uint64_t bits = 23, uint64_t seed = 0) :
                _bits(bits)
            {
                linbox_check(bits >1);
                if (! seed)
                {
                    //seed = BaseTimer::seed();
                    seed = 818466634;
                    __builtin_printf("seed=%llu\n", (unsigned long long) seed);
                }

[user202729]

If anyone see something similar in the future: to stress test quickly, do something like

p = m.charpoly()
while True:
    m._clear_cache()
    q = m.charpoly()
    assert q == p

together with patching linbox to print out the seed.

[user202729]

I tried running @videlec 's script above and get

sage: test_charpoly()
dim=3 density=0.7778 nbits(height)=2
flint        0.0479
linbox       0.5193
generic      0.0553

dim=4 density=0.5000 nbits(height)=2
flint        0.0830
linbox       0.5498
generic      0.0508

dim=5 density=0.5600 nbits(height)=2
flint        0.0850
linbox       0.5522
generic      0.0530

dim=6 density=0.6389 nbits(height)=2
flint        0.0886
linbox       0.5644
generic      0.0954

dim=10 density=0.3900 nbits(height)=2
flint        0.1090
linbox       0.6863
generic      0.2985

dim=10 density=1.0000 nbits(height)=2
flint        0.1100
linbox       0.6361
generic      0.7016

dim=30 density=0.0900 nbits(height)=2
flint        0.5770
linbox       4.2802
generic      8.9204

dim=30 density=0.3833 nbits(height)=2
flint        0.6313
linbox       1.6566
generic      36.1790

dim=30 density=1.0000 nbits(height)=2
flint        0.6399
linbox       1.7095
generic      71.5537

flint wins in almost all cases, only lose to generic in a few cases, although these looks negligible.