CoCalc -- 18.783 Lecture 3 Root finding in finite fields.ipynb

⁶⁷ views
ubuntu2404

Kernel: SageMath 10.7

Problem:	Given a polynomial $f\in \mathbb F_q[x]$ , determine the $\mathbb F_q$ -rational roots of $f$ .
Solution:	Rabin's algorithm!

In [3]:

p=61; Fp=GF(p); R.<x>=PolynomialRing(Fp); f = x^8-2*x+5

In [4]:

f.roots()

Out[4]:

[(44, 1), (39, 1), (15, 1)]

In [5]:

f.factor()

Out[5]:

(x + 17) * (x + 22) * (x + 46) * (x^2 + 46*x + 1) * (x^3 + 52*x^2 + 41*x + 33)

In [7]:

f.change_ring(GF(p^6)).roots()

Out[7]:

[(44, 1),
 (39, 1),
 (15, 1),
 (46*z6^5 + 32*z6^4 + 56*z6^3 + 22*z6^2 + 48*z6 + 44, 1),
 (7*z6^5 + 35*z6^4 + 2*z6^3 + 4*z6 + 37, 1),
 (15*z6^5 + 29*z6^4 + 5*z6^3 + 39*z6^2 + 13*z6 + 32, 1),
 (28*z6^5 + 30*z6^4 + 20*z6^3 + 26*z6^2 + 21*z6 + 30, 1),
 (26*z6^5 + 57*z6^4 + 39*z6^3 + 35*z6^2 + 36*z6 + 3, 1)]

Theorem:	$\mathbb F_{p^n} = \{\alpha\in \overline{\mathbb F}_p: \alpha^{p^n}=\alpha\}$ , in other words, we have $x^{p^n}-x=\prod_{\alpha\in \mathbb F_{p^n}}(x-\alpha)$
Proof:	For $n=1$ this follows from Fermat's little theorem. If we identify $\mathbb F_p\simeq \mathbb Z/p\mathbb Z$ , we have $a^{p-1}=1\bmod p$ for $a\in [0,p-1]$ , so the $p$ elements of $\mathbb Z/p\mathbb Z$ are precisely the roots of $x^p-x$ . For $n>1$ this is true by definition: $\mathbb F_{p^n}:=\overline{\mathbb{F}}_{p}^{\langle \alpha \mapsto \alpha^{p^n}\rangle}$ is the fixed field of the $p^n$ -power Frobenius automorphism of $\overline{\mathbb F}_p$ , equivalently, the splitting field of $x^{p^n}-x$ over $\mathbb F_p$ .
Corollary:	For any prime power $q$ we have $\mathbb F_{q^n} = \{\alpha\in \overline{\mathbb F}_q: \alpha^{q^n}=\alpha\}$ and $x^{q^n}-x=\prod_{\alpha\in \mathbb F_{q^n}}(x-\alpha)$ .

In [8]:

product([x-i for i in Fp]) == x^p-x

Out[8]:

True

In [9]:

product([x.change_ring(GF(p^2))-i for i in GF(p^2)]) == x^(p^2)-x

Out[9]:

True

Corollary:

Let

f\in \mathbb F_q[x]

be nonzero and let

S:=\{\alpha\in \mathbb F_{q^n}:f(\alpha)=0\}

. Then

\gcd(f(x),x^{q^n}-x)=\prod_{\alpha\in S}(x-\alpha).

In particular,

\#S=\deg(\gcd(f(x),x^{q^n}-x)

As usual,

\gcd(a,b)

denotes the unique monic representative of the nonzero principal ideal

(a,b)

\mathbb F_q[x]

In [10]:

%time gcd(f,x^p-x).degree()

Out[10]:

CPU times: user 53 µs, sys: 10 µs, total: 63 µs
Wall time: 66.8 µs

3

In [11]:

%time gcd(f,x^(p^2)-x).degree()

Out[11]:

CPU times: user 558 µs, sys: 0 ns, total: 558 µs
Wall time: 2.31 ms

5

In [13]:

%time gcd(f,x^(p^3)-x).degree()

Out[13]:

CPU times: user 29.4 ms, sys: 146 µs, total: 29.6 ms
Wall time: 29.3 ms

6

In [14]:

%time gcd(f,x^(p^6)-x).degree()

Out[14]:

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
File <timed eval>:1
File /ext/sage/10.7/src/sage/rings/polynomial/polynomial_zmod_flint.pyx:534, in sage.rings.polynomial.polynomial_zmod_flint.Polynomial_zmod_flint.__pow__()
    532             return (<Polynomial_zmod_flint>self)._powmod_bigexp(exp, modulus)
    533 
--> 534     return Polynomial_template.__pow__(self, exp, modulus)
    535 
    536 cdef Polynomial _powmod_bigexp(Polynomial_zmod_flint self, Integer exp, Polynomial_zmod_flint m):
File /ext/sage/10.7/src/sage/rings/polynomial/polynomial_template.pxi:658, in sage.rings.polynomial.polynomial_zmod_flint.Polynomial_template.__pow__()
    656 
    657         if modulus is None:
--> 658             celement_pow(&r.x, &(<Polynomial_template>self).x, e, NULL, (<Polynomial_template>self)._cparent)
    659         else:
    660             if parent is not (<Polynomial_template>modulus)._parent and parent != (<Polynomial_template>modulus)._parent:
File /ext/sage/10.7/src/sage/libs/flint/nmod_poly_linkage.pxi:548, in sage.rings.polynomial.polynomial_zmod_flint.celement_pow()
    546     sig_off()
    547 else:
--> 548     sig_on()
    549     nmod_poly_pow(res, x, e)
    550     sig_off()
KeyboardInterrupt: 

Key point:

Make sure you exponentiate in the right ring!

The first step in computing

\gcd(f(x),x^q-x)

is to compute

x^q-x\bmod f

(note that

q

may be exponentially large, but

\deg f

certainly won't be).
If

\pi_f\colon \mathbb{F}_q[x]\to\mathbb{F}_q[x]/(f)

is the ring homomorphism defined by

x\mapsto x\bmod f

, then

\pi_f(x^q-x)=\pi_f(x^q)-\pi_f(x)=\pi_f(x)^q-\pi_f(x),

so the order of operations makes no mathematical difference, but the RHS can be computed exponentially faster than either the middle or the LHS!

The complexity of computing

x^q\bmod f

with Euclidean division is quasi-linear in

q

,
The complexity of computing

\pi_f(x)^q

using binary exponentiation is quasi-linear in $\log q$ .

In [15]:

pif.<z>=R.quotient(f)

In [16]:

%time print(pif(x)^(p^5)-pif(x))

Out[16]:

4*z^7 + 3*z^6 + 47*z^5 + 11*z^4 + 14*z^3 + 7*z^2 + 54*z + 54
CPU times: user 0 ns, sys: 3.22 ms, total: 3.22 ms
Wall time: 3.84 ms

In [17]:

%time print(pif(x)^(p^10)-pif(x))

Out[17]:

17*z^7 + 8*z^6 + 54*z^5 + 24*z^4 + 56*z^3 + 56*z^2 + 39*z + 16
CPU times: user 875 µs, sys: 0 ns, total: 875 µs
Wall time: 883 µs

In [18]:

p=2^255-19;  F=GF(p); R.<x>=PolynomialRing(F); f=x^8-2*x+5; pif.<z>=R.quotient(f);
%time print(f.factor())

Out[18]:

(x + 51027038539503343326764519138825597294378744664834164009924300106595537523144) * (x^2 + 50611585019097526067500995121554429372177258784480155764381858723769139594768*x + 38688838658452246841078006525201819047008644554995237691159504617557701189662) * (x^5 + 14153465678715326029305470748307881186713981216326244265151425177548452521986*x^4 + 30970415661099139061786692430821851662201464041579196405043008572374691402945*x^3 + 499132331909190267118163511147975581363888184276674106555956712514448901935*x^2 + 8376079156473006437715373698581520278472986811325551025499711711353793595153*x + 26683694078622147531159282359840405681456757722151517437587006163097167872535)
CPU times: user 81.1 ms, sys: 115 µs, total: 81.2 ms
Wall time: 113 ms

In [19]:

%time print(pif(x)^p-pif(x))
%time print(pif(x)^(p^2)-pif(x))
%time print(pif(x)^(p^5)-pif(x))
%time print(pif(x)^(p^10)-pif(x))

Out[19]:

56912153353292371673289258375060596419750666264775183997983451927556626280597*z^7 + 54126375986844650798941255777254763470210137102295110078507311106073176629557*z^6 + 26370191539067697445434004938354527613678587940042130721275548622177241819874*z^5 + 33849979231143270125482603335989690944151184574514800928473029898163357518769*z^4 + 14890534162725728972685758677938903637622565981985966297127122082293952614007*z^3 + 32092641922173492064244277474385096364212011307267240087973213220570821403539*z^2 + 19373733151276492061680343909784983942865853287652217931214007763332010881723*z + 1721809714390020567368409726245875637047570619078921538054345133699892443866
CPU times: user 29.2 ms, sys: 0 ns, total: 29.2 ms
Wall time: 28.1 ms
21481264858722799008137586088768265748406991549253112395949112392374122246421*z^7 + 30454705691416756665559836626909501648854654660566544683647698737233552665571*z^6 + 30071754228606174566974412478038911889313403236686985894296111795798760887565*z^5 + 6626651469030646591403975846388523563961268257092576649285317646742810445354*z^4 + 13734065835341445050380685382411230286132530737420188776974535208986900394204*z^3 + 3733579156307926464850917376908926679146433210878375104939753916552117231199*z^2 + 21130108819928400076858645456221967025028186384174163686091706670275540557871*z + 5829819961228675019598343722913500385264009038422264821885647316312709683725
CPU times: user 25.6 ms, sys: 0 ns, total: 25.6 ms
Wall time: 25.4 ms
39852906293622564262290251482505193106960838840295569642385992100007716813453*z^7 + 49031094070818215186020428043541136580072121910298010885525896127940080431690*z^6 + 8824657334656670085288568905347554613799480630449365182254624336871003384652*z^5 + 2525490141248154954439090036807371260760412242530540629724967806510296152832*z^4 + 2688011256477899070540064756639638331871251655577944155526446365433297330662*z^3 + 54245435696182649328554527698021331561801752842463385322324123351693497741625*z^2 + 5135820341176756270775854908699968559990922783283432765689944150468863154209*z + 31575492068812183536616671788217831434429768611918246660349899831910484011368
CPU times: user 54.6 ms, sys: 0 ns, total: 54.6 ms
Wall time: 54.6 ms
0
CPU times: user 92 ms, sys: 0 ns, total: 92 ms
Wall time: 92.1 ms

In [20]:

def distinct_root_poly(f):
    q = f.base_ring().cardinality()
    pif = f.parent().quotient(f)
    return gcd(f,(pif(x)^q).lift()-x)

def count_distinct_roots(f):
    return distinct_root_poly(f).degree()

In [23]:

p=61;  F=GF(p); R.<x>=PolynomialRing(F); f=x^8-2*x+5; pif=R.quotient(f)
for n in range(1,7):
    print("%d distinct roots over F_{%d^%d}" % (count_distinct_roots(f.change_ring(GF(p^n))),p,n))

Out[23]:

distinct roots over F_{61^1}
distinct roots over F_{61^2}
distinct roots over F_{61^3}
distinct roots over F_{61^4}
distinct roots over F_{61^5}
distinct roots over F_{61^6}

Okay, we can count (distinct) roots, but what if we actually want to know what they are?

In [24]:

gcd(f,(pif(x)^p-x).lift())

Out[24]:

x^3 + 24*x^2 + 33*x + 2

In [25]:

s = (p-1) // 2 # we assume p (or more generally q) is odd
gcd(f,(pif(x)^s-1).lift())

Out[25]:

x^2 + 7*x + 36

In [26]:

print(gcd(f,(pif(x)^p-x).lift())/gcd(f,(pif(x)^s-1).lift()))

Out[26]:

x + 17

Excellent, we managed to split the roots of $f$ by picking out two that are roots of $x^s-1$ , where $s=(q-1)/2$ (we are assuming $q$ is odd).
What are the roots of $x^s-1$ anyway?

In [27]:

print(sorted((x^s-1).roots()))

Out[27]:

[(1, 1), (3, 1), (4, 1), (5, 1), (9, 1), (12, 1), (13, 1), (14, 1), (15, 1), (16, 1), (19, 1), (20, 1), (22, 1), (25, 1), (27, 1), (34, 1), (36, 1), (39, 1), (41, 1), (42, 1), (45, 1), (46, 1), (47, 1), (48, 1), (49, 1), (52, 1), (56, 1), (57, 1), (58, 1), (60, 1)]

My keen sense of pattern recognition (or a recollection of Euler's criterion) tells me these are the roots of $x^s-1$ are precisely the squares (quadratic residues) in $\mathbb F_p^\times.$

In [28]:

product([x-a^2 for a in F if a < p-a]) == (x^s-1)

Out[28]:

True

But if $f$ has multiple roots that are quadratic residues (or maybe none), then what?

Rabin: Instead of $x^s-1$ , use $(x+\delta)^s-1$ for a random $\delta\in \mathbb F_q$ .

Definition:	Let us say that $\alpha,\beta\in \mathbb F_q$ are of different type if both are nonzero and $\alpha^s\ne \beta^s$ , in other words, one is a quadratic residue and one is not.
Theorem:	For every pair of distinct $\alpha,\beta\in \mathbb F_q$ we have $\#\{\delta\in \mathbb F_q:\alpha+\delta \text{ and } \beta+\delta \text{ are of different type }\} = \frac{q-1}{2}.$
Proof:	See Section 3.11 in the lecture notes for a short easy proof.
Corollary:	Let $f\in \mathbb F_q[x]$ , let $g(x)=\gcd(f,x^q-x)$ , let $s=(q-1)/2\in \Z$ , and let $\delta$ be a uniformly random element of $\mathbb F_q$ . If $g$ has more than one nonzero root, then $\Pr[0<\deg\gcd(g,(x+\delta)^s-1)<\deg g] \ge \frac{1}{2}.$

In [29]:

def rational_root(f):
    """Returns a rational root of f or None if f has no rational roots."""
    q = f.base_ring().cardinality()
    assert q % 2 == 1
    s = (q-1) // 2 # Note // yields an integer (whereas / would yield a rational)
    x = f.parent().gen()
    if f[0]==0:
        return f.base_ring()(0)
    g = distinct_root_poly(f)
    print("Initial g = gcd(f,x^q-x) = %s" % g)
    pig = g.parent().quotient(g)
    while g.degree() > 1:
        delta = g.base_ring().random_element()
        print("delta = %s" % delta)
        h = gcd(g, ((pig(x)-delta)^s-1).lift())
        if 0 < h.degree() < g.degree():
            g = h if 2*h.degree() <= g.degree() else g // h # as above, we use // to get a polynomial
            pig = g.parent().quotient(g)
            print("Better g = %s" % g)
    return -g[0] if g.degree() == 1 else None

Theorem:

Let

n=\log q

and

d=\deg f

. The expected running time of

\texttt{rational\_root}

\tilde O(n^2d)

bit operations.

Here the soft-O notation means up to a poly-logarithmic factors. See the lecture notes for a more precise bound and a proof.

In [30]:

for i in range(5):
    r = rational_root(f)
    print("Found root %s\n" % r)
    assert f(r) == 0

Out[30]:

Initial g = gcd(f,x^q-x) = x^3 + 24*x^2 + 33*x + 2
delta = 14
Better g = x + 17
Found root 44

Initial g = gcd(f,x^q-x) = x^3 + 24*x^2 + 33*x + 2
delta = 5
Better g = x + 46
Found root 15

Initial g = gcd(f,x^q-x) = x^3 + 24*x^2 + 33*x + 2
delta = 34
Better g = x + 17
Found root 44

Initial g = gcd(f,x^q-x) = x^3 + 24*x^2 + 33*x + 2
delta = 57
Better g = x + 22
Found root 39

Initial g = gcd(f,x^q-x) = x^3 + 24*x^2 + 33*x + 2
delta = 10
Better g = x + 22
Found root 39

In [31]:

h = product([x-i for i in range(1,20)])
for i in range(5):
    r = rational_root(h)
    print("Found root %s\n" % r)
    assert h(r) == 0

Out[31]:

Initial g = gcd(f,x^q-x) = x^19 + 54*x^18 + 40*x^17 + x^16 + 58*x^15 + 21*x^14 + 10*x^13 + 37*x^12 + 36*x^11 + 4*x^10 + 41*x^9 + 41*x^8 + 38*x^7 + 10*x^6 + 7*x^5 + 7*x^4 + 28*x^3 + 57*x^2 + 39*x + 19
delta = 25
Better g = x^9 + 37*x^8 + 31*x^7 + 7*x^6 + 11*x^5 + 54*x^4 + 56*x^3 + 16*x^2 + 34*x + 17
delta = 18
Better g = x^4 + 12*x^3 + 36*x^2 + 25*x + 14
delta = 13
Better g = x + 50
Found root 11

Initial g = gcd(f,x^q-x) = x^19 + 54*x^18 + 40*x^17 + x^16 + 58*x^15 + 21*x^14 + 10*x^13 + 37*x^12 + 36*x^11 + 4*x^10 + 41*x^9 + 41*x^8 + 38*x^7 + 10*x^6 + 7*x^5 + 7*x^4 + 28*x^3 + 57*x^2 + 39*x + 19
delta = 37
Better g = x^7 + 46*x^6 + 24*x^5 + 40*x^4 + 17*x^3 + 55*x^2 + 32*x + 29
delta = 19
Better g = x^3 + 31*x^2 + 50*x + 40
delta = 13
delta = 38
delta = 47
Better g = x + 60
Found root 1

Initial g = gcd(f,x^q-x) = x^19 + 54*x^18 + 40*x^17 + x^16 + 58*x^15 + 21*x^14 + 10*x^13 + 37*x^12 + 36*x^11 + 4*x^10 + 41*x^9 + 41*x^8 + 38*x^7 + 10*x^6 + 7*x^5 + 7*x^4 + 28*x^3 + 57*x^2 + 39*x + 19
delta = 56
Better g = x^9 + 27*x^8 + 38*x^7 + 55*x^6 + 19*x^5 + 48*x^4 + 52*x^3 + 36*x^2 + 54
delta = 52
Better g = x^4 + 29*x^3 + 3*x^2 + 32*x + 30
delta = 41
Better g = x + 54
Found root 7

Initial g = gcd(f,x^q-x) = x^19 + 54*x^18 + 40*x^17 + x^16 + 58*x^15 + 21*x^14 + 10*x^13 + 37*x^12 + 36*x^11 + 4*x^10 + 41*x^9 + 41*x^8 + 38*x^7 + 10*x^6 + 7*x^5 + 7*x^4 + 28*x^3 + 57*x^2 + 39*x + 19
delta = 2
Better g = x^8 + 45*x^7 + 55*x^6 + 44*x^5 + 22*x^4 + 49*x^3 + 47*x^2 + 60*x + 21
delta = 41
Better g = x^2 + 40*x + 38
delta = 35
Better g = x + 42
Found root 19

Initial g = gcd(f,x^q-x) = x^19 + 54*x^18 + 40*x^17 + x^16 + 58*x^15 + 21*x^14 + 10*x^13 + 37*x^12 + 36*x^11 + 4*x^10 + 41*x^9 + 41*x^8 + 38*x^7 + 10*x^6 + 7*x^5 + 7*x^4 + 28*x^3 + 57*x^2 + 39*x + 19
delta = 26
Better g = x^9 + 27*x^8 + 56*x^7 + 4*x^6 + 14*x^5 + 54*x^4 + 59*x^3 + 50*x^2 + 43*x + 40
delta = 23
Better g = x^4 + 23*x^3 + 40*x^2 + 23*x + 21
delta = 24
Better g = x + 45
Found root 16

In [32]:

p=2^61-1;  F=GF(p); R.<x>=PolynomialRing(F);  f=x^8-2*x+5; h = product([x-i for i in range(1,20)])
for i in range(3):
    %time print("Found root %s\n" % rational_root(h))

Out[32]:

Initial g = gcd(f,x^q-x) = x^19 + 2305843009213693761*x^18 + 16815*x^17 + 2305843009212773401*x^16 + 34916946*x^15 + 2305843008239752051*x^14 + 20692933630*x^13 + 2305842666961182051*x^12 + 4465226757381*x^11 + 2305796728565942041*x^10 + 381922055502195*x^9 + 2303339150458226401*x^8 + 12953636989943896*x^7 + 2253582105851181231*x^6 + 161429736530118960*x^5 + 1934458221868465951*x^4 + 610116075740491776*x^3 + 1637233278872540671*x^2 + 431565146817638400*x + 2184197908804861951
delta = 1405864291122341607
Better g = x^7 + 2305843009213693888*x^6 + 1660*x^5 + 2305843009213670221*x^4 + 198949*x^3 + 2305843009212714784*x^2 + 2621790*x + 2305843009210745551
delta = 1103017449276665091
Better g = x^3 + 2305843009213693922*x^2 + 255*x + 2305843009213693276
delta = 386511100125869182
Better g = x + 2305843009213693942
Found root 9

CPU times: user 11.7 ms, sys: 0 ns, total: 11.7 ms
Wall time: 11.1 ms
Initial g = gcd(f,x^q-x) = x^19 + 2305843009213693761*x^18 + 16815*x^17 + 2305843009212773401*x^16 + 34916946*x^15 + 2305843008239752051*x^14 + 20692933630*x^13 + 2305842666961182051*x^12 + 4465226757381*x^11 + 2305796728565942041*x^10 + 381922055502195*x^9 + 2303339150458226401*x^8 + 12953636989943896*x^7 + 2253582105851181231*x^6 + 161429736530118960*x^5 + 1934458221868465951*x^4 + 610116075740491776*x^3 + 1637233278872540671*x^2 + 431565146817638400*x + 2184197908804861951
delta = 532641917256937529
Better g = x^5 + 2305843009213693910*x^4 + 610*x^3 + 2305843009213689911*x^2 + 11584*x + 2305843009213682687
delta = 1312713933391314249
Better g = x^2 + 2305843009213693927*x + 128
delta = 181866373881774927
delta = 368701929104282100
Better g = x + 2305843009213693935
Found root 16

CPU times: user 12.4 ms, sys: 4 µs, total: 12.4 ms
Wall time: 12.1 ms
Initial g = gcd(f,x^q-x) = x^19 + 2305843009213693761*x^18 + 16815*x^17 + 2305843009212773401*x^16 + 34916946*x^15 + 2305843008239752051*x^14 + 20692933630*x^13 + 2305842666961182051*x^12 + 4465226757381*x^11 + 2305796728565942041*x^10 + 381922055502195*x^9 + 2303339150458226401*x^8 + 12953636989943896*x^7 + 2253582105851181231*x^6 + 161429736530118960*x^5 + 1934458221868465951*x^4 + 610116075740491776*x^3 + 1637233278872540671*x^2 + 431565146817638400*x + 2184197908804861951
delta = 1200605053011711393
Better g = x^8 + 2305843009213693874*x^7 + 2472*x^6 + 2305843009213650981*x^5 + 438989*x^4 + 2305843009211021938*x^3 + 9340578*x^2 + 2305843009196800051*x + 12004200
delta = 1073587239949998354
Better g = x^4 + 2305843009213693898*x^3 + 1009*x^2 + 2305843009213685884*x + 22230
delta = 984083528923832701
Better g = x + 2305843009213693932
Found root 19

CPU times: user 11.9 ms, sys: 0 ns, total: 11.9 ms
Wall time: 11.5 ms

What if we want all the rational roots?

In [33]:

def distinct_rational_roots(f,recurse=False):
    """Returns a list of all the distinct rational roots of f."""
    q = f.base_ring().cardinality()
    assert q % 2 == 1
    s = (q-1) // 2
    x = f.parent().gen()
    roots = []
    if recurse:
        g = f
    else:
        g = distinct_root_poly(f)
        if g(0) == 0:
            roots = [g.base_ring()(0)]
            g = g // x
    if g.degree() <= 1:
        return [-g[0]] if g.degree() == 1 else []
    while True:
        delta = g.base_ring().random_element()
        pig = g.parent().quotient(g)
        h = gcd(g, ((pig(x)-delta)^s-1).lift())
        if 0 < h.degree() < g.degree():
            roots += distinct_rational_roots(h, True)
            roots += distinct_rational_roots(g//h, True)
            return roots

Theorem:

Let

n=\log q

and

d=\deg f

. The expected running time of

\texttt{distinct\_rational\_root}

\tilde O(n^2d)

bit operations.

This looks the same as the bound for $\texttt{rational\_root}$ , but in fact there is an extra factor of $\log d$ hidden in the soft-O notation.

In [34]:

%time distinct_rational_roots(f)

Out[34]:

CPU times: user 4.19 ms, sys: 22 µs, total: 4.21 ms
Wall time: 4.17 ms

[1046353483181369856]

In [35]:

%time distinct_rational_roots(product([x-F.random_element() for i in range(10)]))

Out[35]:

CPU times: user 36.9 ms, sys: 0 ns, total: 36.9 ms
Wall time: 35.7 ms

[1358408883196038487,
 1548219673292526894,
 1813132387378194887,
 1737755998770690254,
 495898262673194171,
 631460305450662880,
 456341232373901345,
 1457055548390529253,
 462954065793106284,
 1445898814544229634]

What if we want all the rational roots and their multiplicities?

Yun: Compute the squarefree factorization!

Definition:	The squarefree factorization of a monic $f\in \mathbb F_q[x]$ is the unique list of monic squarefree polynomials $g_1,\ldots,g_m$ with $g_m\ne 1$ for which $f=g_1g_2^2\cdots g_m^m$ .
Key fact:	If $f=f_1\cdots f_n$ is the factorization of a monic $f\in \mathbb F_q[x]$ into irreducibles then $f_i=f_j$ for some $i\ne j$ if and only if $f_i\|\gcd(f,f')$ . This holds in any perfect field.

In [36]:

def squarefree_factorization(f):
    """Yun's algorithm to compute a list of polynomials g_1,..,g_m such that f=lc(f)g_1g_2^2...g_m^m with g_m != 1, assuming deg(f) < char(Fq).  We include g_0=1 for convenience."""
    assert f.degree() < f.base_ring().characteristic()
    if not f.is_monic:
        f = f.monic()
    df = f.derivative()
    u = gcd(f,df)
    v, w = f//u, df//u
    gs =[f.parent()(1)]
    while True:
        g = gcd(v, w-v.derivative()); gs.append(g)
        if v.degree() == g.degree():
            return gs
        v, w = v//g, (w-v.derivative())//g

Theorem:

Let

n=\log q

and

d=\deg f

. The running time of

\texttt{squarefree\_factorization}

\tilde O(nd)

bit operations.

This is negligible compared to the complexity of $\texttt{distinct\_rational\_roots}$ , so we might as well call $\texttt{squarefree\_factorization}$ first.
Even if we only want distinct roots, this will actually save time if $f$ is not squarefree.

In [37]:

h = product([(x^2+F.random_element()*x+F.random_element())^i for i in range(1,10)])
%time gs = squarefree_factorization(h)
assert h == product([gs[i]^i for i in range(len(gs))])
gs

Out[37]:

CPU times: user 2.19 ms, sys: 0 ns, total: 2.19 ms
Wall time: 1.98 ms

[1,
 x^2 + 1355173467251623258*x + 1735945799132775481,
 x^2 + 957635963191288544*x + 1860599609950326465,
 x^2 + 696390961420961863*x + 310902225905322874,
 x^2 + 176377541039972260*x + 216672670731559675,
 x^2 + 125144412738599419*x + 1371124668291939569,
 x^2 + 457549749373420095*x + 1073244427331095582,
 x^2 + 205794091840346165*x + 2119525877933536593,
 x^2 + 1005355233351380146*x + 2020926823563273184,
 x^2 + 1594893690500222301*x + 385338769927047306]

In [38]:

def rational_roots(f):
    """Returns a list of tuples identifying the rational roots of f and their multiplicities."""
    gs = squarefree_factorization(f)
    return reduce(lambda x, y: x+y,[[(a,i) for a in distinct_rational_roots(gs[i])] for i in range(1,len(gs))])

In [39]:

h = product([(x-F.random_element())^i for i in range(1,6)])
%time rational_roots(h)

Out[39]:

CPU times: user 11.9 ms, sys: 0 ns, total: 11.9 ms
Wall time: 11.5 ms

[(2107785688068906530, 1),
 (2091242374012486761, 2),
 (1578399650426587107, 3),
 (533312910932319669, 4),
 (725310156226374005, 5)]

What if we want to factor $f$ into irreducibles?

Rabin: You can do that with my root-finding algorithm.

Cantor-Zassenhaus: Yes, but that requires working in extension fields, and with our algorithm there is no need!

The first step is to compute the distinct degree factorization

Definition:	The distinct degree factorization of a monic $f\in \mathbb F_q[x]$ is the unique list of polynomials $g_1,\ldots,g_d$ with $f=g_1\cdots g_d$ such that each $g_i$ is a (possibly empty) product of distinct monic irreducible polynomials of degree $i$ , for $1\le i\le d=\deg f$ .
Key point:	This can be done deterministically by taking succesive gcds with $x^{q^i}-x$ .

In [40]:

def distinct_degree_poly(f,i):
    q = f.base_ring().cardinality()
    pif = f.parent().quotient(f)
    return gcd(f,(pif(x)^(q^i)).lift()-x)

def distinct_degree_factorization(f, squarefree=False):
    """
        Computes the distinct degree factorization of f as list gs with gs[0]=1 and gs[i]=g_i such that f=prod_i gs[i]^i.
        The optional parameter squarefree inidicates whether f is known to be squarefree (usually one computes the squarefree factorization first).
    """
    if not squarefree:
        return naive_distinct_degree_factorization(f)
    d = f.degree()
    gs = [f.parent()(1) for i in range(d+1)]
    for i in range(1,d+1):
        if 2*i > f.degree():
            gs[f.degree()] = f
            break
        gs[i] = distinct_degree_poly(f,i)
        f = f // gs[i]
    return gs

def naive_distinct_degree_factorization(f):
    """
        Computed distinct degree factorization of arbitrary f by combining distinct degree factorization of polys in square free factorization.
        There is no reason to ever use this, it is provided just for the sake of completeness
    """
    d = f.degree()
    gs = [f.parent()(1) for i in range(d+1)]
    hs = [distinct_degree_factorization(g, squarefree=True) for g in squarefree_factorization(f)]
    # combine results (this is silly, we are undoing work here)
    for i in range(1,len(hs)):
        for j in range(1,len(hs[i])):
            gs[j] *= hs[i][j]^i
    return gs

def facpat(f):
    """Returns a dictionary {d:i} whose keys are the degrees of the irreducible factors of f and whose values are the number of irreducible factors of that degree."""
    X = {}

    gs = squarefree_factorization(f)
    for i in range(1,len(gs)):
        gis = distinct_degree_factorization(gs[i],squarefree=True)
        for j in range(1,len(gis)):
            if gis[j].degree() > 0:
                X[j] = X.get(j,0) + i * (gis[j].degree() // j)
    return X

In [41]:

%time facpat(f)

Out[41]:

CPU times: user 7.33 ms, sys: 0 ns, total: 7.33 ms
Wall time: 7.26 ms

{1: 1, 2: 2, 3: 1}

In [42]:

h = (x^9-1)^2*(x^32-1)
assert h == product(distinct_degree_factorization(h))
%time print(facpat(h))
print([(r[0].degree(),r[1]) for r in h.factor()])

Out[42]:

{1: 20, 2: 15}
CPU times: user 4.31 ms, sys: 36 µs, total: 4.35 ms
Wall time: 4.26 ms
[(1, 1), (1, 2), (1, 2), (1, 2), (1, 2), (1, 2), (1, 2), (1, 2), (1, 2), (1, 3), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1), (2, 1)]

The key innovation introduced by Cantor-Zassenhaus is that when searching for irreducible factors $g$ of $f$ of degree $j$ , rather than computing $\gcd(f,(x+\delta)^s-1))$ using a random $\delta\in \mathbb F_{q^j}$ and computing the gcd in $\mathbb F_{q^j}[x]$ (which is how one would go about finding a root of $g$ using Rabin's algorithm), it is better to compute $\gcd(f,u^s-1)$ in $\mathbb F_q[x]$ using a random $u\in \mathbb F_q[x]$ coprime to $f$ with $\deg u < \deg f$ . This faster because we are working the ring $\mathbb F_q[x]/(f)$ rather than $\mathbb F_{q^j}[x]/(f)$ (so the bit-size of each element is smaller by a factor of $j$ ), and the probability of success is the same.

Theorem:	Let $f\in \mathbb F_q[x]$ be the product of $r\ge 2$ distinct monic irreducible polynomials of degree $j$ and let $u\in \mathbb F_q[j]$ be a random polynomial coprime to $f$ with $\deg u < \deg f$ . Then $\Pr[0<\deg\gcd(f,u^s-1)<\deg f] = 2^{1-r} \ge \frac{1}{2}.$
Proof:	This follows from applying the Chinese Remainder Theorem to $\mathbb F_q[x]/(f)\simeq \mathbb F_{q^j}^r$ . See Section 3.12 of the lecture notes for details.

In [43]:

def equal_degree_factorization(f,j):
    """Completely factors a monic polynomial f in Fq[x] that is known to be the product of distinct irreducible polynomials of degree j"""
    d = f.degree()
    assert d % j == 0
    if d == 0:
        return []
    elif d == j:
        return [f]
    Fq = f.base_ring(); q = Fq.cardinality()
    assert q%2 == 1
    s = (q^j-1) // 2
    x = f.parent().gen()
    pif = f.parent().quotient(f)
    while True:
        # generate a random polynomial of degree < deg(f), note f.parent()(v) converts a vector of Fq coefficients to a polynomial
        u = f.parent()([Fq.random_element() for i in range(d)])
        h = gcd(f,u)
        if 0 < h.degree() < d:
            return equal_degree_factorization(h,j) + equal_degree_factorization(f//h,j)
        # at this point we know u is coprime to f
        h = gcd(f,(pif(u)^s-1).lift())
        if 0 < h.degree() < d:
            return equal_degree_factorization(h,j) + equal_degree_factorization(f//h,j)

def factorization(f):
    """Computes the complete factorization of a non-constant f in Fq[x] using the Cantor-Zassenhaaus algorithm (assumes q is odd)."""
    assert f.degree() > 0
    if not f.is_monic:
        f = f.monic()
    x = f.parent().gen()
    res = []
    gs = squarefree_factorization(f)
    for i in range(1,len(gs)):
        gis = distinct_degree_factorization(gs[i],squarefree=True)
        for j in range(1,len(gis)):
            res += [(g,i) for g in equal_degree_factorization(gis[j],j)]
    # sort output to make it canonical (independent of random choices)
    return sorted(res)

In [44]:

assert factorization(f) == sorted(list(f.factor()))
print(factorization(f))
%timeit -n 10 factorization(f)
%timeit -n 10 f.factor()

Out[44]:

[(x + 1259489526032324095, 1), (x^2 + 285383323066146333*x + 1665318937494081999, 1), (x^2 + 794947353545739013*x + 1355389809259494196, 1), (x^3 + 2271865815783178461*x^2 + 2223302821568905349*x + 726924123724355165, 1)]
6.47 ms ± 2.73 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
309 µs ± 52.8 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [45]:

assert factorization(h) == sorted(list((h).factor()))
print(factorization(h))
%timeit -n 1 -r 5 factorization(h)
%timeit -n 1 -r 5 h.factor()

Out[45]:

[(x + 1, 1), (x + 52855892431037422, 2), (x + 541562117690345921, 2), (x + 636260618972345636, 2), (x + 1202998424213388074, 2), (x + 1669582390241348316, 2), (x + 1672775772132701819, 2), (x + 1711424999092310608, 2), (x + 1735911822081298009, 2), (x + 2305843009213693950, 3), (x^2 + 1, 1), (x^2 + 2147483648*x + 1, 1), (x^2 + 44054674105924332*x + 1, 1), (x^2 + 166609066672189134*x + 1, 1), (x^2 + 658669255311844864*x + 1, 1), (x^2 + 826015582034181567*x + 1, 1), (x^2 + 911195041853779459*x + 1, 1), (x^2 + 1080792493261747995*x + 1, 1), (x^2 + 1225050515951945956*x + 1, 1), (x^2 + 1394647967359914492*x + 1, 1), (x^2 + 1479827427179512384*x + 1, 1), (x^2 + 1647173753901849087*x + 1, 1), (x^2 + 2139233942541504817*x + 1, 1), (x^2 + 2261788335107769619*x + 1, 1), (x^2 + 2305843007066210303*x + 1, 1)]
57.7 ms ± 10.7 ms per loop (mean ± std. dev. of 5 runs, 1 loop each)
3.15 ms ± 69.1 µs per loop (mean ± std. dev. of 5 runs, 1 loop each)

Theorem:

Let

n=\log q

and

d=\deg f

. The expected running time of

\texttt{factorization}

\tilde O(n^2d^2)

bit operations.

There are other factorization methods based on linear algebra and modular composition that improve the dependence on $d$ .
Currently the best asymptotic bound is achieved by the Kedlaya-Umans algorithm with an expected running time of $\tilde O(d^{3/2}n + dn^2)$ .

In [46]:

p=2^255-19;  F=GF(p); R.<x>=PolynomialRing(F);  f=x^8-2*x+5; h = product([x-i for i in range(1,20)])
assert factorization(f^3+f^2+1) == sorted(list((f^3+f^2+1).factor()))
%time print(factorization(f))

Out[46]:

[(x + 51027038539503343326764519138825597294378744664834164009924300106595537523144, 1), (x^2 + 50611585019097526067500995121554429372177258784480155764381858723769139594768*x + 38688838658452246841078006525201819047008644554995237691159504617557701189662, 1), (x^5 + 14153465678715326029305470748307881186713981216326244265151425177548452521986*x^4 + 30970415661099139061786692430821851662201464041579196405043008572374691402945*x^3 + 499132331909190267118163511147975581363888184276674106555956712514448901935*x^2 + 8376079156473006437715373698581520278472986811325551025499711711353793595153*x + 26683694078622147531159282359840405681456757722151517437587006163097167872535, 1)]
CPU times: user 170 ms, sys: 182 µs, total: 170 ms
Wall time: 170 ms

In [0]:

Product

Resources

Company