Hensel's Lifting Method

近日尝试用sagemath自带的roots()求解\(f(x) \equiv 0 \pmod{p^a}\)，发现速度极慢。了解到Hensel's Lemma可以解决这个问题，但我在sagemath的申必文档中并未搜索到该方法的直接实现。故自己写了一个。

简述

Hensel's lemma

Hensel’s Lemma 用于将 \(f(x) \equiv 0 \pmod {p}\) 的解提升至\(Zmod(p^a)\). (p为素数)


下面的算法称之为 Hensel’s lifting method :

首先，已知\(x \equiv r_{k-1}\) 为 \(f(x) \equiv 0 \pmod {p^{k-1}}\) 的解。

接下来将数域由\(GF(p^{k-1})\)提升至\(GF(p^k)\)，记 \(r_k\) 为 \(f(x) \equiv 0 \pmod {p^{k}}\) 的解

  • 若 \(f'(r_{k-1}) \not\equiv 0\), 则有唯一的\(t\)使得 \(f(r_{k-1} + tp^{k-1}) \equiv 0 \pmod{p^k}\) 成立，此时\(r_k = r_{k-1} - \displaystyle\frac{f(r_{k-1})}{f'(r_{k-1})}\)
  • 若 \(f'(r_{k-1}) \equiv 0\), 且 \(f(r_{k-1}) \equiv 0\)，则t为任意整数，\(r_k = r_{k-1}\)
  • 若 \(f'(r_{k-1}) \equiv 0\), 且 \(f(r_{k-1}) \not\equiv 0\)，则没有满足条件的整数t，无法利用\(r_{k-1}\)求解

以上步骤循环，可以将\(f(x) \equiv 0\) 在 \(GF(p^a)\)上的求解规约到 \(GF(p)\)


注意到\(\displaystyle\frac{f(r_{k-1})}{f'(r_{k-1})}\)未必能直接产生整数，将计算放到p-adic ring上更方便。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def hensel_solve(f, p, r):
    """
    Solves polynomial roots in the ring Zmod(p**r) using Hensel's lifting method.

    Parameters:
    f (polynomial): The polynomial equation.
    p (int): A prime number.
    r (int): The exponent.
    
    Raises:
    ValueError: If p is not a prime number or if f has no roots.
    """
    if not is_prime(p):
        raise ValueError("p must be a prime")
    f = f.change_ring(Zp(p))
    F = f.change_ring(Zmod(pow(p,r)))
    P = Zp(p,max(30, r))
    Fd = derivative(F)
    origin_roots = f.roots()
    if not len(origin_roots):
        raise ValueError("f has no roots")
    ans = set()
    for x in origin_roots:
        x_k = ZZ(x[0])
        flag = 0
        for k in range(1,r):
            if Fd(x_k) == P(0):
                if Zmod(pow(p,r))(f(x_k)) == 0:
                    continue
                else:
                    flag = 1
                    break
            else:
                x_k = Zmod(pow(p,r))(P(x_k)-P(F(x_k))/P(Fd(x_k)))
        if not flag:
            ans.update({x_k})
    return list(ans)

发布在了pypi上，包名是sageball，以后有空继续加东西 :)

用例：

1
pip install sageball
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
sage: from sageball import *
sage: from Crypto.Util.number import *
sage: p = getPrime(256)
sage: R.<y> = PolynomialRing(Zmod(p**20))
sage: f = y^3 + 88*y^2 - 99999
sage: f.change_ring(Zmod(p)).roots()
[(69308405488968049481748767212436730722801895792767406183685348133444233479868,
  1)]
sage: hensel_solve(f, p, 20)
[37703239737291277027365421519817323510865478112086210957276760542087052163016132599406943585036792976177249575532822533382718659736377573208253677268888201655419242055704276470729367796428426925172977074351292389204307294573364430803709392333960546348414634462190043821875241702966811101872263679102636643307528125089896228460596363355183960529903004848582857601840357174398281039729957164593085194467208379329098387709805523851043891922678228619680505251195773703511681160617319758877015996276560753899603002631615128325997643360635074115199442764683911096103905062585582942055722166181353528675325076096297244626026449521382868966037541063655163348755341756755617940070562452441891408449830499203414189319262630366583841725437180190345364849360870085927141091486810600074199897251843646588089235927963654342587841938352499246682037892369415183228114192308636949277636053376129783325911050681522271118165681832491205844031365819303630911196836805570188269207857049649045106759562098511078380779871363209078702905615571385458121535111843353457641165749407944240106157067093368333261474712361323801567878571827355353862484780988204946203231824057103525010377318719780405547502696530468817018627278660067985068655748252348692956763839842442577386154738340668457963492309894080502092443676751298594398480522893553162155245967675698842596912463075386927752665256904505233043287372778132982200372598578779053114623191692012462398153277118305607887829172213095033538161975762360571566038379832338161703977843510412337910651837461493861836641371994]
sage: f(3770323973729127702736542151981732351086547811208621095727676054208705216301613259940694358503679297617724957553
....: 282253338271865973637757320825367726888820165541924205570427647072936779642842692517297707435129238920430729457336
....: 443080370939233396054634841463446219004382187524170296681110187226367910263664330752812508989622846059636335518396
....: 052990300484858285760184035717439828103972995716459308519446720837932909838770980552385104389192267822861968050525
....: 119577370351168116061731975887701599627656075389960300263161512832599764336063507411519944276468391109610390506258
....: 558294205572216618135352867532507609629724462602644952138286896603754106365516334875534175675561794007056245244189
....: 140844983049920341418931926263036658384172543718019034536484936087008592714109148681060007419989725184364658808923
....: 592796365434258784193835249924668203789236941518322811419230863694927763605337612978332591105068152227111816568183
....: 249120584403136581930363091119683680557018826920785704964904510675956209851107838077987136320907870290561557138545
....: 812153511184335345764116574940794424010615706709336833326147471236132380156787857182735535386248478098820494620323
....: 182405710352501037731871978040554750269653046881701862727866006798506865574825234869295676383984244257738615473834
....: 066845796349230989408050209244367675129859439848052289355316215524596767569884259691246307538692775266525690450523
....: 304328737277813298220037259857877905311462319169201246239815327711830560788782917221309503353816197576236057156603
....: 8379832338161703977843510412337910651837461493861836641371994)
0