3 files changed, 199 insertions, 69 deletions

diff --git a/pyecsca/ec/divpoly.py b/pyecsca/ec/divpoly.py
index 6621b76..63024cb 100644
--- a/pyecsca/ec/divpoly.py
+++ b/pyecsca/ec/divpoly.py
@@ -11,6 +11,14 @@ from .curve import EllipticCurve
 from .mod import Mod
 from .model import ShortWeierstrassModel
 
+has_pari = False
+try:
+    import cypari2
+
+    has_pari = True
+except ImportError:
+    cypari2 = None
+
 
 def values(*ns: int) -> Mapping[int, Tuple[int, ...]]:
     done: Set[int] = set()
@@ -102,10 +110,12 @@ def b_invariants(curve: EllipticCurve) -> Tuple[Mod, ...]:
     """
     if isinstance(curve.model, ShortWeierstrassModel):
         a1, a2, a3, a4, a6 = a_invariants(curve)
-        return (a1 * a1 + 4 * a2,
-                a1 * a3 + 2 * a4,
-                a3 ** 2 + 4 * a6,
-                a1 ** 2 * a6 + 4 * a2 * a6 - a1 * a3 * a4 + a2 * a3 ** 2 - a4 ** 2)
+        return (
+            a1 * a1 + 4 * a2,
+            a1 * a3 + 2 * a4,
+            a3**2 + 4 * a6,
+            a1**2 * a6 + 4 * a2 * a6 - a1 * a3 * a4 + a2 * a3**2 - a4**2,
+        )
     else:
         raise NotImplementedError
 
@@ -142,7 +152,7 @@ def divpoly0(curve: EllipticCurve, *ns: int) -> Mapping[int, Poly]:
         if i == -2:
             val = mem[-1] ** 2
         elif i == -1:
-            val = Kx(4) * x ** 3 + b2 * x ** 2 + Kx(2) * b4 * x + b6
+            val = Kx(4) * x**3 + b2 * x**2 + Kx(2) * b4 * x + b6
         elif i == 0:
             val = Kx(0)
         elif i < 0:
@@ -150,9 +160,11 @@ def divpoly0(curve: EllipticCurve, *ns: int) -> Mapping[int, Poly]:
         elif i in (1, 2):
             val = Kx(1)
         elif i == 3:
-            val = Kx(3) * x ** 4 + b2 * x ** 3 + Kx(3) * b4 * x ** 2 + Kx(3) * b6 * x + b8
+            val = (
+                Kx(3) * x**4 + b2 * x**3 + Kx(3) * b4 * x**2 + Kx(3) * b6 * x + b8
+            )
         elif i == 4:
-            val = -mem[-2] + (Kx(6) * x ** 2 + b2 * x + b4) * mem[3]
+            val = -mem[-2] + (Kx(6) * x**2 + b2 * x + b4) * mem[3]
         elif i % 2 == 0:
             m = (i - 2) // 2
             val = mem[m + 1] * (mem[m + 3] * mem[m] ** 2 - mem[m - 1] * mem[m + 2] ** 2)
@@ -203,7 +215,9 @@ def divpoly(curve: EllipticCurve, n: int, two_torsion_multiplicity: int = 2) ->
 
 
 @public
-def mult_by_n(curve: EllipticCurve, n: int, x_only: bool = False) -> Tuple[Tuple[Poly, Poly], Optional[Tuple[Poly, Poly]]]:
+def mult_by_n(
+    curve: EllipticCurve, n: int, x_only: bool = False
+) -> Tuple[Tuple[Poly, Poly], Optional[Tuple[Poly, Poly]]]:
     """
     Compute the multiplication-by-n map on an elliptic curve.
 
@@ -263,18 +277,43 @@ def mult_by_n(curve: EllipticCurve, n: int, x_only: bool = False) -> Tuple[Tuple
     mxd_full_denom = mxd_dn_denom
 
     # > a1*mx
-    a1mx_num = (Kxy(a1) * mx[0])
+    a1mx_num = Kxy(a1) * mx[0]
     a1mx_denom = mx[1]  # noqa
 
     # > a3
-    a3_num = (Kxy(a3) * mx[1])
+    a3_num = Kxy(a3) * mx[1]
     a3_denom = mx[1]  # noqa
 
     # The mx.derivative part has a different denominator, basically mx[1]^2 * m
-    # so the rest needs to be multiplied by this factor when subtracitng.
+    # so the rest needs to be multiplied by this factor when subtracting.
     mxd_fact = mx[1] * n
 
-    my_num = (mxd_full_num - a1mx_num * mxd_fact - a3_num * mxd_fact)
+    my_num = mxd_full_num - a1mx_num * mxd_fact - a3_num * mxd_fact
     my_denom = mxd_full_denom * Kxy(2)
     my = (my_num, my_denom)
     return mx, my
+
+
+if has_pari:
+
+    def mult_by_n_pari(curve: EllipticCurve, order: int, n: int):
+        pari = cypari2.Pari()
+        p = pari(curve.prime)
+        a = pari.Mod(curve.parameters["a"], p)
+        b = pari.Mod(curve.parameters["b"], p)
+        E = pari.ellinit([a, b])
+        E[15][0] = pari(order)
+        while True:
+            try:
+                mx = pari.ellxn(E, n)
+                break
+            except cypari2.PariError as e:
+                if e.errnum() == 17:  # out of stack memory
+                    pari.allocatemem(0)
+                else:
+                    raise e
+        x = symbols("x")
+        K = FF(curve.prime)
+        mx_num = Poly([int(coeff) for coeff in reversed(mx[0])], x, domain=K)
+        mx_denom = Poly([int(coeff) for coeff in reversed(mx[1])], x, domain=K)
+        return mx_num, mx_denom
diff --git a/pyproject.toml b/pyproject.toml
index 48ad351..089c6c2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -61,6 +61,7 @@
 "smartcard" = ["pyscard"]
 "leia" = ["smartleia"]
 "gmp" = ["gmpy2"]
+"pari" = ["cysignals", "cypari2"]
 "dev" = ["mypy", "flake8", "interrogate", "pyinstrument", "black", "types-setuptools"]
 "test" = ["pytest>=7.0.0", "coverage", "pytest-cov", "pytest-sugar"]
 "doc" = ["sphinx", "sphinx-autodoc-typehints", "nbsphinx", "sphinx-paramlinks", "sphinx_design"]
diff --git a/test/ec/test_divpoly.py b/test/ec/test_divpoly.py
index 07eed76..065e1ca 100644
--- a/test/ec/test_divpoly.py
+++ b/test/ec/test_divpoly.py
@@ -1,4 +1,6 @@
 import json
+
+import pytest
 from importlib_resources import files
 
 import test.data.divpoly
@@ -8,47 +10,90 @@ from pyecsca.ec.divpoly import a_invariants, b_invariants, divpoly0, divpoly, mu
 
 def test_ainvs(secp128r1):
     ainvs = a_invariants(secp128r1.curve)
-    assert ainvs == (0, 0, 0, 340282366762482138434845932244680310780, 308990863222245658030922601041482374867)
+    assert ainvs == (
+        0,
+        0,
+        0,
+        340282366762482138434845932244680310780,
+        308990863222245658030922601041482374867,
+    )
 
 
 def test_binvs(secp128r1):
     binvs = b_invariants(secp128r1.curve)
-    assert binvs == (0, 340282366762482138434845932244680310777, 215116352601536216819152607431888567119,
-                     340282366762482138434845932244680310774)
+    assert binvs == (
+        0,
+        340282366762482138434845932244680310777,
+        215116352601536216819152607431888567119,
+        340282366762482138434845932244680310774,
+    )
 
 
 def test_divpoly0(secp128r1):
     # Data from sagemath
-    coeffs = [11, 0, 340282366762482138434845932244680302401, 211962053797180672439257756222135086642,
-              340282366762482138434845932244678441564, 115415922367823003571854983213102698477,
-              152803211743444076787231275062278784385, 68540219804769369063918923691867278088,
-              43207172520353703997069627419519708522, 83208285732019037267730920881743782729,
-              93286967763556583502947234289842152563, 324950611928652823046744874201355360259,
-              244242343224213805514200367379671854852, 307096814154284337284845014037169929735,
-              180946781765592277412990188457219828893, 301253861469456022084288029442105687698,
-              58053323975526190296189278379252064657, 224437885189054146208302696540070489578,
-              281987318191429654256483850017931541622, 21449216018131966691124843738286677726,
-              10958264881628724646042625283328121348, 104868338562600481545003572552335444641,
-              127205813185570107009206143413997395181, 116865717360861207318274706645935808417,
-              281460458922812844939222119784601506753, 336607098463310980140968249747513775735,
-              304486486784143285234063826161805094682, 194935097339732797131694429642153881938,
-              193523171473792085604518744912658246509, 204844449336357293979832621297234119270,
-              244481753281744913785581086721299830802, 46816299473081369405217767361380254657,
-              303070923752707405164354702252828590781, 222516549119176621389776816552836322766,
-              292006660232236762950883960515487362063, 53617127992846936725441702182362940200,
-              242498306026562585655027965022211017540, 25039963304689451659955607939868533124,
-              328580435950647191774558154445103295305, 24226614081978788956695324769468902511,
-              147945052666123617872720080832548744564, 287190187011075399698210761813202261601,
-              117131681517270554750959286838283723521, 35018410385280384289320020556813474742,
-              83939964512240352730304831725346032711, 147219996946006689656600631222993527180,
-              280430477096741745234510250577626566690, 32753113267385981127807026368593329576,
-              105134319561523011785486683031223863934, 206456116679151691099661865534540095270,
-              116180470443213022739312068090342951131, 245850120846480965440408943459023315919,
-              45805943896736805301879725516256422457, 226777421435695229777151315574975350291,
-              283680841707610526659029980964566557627, 53168487339451866167506032177471934158,
-              69212302225932892622760219621519562036, 183916411340675637978873336955593385541,
-              119478537598919956688656337369481692789, 234767298887335988751880131162396819780,
-              218412162101425422347176804186940045781]
+    coeffs = [
+        11,
+        0,
+        340282366762482138434845932244680302401,
+        211962053797180672439257756222135086642,
+        340282366762482138434845932244678441564,
+        115415922367823003571854983213102698477,
+        152803211743444076787231275062278784385,
+        68540219804769369063918923691867278088,
+        43207172520353703997069627419519708522,
+        83208285732019037267730920881743782729,
+        93286967763556583502947234289842152563,
+        324950611928652823046744874201355360259,
+        244242343224213805514200367379671854852,
+        307096814154284337284845014037169929735,
+        180946781765592277412990188457219828893,
+        301253861469456022084288029442105687698,
+        58053323975526190296189278379252064657,
+        224437885189054146208302696540070489578,
+        281987318191429654256483850017931541622,
+        21449216018131966691124843738286677726,
+        10958264881628724646042625283328121348,
+        104868338562600481545003572552335444641,
+        127205813185570107009206143413997395181,
+        116865717360861207318274706645935808417,
+        281460458922812844939222119784601506753,
+        336607098463310980140968249747513775735,
+        304486486784143285234063826161805094682,
+        194935097339732797131694429642153881938,
+        193523171473792085604518744912658246509,
+        204844449336357293979832621297234119270,
+        244481753281744913785581086721299830802,
+        46816299473081369405217767361380254657,
+        303070923752707405164354702252828590781,
+        222516549119176621389776816552836322766,
+        292006660232236762950883960515487362063,
+        53617127992846936725441702182362940200,
+        242498306026562585655027965022211017540,
+        25039963304689451659955607939868533124,
+        328580435950647191774558154445103295305,
+        24226614081978788956695324769468902511,
+        147945052666123617872720080832548744564,
+        287190187011075399698210761813202261601,
+        117131681517270554750959286838283723521,
+        35018410385280384289320020556813474742,
+        83939964512240352730304831725346032711,
+        147219996946006689656600631222993527180,
+        280430477096741745234510250577626566690,
+        32753113267385981127807026368593329576,
+        105134319561523011785486683031223863934,
+        206456116679151691099661865534540095270,
+        116180470443213022739312068090342951131,
+        245850120846480965440408943459023315919,
+        45805943896736805301879725516256422457,
+        226777421435695229777151315574975350291,
+        283680841707610526659029980964566557627,
+        53168487339451866167506032177471934158,
+        69212302225932892622760219621519562036,
+        183916411340675637978873336955593385541,
+        119478537598919956688656337369481692789,
+        234767298887335988751880131162396819780,
+        218412162101425422347176804186940045781,
+    ]
     K = FF(secp128r1.curve.prime)
     poly = divpoly0(secp128r1.curve, 11)[11]
     computed = list(map(K, poly.all_coeffs()))
@@ -58,18 +103,34 @@ def test_divpoly0(secp128r1):
 def test_divpoly(secp128r1):
     # Data from sagemath
     K = FF(secp128r1.curve.prime)
-    coeffs_0 = {(0,): K(16020440675387382717114730680672549016), (1,): K(269851015321770885610377847857290470365),
-                (2,): K(340282366762482138434845932244680310693), (3,): K(109469325440469337582450480850803806492),
-                (4,): K(340282366762482138434845932244680310753), (6,): K(2)}
+    coeffs_0 = {
+        (0,): K(16020440675387382717114730680672549016),
+        (1,): K(269851015321770885610377847857290470365),
+        (2,): K(340282366762482138434845932244680310693),
+        (3,): K(109469325440469337582450480850803806492),
+        (4,): K(340282366762482138434845932244680310753),
+        (6,): K(2),
+    }
     assert divpoly(secp128r1.curve, 4, 0).as_dict() == coeffs_0
-    coeffs_1 = {(6, 1): K(4), (4, 1): K(340282366762482138434845932244680310723),
-                (3, 1): K(218938650880938675164900961701607612984), (2, 1): K(340282366762482138434845932244680310603),
-                (1, 1): K(199419663881059632785909763469900629947), (0, 1): K(32040881350774765434229461361345098032)}
+    coeffs_1 = {
+        (6, 1): K(4),
+        (4, 1): K(340282366762482138434845932244680310723),
+        (3, 1): K(218938650880938675164900961701607612984),
+        (2, 1): K(340282366762482138434845932244680310603),
+        (1, 1): K(199419663881059632785909763469900629947),
+        (0, 1): K(32040881350774765434229461361345098032),
+    }
     assert divpoly(secp128r1.curve, 4, 1).as_dict() == coeffs_1
-    coeffs_2 = {(9,): K(8), (7,): K(340282366762482138434845932244680310639),
-                (6,): K(187545273439985507098415273777631738640), (4,): K(117928913205007755574446043156465405646),
-                (3,): K(244159722710157842132157548160645018307), (2,): K(200234655086793134086408617236124137371),
-                (1,): K(51914434605509249526780779992574428819), (0,): K(60581150995923875019702403440670701629)}
+    coeffs_2 = {
+        (9,): K(8),
+        (7,): K(340282366762482138434845932244680310639),
+        (6,): K(187545273439985507098415273777631738640),
+        (4,): K(117928913205007755574446043156465405646),
+        (3,): K(244159722710157842132157548160645018307),
+        (2,): K(200234655086793134086408617236124137371),
+        (1,): K(51914434605509249526780779992574428819),
+        (0,): K(60581150995923875019702403440670701629),
+    }
     assert divpoly(secp128r1.curve, 4, 2).as_dict() == coeffs_2
 
 
@@ -77,16 +138,28 @@ def test_mult_by_n(secp128r1):
     # Data from sagemath
     K = FF(secp128r1.curve.prime)
     coeffs_mx_num = [1, 0, 6, 250332028321891843231386649625583487328, 9]
-    coeffs_mx_denom = [4, 0, 340282366762482138434845932244680310771, 215116352601536216819152607431888567119]
-    coeffs_my_num = {(6, 1): K(8), (4, 1): K(340282366762482138434845932244680310663),
-                     (3, 1): K(97594934999395211894955991158534915185),
-                     (2, 1): K(340282366762482138434845932244680310423),
-                     (1, 1): K(58556960999637127136973594695120949111),
-                     (0, 1): K(64081762701549530868458922722690196064)}
-    coeffs_my_denom = {(6, 0): K(64), (4, 0): K(340282366762482138434845932244680310399),
-                       (3, 0): K(78075947999516169515964792926827932148), (2, 0): K(576),
-                       (1, 0): K(106054522763933629886951553464196514339),
-                       (0, 0): K(276200604060932607566387009521990114935)}
+    coeffs_mx_denom = [
+        4,
+        0,
+        340282366762482138434845932244680310771,
+        215116352601536216819152607431888567119,
+    ]
+    coeffs_my_num = {
+        (6, 1): K(8),
+        (4, 1): K(340282366762482138434845932244680310663),
+        (3, 1): K(97594934999395211894955991158534915185),
+        (2, 1): K(340282366762482138434845932244680310423),
+        (1, 1): K(58556960999637127136973594695120949111),
+        (0, 1): K(64081762701549530868458922722690196064),
+    }
+    coeffs_my_denom = {
+        (6, 0): K(64),
+        (4, 0): K(340282366762482138434845932244680310399),
+        (3, 0): K(78075947999516169515964792926827932148),
+        (2, 0): K(576),
+        (1, 0): K(106054522763933629886951553464196514339),
+        (0, 0): K(276200604060932607566387009521990114935),
+    }
     mx, my = mult_by_n(secp128r1.curve, 2)
     mx_num, mx_denom = mx
     assert coeffs_mx_num == list(map(K, mx_num.all_coeffs()))
@@ -101,12 +174,29 @@ def test_mult_by_n_large(secp128r1):
     mx, my = mult_by_n(secp128r1.curve, 21)
     with files(test.data.divpoly).joinpath("mult_21.json").open("r") as f:
         sage_data = json.load(f)
-        sage_data["mx"][0] = {eval(key): K(val) for key, val in sage_data["mx"][0].items()}
-        sage_data["mx"][1] = {eval(key): K(val) for key, val in sage_data["mx"][1].items()}
-        sage_data["my"][0] = {eval(key): K(val) for key, val in sage_data["my"][0].items()}
-        sage_data["my"][1] = {eval(key): K(val) for key, val in sage_data["my"][1].items()}
+        sage_data["mx"][0] = {
+            eval(key): K(val) for key, val in sage_data["mx"][0].items()
+        }
+        sage_data["mx"][1] = {
+            eval(key): K(val) for key, val in sage_data["mx"][1].items()
+        }
+        sage_data["my"][0] = {
+            eval(key): K(val) for key, val in sage_data["my"][0].items()
+        }
+        sage_data["my"][1] = {
+            eval(key): K(val) for key, val in sage_data["my"][1].items()
+        }
 
         assert mx[0].as_dict() == sage_data["mx"][0]
         assert mx[1].as_dict() == sage_data["mx"][1]
         assert my[0].as_dict() == sage_data["my"][0]
         assert my[1].as_dict() == sage_data["my"][1]
+
+
+def test_mult_by_n_pari(secp128r1):
+    _ = pytest.importorskip("cypari2")
+    from pyecsca.ec.divpoly import mult_by_n_pari
+
+    mx_pari = mult_by_n_pari(secp128r1.curve, secp128r1.order, 10)
+    mx_our, _ = mult_by_n(secp128r1.curve, 10, x_only=True)
+    assert mx_pari == mx_our