libFES : Fast Exhaustive Search for Polynomial Systems over F2

Charles Bouillaguet set up a nice shiny website for libFES the library for exhaustive search on polynomial systems over \mathbb{F}_2. The library has a Sage interface, so it’s easy to get started. It’s also integrated in Charles’ upcoming one-stop boolean system solving patch.

He calls his benchmarketing “bragging rights” … and boy has he earned those rights! Check it out, libFES is fast!

Sage 4.8 is out

If you care about the stuff I care about (and why else would you read this blog?) you might get excited about a few changes in Sage.

Efficient linear algebra for \mathbb{F}_{2^e}

The very first non-trivial patch I ever produced for Sage was about interfacing with NTL for dense linear algebra over \mathbb{F}_{2^e} (I was interested in algebraic attacks against AES at the time). Here’s William’s reply:

Your NTL patch worked perfectly for me first try. I tried more benchmarks (on Pentium-M 1.8Ghz).

[…]

This is pretty good; vastly better than what’s was in SAGE by default, and way better than PARI. Note that MAGMA is much faster though (nearly 8 times faster):

[…]

MAGMA uses (1) […] and (2) a totally different algorithm for computing the echelon form. […] As far as I know, the MAGMA method is not implemented anywhere in the open source world But I’d love to be wrong about that… or even remedy that.

Well, that was 2006. Fast forward to the year 2011 and we get the following timings for computing the reduced row echelon form of a 1,000 x 1,000 matrix over \mathbb{F}_{256}: Sage 4.7.2 takes 36.53 seconds , NTL 5.4.2 takes 31.06 seconds and Magma 2.15 does it in 0.87 seconds. So essentially, the situation didn’t change at all for the better.

With Sage 4.8 this situation changes dramatically  and we get that Sage performs this computation in 0.08 seconds, that’s 450 times faster than Sage 4.7.2. This is because M4RIE was merged in Sage 4.8. Hence, Sage is now (in some cases by far) fastest system to do linear algebra with dense matrices over \mathbb{F}_{2^e} for 1 \leq e \leq 8  and usually also for 9 \leq e \leq 10.

Efficient linear algebra for \mathbb{F}_{p}

One can tell a similar story for \mathbb{F}_p for, say, small to medium sized primes p. In Sage 4.7.2 it took 1.12* seconds to multiply two 1,000 x 1,000 matrices over \mathbb{F}_{251} (although you always had the option to call LinBox explicitly which was way faster but took more memory). With Sage 4.8 the same computation takes *0.16 seconds. For comparison, Magma 2.15 takes 0.22 seconds. So here again Sage moved from poor performance to best in class performance between 4.7.2 and 4.8 simply by making proper use of available libraries.

Viable Alternative yet?

Overall, the story for dense linear algebra in Sage for small finite fields \mathbb{F}_q  is as follows.

q Implementation Comments
2 M4RI Fastest implementation or equal performance depending on platform
3,5,7 \dots LinBox Decent performance, but faster implementations are known in the literature. Also, Magma is a bit faster on my machine.
prime < 2^{23} LinBox Fastest implementation or equal performance depending on platform.
2^e for 2 \leq e \leq 8 M4RIE Fastest
p^e for p>2 Generic Very poor performance, but some work is being done.

So, once we fix that last row Sage finally achieves “viable alternative” quality when it comes to dense linear algebra over \mathbb{F}_{q} if q is q < 2^{16}.

GAP for dense linear algebra over GF(2^e)

For some strange reason it hardly every occurs to me to benchmark GAP when it comes to linear algebra over GF(2^e). Turns out, that’s a really bad idea as it seems to be the best open-source implementation for such things. For example, computing the reduced row echelon form of a 4,000 x 4,000 random matrix over \mathbb{F}_{16} takes 16.94s in GAP 4.4.12 but 3267s in Sage and 1005s in NTL 5.4.2. So much much faster than other established open-source projects out there.

I hope libgap becomes a reality at some point in the future.

PS: For comparison, Magma does the same task in in 4.67s and M4RIE in 0.71s.

Sage, Degrees and Gröbner Bases

For one reason on another I found myself computing Gröbner bases over fields which are not \mathbb{F}_2.  Thus, my interest in Sage’s interface to Singular and Magma increased quite a bit again recently (before I was mainly using PolyBoRi). Hence I wrote two patches which need review (hint, hint):

#10331

A natural question when it comes to Gröbner basis computations is of course how much resources a particular input system is going to require. For many systems (in fact it is conjectured that it holds for most systems) it turns out one can estimate the degree which will be reached during the computation by computing the degree of semi-regularity of the homogenisation of the system, i.e. one computes the first non-zero coefficient of the following power series:

\sum_{k\geq 0} c_k z^k = \frac{\prod_{i=0}^{m-1} (1-z^{d_i})}{(1-z)^n}

where d_i are the respective degrees of the input polynomials. While there is a Magma script readily available for computing the degree of semi-regularity (by Luk Bettale) we didn’t have it for Sage.

sage: sr = mq.SR(1,2,1,4)
sage: F,s = sr.polynomial_system()
sage: I = F.ideal()
sage: I.degree_of_semi_regularity()
3

Thus, we expect to only reach degree three to compute a Gröbner basis  for I:

sage: I = F.ideal()
sage: _ = I.groebner_basis('singular:slimgb',prot='sage')
Leading term degree: 1.
Leading term degree: 1. Critical pairs: 56.
Leading term degree: 2.
Leading term degree: 2. Critical pairs: 69.
Leading term degree: 3.
Leading term degree: 3. Critical pairs: 461.
Leading term degree: 3.
Leading term degree: 3. Critical pairs: 535.
Leading term degree: 2.
Leading term degree: 2. Critical pairs: 423.
Leading term degree: 2.
Leading term degree: 2. Critical pairs: 367.
Leading term degree: 3.
Leading term degree: 3. Critical pairs: 0.
Leading term degree: 1.
Leading term degree: 1. Critical pairs: 39.
Leading term degree: 1. Critical pairs: 38.
Leading term degree: 1. Critical pairs: 37.
Leading term degree: 1. Critical pairs: 36.
Leading term degree: 1. Critical pairs: 35.
Leading term degree: 1. Critical pairs: 34.
Leading term degree: 1. Critical pairs: 33.
Leading term degree: 1. Critical pairs: 32.
Leading term degree: 1. Critical pairs: 31.
Leading term degree: 1. Critical pairs: 30.
Leading term degree: 1. Critical pairs: 29.
Leading term degree: 1. Critical pairs: 28.
Leading term degree: 1. Critical pairs: 27.
Leading term degree: 1. Critical pairs: 26.
Leading term degree: 1. Critical pairs: 25.
Leading term degree: 1. Critical pairs: 24.
Leading term degree: 1. Critical pairs: 23.
Leading term degree: 1. Critical pairs: 22.
Leading term degree: 1. Critical pairs: 21.
Leading term degree: 1. Critical pairs: 20.
Leading term degree: 1. Critical pairs: 19.
Leading term degree: 1. Critical pairs: 18.
Leading term degree: 1. Critical pairs: 17.
Leading term degree: 1. Critical pairs: 16.
Leading term degree: 1. Critical pairs: 15.
Leading term degree: 1. Critical pairs: 14.
Leading term degree: 1. Critical pairs: 13.
Leading term degree: 1. Critical pairs: 12.
Leading term degree: 1. Critical pairs: 11.
Leading term degree: 1. Critical pairs: 10.
Leading term degree: 1. Critical pairs: 9.
Leading term degree: 1. Critical pairs: 8.
Leading term degree: 1. Critical pairs: 7.
Leading term degree: 1. Critical pairs: 6.
Leading term degree: 1. Critical pairs: 5.
Leading term degree: 1. Critical pairs: 4.
Leading term degree: 1. Critical pairs: 3.
Leading term degree: 1. Critical pairs: 2.

Highest degree reached during computation: 3.

Which leads me to my next point.

#10571

It’s nice to see how far a particular computation is progressed and to get some summary information about the  computation in the end. Hence, a new patch which enables live logs from Singular and Magma in Sage.

sage: _ = I.groebner_basis('singular:slimgb',prot=True) # Singular native
1461888602
> def sage584=slimgb(sage581);
def sage584=slimgb(sage581);
1M[23,23](56)2M[56,eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeebbbbb32](69)3M[69,bbb69](461)3M[100,bbbbbb74](535)2M[100,beebbbbbb3](423)2M[23,eeeebbbbbb0](367)3M[28,eeeeebb0](0)
NF:399 product criterion:16403, ext_product criterion:0
[7:3]1(39)s(38)s(37)s(36)s(35)s(34)s(33)s(32)s(31)s(30)s(29)s(28)s(27)s(26)s(25)s(24)s(23)s(22)s(21)s(20)s(19)s(18)s(17)s(16)s(15)s(14)s(13)s(12)s(11)s(10)s(9)s(8)s(7)s(6)s(5)s(4)s(3)s(2)sss
(S:39)---------------------------------------
>

sage: _ = I.groebner_basis('magma', prot=True) # Magma native
Append(~_sage_, 0);
Append(~_sage_, 0);
>>>_sage_[126]:=_sage_[128];
_sage_[126]:=_sage_[128];
>>>Append(~_sage_, 0);
Append(~_sage_, 0);
>>>_sage_[86]:=GroebnerBasis(_sage_[126]);
_sage_[86]:=GroebnerBasis(_sage_[126]);
Homogeneous weights search
Number of variables: 40, nullity: 0
Exact search time: 0.000
********************
FAUGERE F4 ALGORITHM
********************
Coefficient ring: GF(2^4)
Rank: 40
Order: Graded Reverse Lexicographical
NEW hash table
Matrix kind: Zech short
Datum size: 2
No queue sort
Initial length: 80
Inhomogeneous

Initial queue setup time: 0.009
Initial queue length: 48

*******
STEP 1
Basis length: 80, queue length: 48, step degree: 1, num pairs: 8
Basis total mons: 264, average length: 3.300
Number of pair polynomials: 8, at 25 column(s), 0.000
Average length for reductees: 6.00 [8], reductors: 10.00 [8]
Symbolic reduction time: 0.000, column sort time: 0.000
8 + 8 = 16 rows / 25 columns, 32% / 37.641% (8/r)
Before ech memory: 7.9MB
Row sort time: 0.000
0.000 + 0.000 = 0.000 [8]
Delete 1 memory chunk(s); time: 0.000
After ech memory: 7.9MB
Queue insertion time: 0.000
Step 1 time: 0.000, [0.001], mat/total: 0.000/0.009 [0.005], mem: 7.9MB

*******
STEP 2
Basis length: 88, queue length: 48, step degree: 2, num pairs: 32
Basis total mons: 340, average length: 3.864
Number of pair polynomials: 32, at 169 column(s), 0.000
Average length for reductees: 3.88 [32], reductors: 7.12 [192]
Symbolic reduction time: 0.000, column sort time: 0.000
32 + 192 = 224 rows / 293 columns, 2.2733% / 5.8884% (6.6607/r)
Before ech memory: 7.9MB
Row sort time: 0.000
0.000 + 0.000 = 0.000 [8]
Delete 1 memory chunk(s); time: 0.000
After ech memory: 7.9MB
Queue insertion time: 0.000
Step 2 time: 0.000, [0.003], mat/total: 0.000/0.009 [0.008], mem: 7.9MB

*******
STEP 3
Basis length: 96, queue length: 69, step degree: 3, num pairs: 69
Basis total mons: 472, average length: 4.917
Number of pair polynomials: 69, at 540 column(s), 0.000
Average length for reductees: 13.20 [69], reductors: 5.50 [343]
Symbolic reduction time: 0.000, column sort time: 0.000
69 + 343 = 412 rows / 719 columns, 0.94387% / 2.4223% (6.7864/r), bv
Before ech memory: 7.9MB
Row sort time: 0.000
0.000 + 0.000 = 0.000 [69]
Delete 1 memory chunk(s); time: 0.000
After ech memory: 7.9MB
Queue insertion time: 0.010
Step 3 time: 0.010, [0.015], mat/total: 0.000/0.019 [0.023], mem: 7.9MB

*******
STEP 4
Basis length: 165, queue length: 736, step degree: 3, num pairs: 461
Basis total mons: 1368, average length: 8.291
Number of pair polynomials: 461, at 802 column(s), 0.000
Average length for reductees: 10.35 [461], reductors: 7.49 [654]
Symbolic reduction time: 0.000, column sort time: 0.000
461 + 654 = 1115 rows / 804 columns, 1.0789% / 2.7301% (8.6744/r)
Before ech memory: 7.9MB
Row sort time: 0.000
0.010 + 0.000 = 0.010 [129]
Delete 1 memory chunk(s); time: 0.000
Number of unused reductors: 2
After ech memory: 7.9MB
Queue insertion time: 0.020
Step 4 time: 0.030, [0.022], mat/total: 0.010/0.049 [0.045], mem: 7.9MB

*******
STEP 5
Basis length: 294, queue length: 2094, step degree: 2, num pairs: 132
Basis total mons: 1669, average length: 5.677
Number of pair polynomials: 132, at 149 column(s), 0.000
Average length for reductees: 2.00 [132], reductors: 5.22 [148]
Symbolic reduction time: 0.000, column sort time: 0.000
132 + 148 = 280 rows / 149 columns, 2.4856% / 7.1242% (3.7036/r)
Before ech memory: 7.9MB
Row sort time: 0.000
0.000 + 0.000 = 0.000 [0]
Delete 1 memory chunk(s); time: 0.000
After ech memory: 7.9MB
Queue insertion time: 0.000
Step 5 time: 0.000, [0.001], mat/total: 0.010/0.049 [0.046], mem: 7.9MB

*******
STEP 6
Basis length: 294, queue length: 1962, step degree: 3, num pairs: 848
Basis total mons: 1669, average length: 5.677
Number of pair polynomials: 57, at 85 column(s), 0.000
Average length for reductees: 2.00 [57], reductors: 2.51 [84]
Symbolic reduction time: 0.000, column sort time: 0.000
57 + 84 = 141 rows / 85 columns, 2.7117% / 8.4127% (2.305/r)
Before ech memory: 7.9MB
Row sort time: 0.000
0.000 + 0.000 = 0.000 [0]
Delete 1 memory chunk(s); time: 0.000
After ech memory: 7.9MB
Queue insertion time: 0.000
Step 6 time: 0.000, [0.000], mat/total: 0.010/0.049 [0.046], mem: 7.9MB

*******
STEP 7
Basis length: 294, queue length: 1114, step degree: 4, num pairs: 1098
Basis total mons: 1669, average length: 5.677
Number of pair polynomials: 0, at 0 column(s), 0.000
Symbolic reduction time: 0.000, column sort time: 0.000
0 + 0 = 0 rows / 0 columns, 0% / 0% (0/r), bv
Before ech memory: 7.9MB
Row sort time: 0.000
0.000 + 0.000 = 0.000 [0]
After ech memory: 7.9MB
Queue insertion time: 0.000
Step 7 time: 0.000, [0.000], mat/total: 0.010/0.049 [0.046], mem: 7.9MB

*******
STEP 8
Basis length: 294, queue length: 16, step degree: 5, num pairs: 16
Basis total mons: 1669, average length: 5.677
Number of pair polynomials: 0, at 0 column(s), 0.000
Symbolic reduction time: 0.000, column sort time: 0.000
0 + 0 = 0 rows / 0 columns, 0% / 0% (0/r), bv
Before ech memory: 7.9MB
Row sort time: 0.000
0.000 + 0.000 = 0.000 [0]
After ech memory: 7.9MB
Queue insertion time: 0.000
Step 8 time: 0.000, [0.000], mat/total: 0.010/0.049 [0.046], mem: 7.9MB

Reduce 294 final polynomial(s) by 294
16 redundant polynomial(s) removed; time: 0.000
Interreduce 40 (out of 294) polynomial(s)
Symbolic reduction time: 0.000
Column sort time: 0.000
40 + 0 = 40 rows / 41 columns, 10.976% / 24.736% (4.5/r)
Row sort time: 0.000
0.000 + 0.000 = 0.000 [40]
Delete 1 memory chunk(s); time: 0.000
Total reduction time: 0.000
Reduction time: 0.000
Final number of polynomials: 278

Number of pairs: 759
Total pair setup time: 0.000
Max num entries matrix: 1115 by 804
Max num rows matrix: 1115 by 804
Total symbolic reduction time: 0.000
Total column sort time: 0.000
Total row sort time: 0.000
Total matrix time: 0.010
Total new polys time: 0.000
Total queue update time: 0.030
Total Faugere F4 time: 0.049, real time: 0.046
>>>_sage_[126]:=0;
_sage_[126]:=0;
>>>

sage: _ = I.groebner_basis('singular:slimgb',prot='sage') # Singular parsed
Leading term degree: 1.
Leading term degree: 1. Critical pairs: 56.
Leading term degree: 2.
Leading term degree: 2. Critical pairs: 69.
Leading term degree: 3.
Leading term degree: 3. Critical pairs: 461.
Leading term degree: 3.
Leading term degree: 3. Critical pairs: 535.
Leading term degree: 2.
Leading term degree: 2. Critical pairs: 423.
Leading term degree: 2.
Leading term degree: 2. Critical pairs: 367.
Leading term degree: 3.
Leading term degree: 3. Critical pairs: 0.
Leading term degree: 1.
Leading term degree: 1. Critical pairs: 39.
Leading term degree: 1. Critical pairs: 38.
Leading term degree: 1. Critical pairs: 37.
Leading term degree: 1. Critical pairs: 36.
Leading term degree: 1. Critical pairs: 35.
Leading term degree: 1. Critical pairs: 34.
Leading term degree: 1. Critical pairs: 33.
Leading term degree: 1. Critical pairs: 32.
Leading term degree: 1. Critical pairs: 31.
Leading term degree: 1. Critical pairs: 30.
Leading term degree: 1. Critical pairs: 29.
Leading term degree: 1. Critical pairs: 28.
Leading term degree: 1. Critical pairs: 27.
Leading term degree: 1. Critical pairs: 26.
Leading term degree: 1. Critical pairs: 25.
Leading term degree: 1. Critical pairs: 24.
Leading term degree: 1. Critical pairs: 23.
Leading term degree: 1. Critical pairs: 22.
Leading term degree: 1. Critical pairs: 21.
Leading term degree: 1. Critical pairs: 20.
Leading term degree: 1. Critical pairs: 19.
Leading term degree: 1. Critical pairs: 18.
Leading term degree: 1. Critical pairs: 17.
Leading term degree: 1. Critical pairs: 16.
Leading term degree: 1. Critical pairs: 15.
Leading term degree: 1. Critical pairs: 14.
Leading term degree: 1. Critical pairs: 13.
Leading term degree: 1. Critical pairs: 12.
Leading term degree: 1. Critical pairs: 11.
Leading term degree: 1. Critical pairs: 10.
Leading term degree: 1. Critical pairs: 9.
Leading term degree: 1. Critical pairs: 8.
Leading term degree: 1. Critical pairs: 7.
Leading term degree: 1. Critical pairs: 6.
Leading term degree: 1. Critical pairs: 5.
Leading term degree: 1. Critical pairs: 4.
Leading term degree: 1. Critical pairs: 3.
Leading term degree: 1. Critical pairs: 2.

Highest degree reached during computation: 3.

Note, that these logs are all live, i.e. the are displayed during the computation as soon as the respective system makes them available.

Both patches are up for review on Sage’s trac server.