Exploration of formulas in open-source ECC libraries

import itertools
import tabulate
import pickle
import numpy as np
import inspect
import tempfile
import sys
import multiprocessing

from copy import deepcopy
from concurrent.futures import ProcessPoolExecutor, as_completed
from contextlib import contextmanager
from importlib import import_module, invalidate_caches
from pathlib import Path

from matplotlib import pyplot as plt
from IPython.display import HTML, display
from tqdm.notebook import tqdm
from importlib_resources import files, as_file

import pyecsca
from pyecsca.ec.model import ShortWeierstrassModel, TwistedEdwardsModel, MontgomeryModel
from pyecsca.ec.formula import AdditionEFDFormula, DoublingEFDFormula, LadderEFDFormula
from pyecsca.ec.params import get_params
from pyecsca.ec.formula.metrics import formula_similarity, formula_similarity_fuzz
from pyecsca.ec.formula.unroll import unroll_formula_expr, unroll_formula
from pyecsca.ec.formula.expand import expand_formula_set


# Allow to use "spawn" multiprocessing method for function defined in a Jupyter notebook.
# https://neuromancer.sk/article/35
@contextmanager
def enable_spawn(func):
    invalidate_caches()
    source = inspect.getsource(func)
    with tempfile.NamedTemporaryFile(suffix=".py", mode="w") as f:
        f.write(source)
        f.flush()
        path = Path(f.name)
        directory = str(path.parent)
        sys.path.append(directory)
        module = import_module(str(path.stem))
        yield getattr(module, func.__name__)
        sys.path.remove(directory)

sw = ShortWeierstrassModel()
mont = MontgomeryModel()
te = TwistedEdwardsModel()

curve25519 = get_params("other", "Curve25519", "xz")
ed25519 = get_params("other", "Ed25519", "extended")
p256_jac3 = get_params("secg", "secp256r1", "jacobian-3")
p256_jac = get_params("secg", "secp256r1", "jacobian")
p256_mod = get_params("secg", "secp256r1", "modified")
p256_proj3 = get_params("secg", "secp256r1", "projective-3")

The formulas

The following formulas were collected from open-source cryptographic libraries and are stored in the tests directory of the pyecsca toolkit.

lib_formula_defs = [
        [
            "add-bc-r1rv76-jac", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp128r1"),
            AdditionEFDFormula,
        ],
        [
            "add-bc-r1rv76-mod", #ok
            ShortWeierstrassModel,
            "modified",
            ("secg", "secp128r1"),
            AdditionEFDFormula,
        ],
        [
            "dbl-bc-r1rv76-jac", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp128r1"),
            DoublingEFDFormula,
        ],
        [
            "dbl-bc-r1rv76-mod", #ok
            ShortWeierstrassModel,
            "modified",
            ("secg", "secp128r1"),
            DoublingEFDFormula,
        ],
        [
            "dbl-bc-r1rv76-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            DoublingEFDFormula,
        ],
        [
            "ladd-bc-r1rv76-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            LadderEFDFormula,
        ],
        [
            "dbl-boringssl-p224", #ok
            ShortWeierstrassModel,
            "jacobian-3",
            ("secg", "secp224r1"),
            DoublingEFDFormula,
        ],
        [
            "add-boringssl-p224", #ok
            ShortWeierstrassModel,
            "jacobian-3",
            ("secg", "secp224r1"),
            AdditionEFDFormula,
        ],
        [
            "add-libressl-v382", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp128r1"),
            AdditionEFDFormula,
        ],
        [
            "dbl-libressl-v382", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp128r1"),
            DoublingEFDFormula,
        ],
        [
            "dbl-secp256k1-v040", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp256k1"),
            DoublingEFDFormula,
        ],
        [
            "add-openssl-z256", #ok
            ShortWeierstrassModel,
            "jacobian-3",
            ("secg", "secp256r1"),
            AdditionEFDFormula,
        ],
        [
            "add-openssl-z256a", #ok
            ShortWeierstrassModel,
            "jacobian-3",
            ("secg", "secp256r1"),
            AdditionEFDFormula,
        ],
        [
            "ladd-openssl-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            LadderEFDFormula,
        ],
        [
            "ladd-hacl-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            LadderEFDFormula,
        ],
        [
            "dbl-hacl-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            DoublingEFDFormula,
        ],
        [
            "dbl-sunec-v21", #ok
            ShortWeierstrassModel,
            "projective-3",
            ("secg", "secp256r1"),
            DoublingEFDFormula,
        ],
        [
            "add-sunec-v21", #ok
            ShortWeierstrassModel,
            "projective-3",
            ("secg", "secp256r1"),
            AdditionEFDFormula,
        ],
        [
            "add-sunec-v21-ed25519", #ok
            TwistedEdwardsModel,
            "extended-1",
            ("other", "Ed25519"),
            AdditionEFDFormula,
        ],
        [
            "dbl-sunec-v21-ed25519", #ok
            TwistedEdwardsModel,
            "extended-1",
            ("other", "Ed25519"),
            DoublingEFDFormula,
        ],
        [
            "ladd-rfc7748", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            LadderEFDFormula,
        ],
        [
            "add-bearssl-v06", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp256r1"),
            AdditionEFDFormula,
        ],
        [
            "dbl-bearssl-v06", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp256r1"),
            DoublingEFDFormula,
        ],
        [
            "add-libgcrypt-v1102", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp256r1"),
            AdditionEFDFormula,
        ],
        [
            "dbl-libgcrypt-v1102", #ok
            ShortWeierstrassModel,
            "jacobian",
            ("secg", "secp256r1"),
            DoublingEFDFormula,
        ],
        [
            "ladd-go-1214", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            LadderEFDFormula,
        ],
        [
            "add-gecc-322", #ok
            ShortWeierstrassModel,
            "jacobian-3",
            ("secg", "secp256r1"),
            AdditionEFDFormula,
        ],
        [
            "dbl-gecc-321", #ok
            ShortWeierstrassModel,
            "jacobian-3",
            ("secg", "secp256r1"),
            DoublingEFDFormula,
        ],
        [
            "ladd-boringssl-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            LadderEFDFormula,
        ],
        [
            "dbl-ipp-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            DoublingEFDFormula,
        ],
        [
            "ladd-botan-x25519", #ok
            MontgomeryModel,
            "xz",
            ("other", "Curve25519"),
            LadderEFDFormula,
        ],
]

Let’s load the formulas now.

lib_formulas = {}
base_path = files(pyecsca).joinpath("../test/data/formulas/")
for formula_def in lib_formula_defs:
    meta_path = base_path / formula_def[0]
    op3_path = base_path / (formula_def[0] + ".op3")
    model = formula_def[1]()
    formula = formula_def[4](meta_path, op3_path, formula_def[0], model.coordinates[formula_def[2]]).to_code()
    lib_formulas[formula_def[0]] = formula

len(lib_formulas)

Now we can setup some code for examining the similarities between the formulas.

def compute_similarities(formulas, curve):
    table = [["One", "Other", "Similarity (output)", "Similarity (IV)"]]

    im_iv = np.zeros((len(formulas), len(formulas)))
    im_out = np.zeros((len(formulas), len(formulas)))
    for formula in formulas:
        if formula.assumptions:
            for assumption_str in formula.assumptions_str:
                lhs, rhs = assumption_str.strip().split(" == ")
                if lhs in formula.inputs:
                    print(f"Warning, formula {formula.name} has assumptions: {assumption_str}")
    for one, other in itertools.product(formulas, formulas):
        i = formulas.index(one)
        j = formulas.index(other)
        if curve is None:
            sim = formula_similarity(one, other)
        else:
            sim = formula_similarity_fuzz(one, other, curve, 100)
        im_iv[i, j] = sim["ivs"]
        im_out[i, j] = sim["output"]
        table.append([one.name, other.name, f"{sim['output']:.2}", f"{sim['ivs']:.2}"])

    return table, im_iv, im_out

def plot_similarities(formulas, im_data, name):
    fig, ax = plt.subplots()
    im = ax.imshow(im_data, vmin=0)
    cbar_ax = fig.add_axes((0.90, 0.11, 0.04, 0.76))
    cbar = fig.colorbar(im, cax=cbar_ax)
    cbar.ax.set_ylabel(f"Similarity ({name})", rotation=-90, va="bottom")

    ax.set_xticks(np.arange(len(formulas)), labels=list(map(lambda f: f.name, formulas)), rotation=90)
    ax.set_yticks(np.arange(len(formulas)), labels=list(map(lambda f: f.name, formulas)), rotation=0)

    for i, one in enumerate(formulas):
        for j, other in enumerate(formulas):
            ax.text(j, i, f"{im_data[i, j]:.2}", ha="center", va="center", color="white")
    plt.show()

def analyze_formulas(formulas, curve=None):
    table, im_iv, im_out = compute_similarities(formulas, curve)
    plot_similarities(formulas, im_iv, "IV")
    plot_similarities(formulas, im_out, "output")

Analysis

xz_ladders = [formula for formula in mont.coordinates["xz"].formulas.values() if formula.name.startswith("ladd") or formula.name.startswith("mladd")] + [
    lib_formulas["ladd-rfc7748"], lib_formulas["ladd-hacl-x25519"], lib_formulas["ladd-openssl-x25519"], lib_formulas["ladd-bc-r1rv76-x25519"],
    lib_formulas["ladd-go-1214"], lib_formulas["ladd-boringssl-x25519"], lib_formulas["ladd-botan-x25519"]]
analyze_formulas(xz_ladders)
analyze_formulas(xz_ladders, curve25519.curve)

xz_dbls = [formula for formula in mont.coordinates["xz"].formulas.values() if formula.name.startswith("dbl")] + [
    lib_formulas["dbl-bc-r1rv76-x25519"], lib_formulas["dbl-hacl-x25519"], lib_formulas["dbl-ipp-x25519"]]
analyze_formulas(xz_dbls)
analyze_formulas(xz_dbls, curve25519.curve)

jac3_adds = [formula for formula in sw.coordinates["jacobian-3"].formulas.values() if formula.name.startswith("add") or formula.name.startswith("madd")] + [
    lib_formulas["add-boringssl-p224"], lib_formulas["add-openssl-z256"], lib_formulas["add-openssl-z256a"],
    lib_formulas["add-gecc-322"]]
jac3_adds_fixed = []
for formula in jac3_adds:
    if formula.coordinate_model != sw.coordinates["jacobian-3"]:
        formula = deepcopy(formula)
        formula.coordinate_model = sw.coordinates["jacobian-3"]
    jac3_adds_fixed.append(formula)
analyze_formulas(jac3_adds_fixed)
analyze_formulas(jac3_adds_fixed, p256_jac3.curve)

jac3_dbls = [formula for formula in sw.coordinates["jacobian-3"].formulas.values() if formula.name.startswith("dbl")] + [
        lib_formulas["dbl-boringssl-p224"], lib_formulas["dbl-gecc-321"]]
jac3_dbls_fixed = []
for formula in jac3_dbls:
    if formula.coordinate_model != sw.coordinates["jacobian-3"]:
        formula = deepcopy(formula)
        formula.coordinate_model = sw.coordinates["jacobian-3"]
    jac3_dbls_fixed.append(formula)
analyze_formulas(jac3_dbls_fixed)
analyze_formulas(jac3_dbls_fixed, p256_jac3.curve)

mod_adds = [formula for formula in sw.coordinates["modified"].formulas.values() if formula.name.startswith("add") or formula.name.startswith("madd")] + [
    lib_formulas["add-bc-r1rv76-mod"]]
analyze_formulas(mod_adds)
analyze_formulas(mod_adds, p256_mod.curve)

mod_dbls = [formula for formula in sw.coordinates["modified"].formulas.values() if formula.name.startswith("dbl")] + [
    lib_formulas["dbl-bc-r1rv76-mod"]]
analyze_formulas(mod_dbls)
analyze_formulas(mod_dbls, p256_mod.curve)

jac_adds = [formula for formula in sw.coordinates["jacobian"].formulas.values() if formula.name.startswith("add")] + [
    lib_formulas["add-bc-r1rv76-jac"], lib_formulas["add-libressl-v382"], lib_formulas["add-bearssl-v06"], lib_formulas["add-libgcrypt-v1102"]]
analyze_formulas(jac_adds)
analyze_formulas(jac_adds, p256_jac.curve)

jac_dbls = [formula for formula in sw.coordinates["jacobian"].formulas.values() if formula.name.startswith("dbl")] + [
    lib_formulas["dbl-secp256k1-v040"], lib_formulas["dbl-libressl-v382"], lib_formulas["dbl-bearssl-v06"], lib_formulas["dbl-libgcrypt-v1102"]]
analyze_formulas(jac_dbls)
analyze_formulas(jac_dbls, p256_jac.curve)

proj3_adds = [formula for formula in sw.coordinates["projective-3"].formulas.values() if formula.name.startswith("add") or formula.name.startswith("madd")] + [
    lib_formulas["add-sunec-v21"]]
analyze_formulas(proj3_adds)
analyze_formulas(proj3_adds, p256_proj3.curve)

proj3_dbls = [formula for formula in sw.coordinates["projective-3"].formulas.values() if formula.name.startswith("dbl")] + [
    lib_formulas["dbl-sunec-v21"]]
analyze_formulas(proj3_dbls)
analyze_formulas(proj3_dbls, p256_proj3.curve)

ext_adds = [formula for formula in te.coordinates["extended-1"].formulas.values() if formula.name.startswith("add") or formula.name.startswith("madd")] + [
    lib_formulas["add-sunec-v21-ed25519"]]
analyze_formulas(ext_adds)

ext_dbls = [formula for formula in te.coordinates["extended-1"].formulas.values() if formula.name.startswith("dbl")] + [
    lib_formulas["dbl-sunec-v21-ed25519"]]
analyze_formulas(ext_dbls)

Expand

We can also expand the set of formulas by applying various transformations (like swapping the order of commutative operations).

Note that this computation is parallelized and you should set an appropriate number of workers.

def expand_out(formulas, name):
    from pyecsca.ec.formula.expand import expand_formula_set
    import pickle

    expanded = expand_formula_set(formulas)
    with open(name, "wb") as f:
        pickle.dump(expanded, f)
    return name, len(expanded)

# 22 is really a maximum usable number here, as there are only 22 "jobs".
context = multiprocessing.get_context("spawn")
with ProcessPoolExecutor(max_workers=22, mp_context=context) as pool, enable_spawn(expand_out) as expand_func:
    futures = []
    args = []
    for coord_name, coords in tqdm(sw.coordinates.items(), desc="Submitting"):
        adds = [formula for formula in coords.formulas.values() if formula.name.startswith("add")]
        lib_adds = [formula for formula in lib_formulas.values() if formula.coordinate_model == coords and formula.name.startswith("add")]
        dbls = [formula for formula in coords.formulas.values() if formula.name.startswith("dbl")]
        lib_dbls = [formula for formula in lib_formulas.values() if formula.coordinate_model == coords and formula.name.startswith("dbl")]
        futures.append(pool.submit(expand_func, adds + lib_adds, f"sw_{coord_name}_adds.pickle"))
        args.append(f"sw_{coord_name}_adds.pickle")
        futures.append(pool.submit(expand_func, dbls + lib_dbls, f"sw_{coord_name}_dbls.pickle"))
        args.append(f"sw_{coord_name}_dbls.pickle")
    for future in tqdm(as_completed(futures), total=len(futures), smoothing=0, desc="Computing"):
        j = futures.index(future)
        arg = args[j]
        error = future.exception()
        if error:
            print(arg, error)
        else:
            res = future.result()
            print(*res)