ZVP-based reverse-engineering

This notebook showcases the ZVP-based reverse-engineering technique for addition formulas.

[ ]:
import io
import numpy as np
import pandas as pd
import holoviews as hv
import xarray as xr
import random
import tabulate
import pickle
import multiprocessing
import inspect
import tempfile
import sys
import re
from collections import Counter
from matplotlib import pyplot as plt
from sympy import FF, ZZ, sympify, symbols, Poly
from contextlib import contextmanager
from importlib import import_module, invalidate_caches
from pathlib import Path
from functools import partial
from itertools import product
from IPython.display import HTML, display
from tqdm.notebook import tqdm, trange
from anytree import RenderTree, LevelOrderIter


from pyecsca.ec.model import ShortWeierstrassModel
from pyecsca.ec.coordinates import AffineCoordinateModel
from pyecsca.ec.curve import EllipticCurve
from pyecsca.ec.params import DomainParameters, load_params_ecgen
from pyecsca.ec.formula import FormulaAction, AdditionFormula, DoublingFormula
from pyecsca.ec.point import Point
from pyecsca.ec.mod import Mod, gcd, SymbolicMod
from pyecsca.ec.context import DefaultContext, local
from pyecsca.ec.mult import LTRMultiplier, AccumulationOrder
from pyecsca.ec.error import NonInvertibleError, UnsatisfiedAssumptionError
from pyecsca.sca.re.tree import Map, Tree
from pyecsca.sca.re.rpa import MultipleContext
from pyecsca.sca.re.zvp import zvp_points, compute_factor_set, unroll_formula,  addition_chain, eliminate_y
from pyecsca.misc.cfg import getconfig
from pyecsca.misc.utils import TaskExecutor

from eval import (eval_tree_symmetric1, eval_tree_asymmetric1, eval_tree_binomial1,
                    success_rate_symmetric, success_rate_asymmetric, success_rate_binomial,
                    query_rate_symmetric, query_rate_asymmetric, query_rate_binomial,
                    amount_rate_symmetric, amount_rate_asymmetric, amount_rate_binomial,
                    precise_rate_symmetric, precise_rate_asymmetric, precise_rate_binomial,
                    success_rate_vs_majority_symmetric, success_rate_vs_majority_asymmetric,
                    success_rate_vs_query_rate_symmetric, load, store)


# 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)

spawn_context = multiprocessing.get_context("spawn")

%matplotlib ipympl
hv.extension("bokeh")
[ ]:
model = ShortWeierstrassModel()
coordsaff = AffineCoordinateModel(model)

Exploration

First lets explore the behavior of addition formulas. The following two cells pick a coordinate model along with some formulas and symbolically unroll a scalar multiplication (assuming a simple LTR multiplier).

[ ]:
which = "jacobian"
coords = model.coordinates[which]
[ ]:
getconfig().ec.mod_implementation = "symbolic"
x, y, z = symbols("x y z")

# A 64-bit prime order curve for testing things out
p = 0xc50de883f0e7b167
field = FF(p)
a = SymbolicMod(Poly(0x4833d7aa73fa6694, x, y, z, domain=field), p)
b = SymbolicMod(Poly(0xa6c44a61c5323f6a, x, y, z, domain=field), p)
gx = SymbolicMod(Poly(0x5fd1f7d38d4f2333, x, y, z, domain=field), p)
gy = SymbolicMod(Poly(0x21f43957d7e20ceb, x, y, z, domain=field), p)
n = 0xc50de885003b80eb
h = 1

infty = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))
g = Point(coords, X=gx, Y=gy, Z=Mod(1, p))

curve = EllipticCurve(model, coords, p, infty, dict(a=a,b=b))
params = DomainParameters(curve, g, n, h)


add = coords.formulas["add-2007-bl"]
dbl = coords.formulas["dbl-2007-bl"]
mult = LTRMultiplier(add, dbl, None, False, AccumulationOrder.PeqRP, True, True)


point = Point(coords,
              X=SymbolicMod(Poly(x, x, y, z, domain=field), params.curve.prime),
              Y=SymbolicMod(Poly(y, x, y, z, domain=field), params.curve.prime),
              Z=SymbolicMod(Poly(z, x, y, z, domain=field), params.curve.prime))
mult.init(params, point)
res = mult.multiply(5)

x_poly = Poly(res.X.x, domain=field)
y_poly = Poly(res.Y.x, domain=field)
z_poly = Poly(res.Z.x, domain=field)
display(x_poly, y_poly, z_poly)

The result is a Point with coordinates that are polynomials in the input coordinates and curve parameters. We now switch back to concrete representation.

[ ]:
cfg = getconfig()
cfg.ec.mod_implementation = "gmp"

Reverse-engineering

Now, lets look at using the ZVP attack for reverse-engineering. First pick 10 curves per group, some random some with \(a \in \{0, -1, -3 \}\). The curves are otherwise not special in any way and just serve to randomize the process, as the existence of ZVP points for a given intermediate value polynomial depends on the curve.

[ ]:
curves = list(map(lambda spec: load_params_ecgen(io.BytesIO(spec.encode()), "affine"), [
    # Random
    """[{"field":{"p":"0xddf438409fc35161"},"a":"0x94d919b72f7dc6d8","b":"0x9f39032abb23f62a","order":"0xddf4383ffa8e6de7","subgroups":[{"x":"0xd5673b3fe176fc6b","y":"0x2d5b0a5bb2141317","order":"0xddf4383ffa8e6de7","cofactor":"0x1","points":[{"x":"0xd5673b3fe176fc6b","y":"0x2d5b0a5bb2141317","order":"0xddf4383ffa8e6de7"}]}]}]""",
    """[{"field":{"p":"0xa42c1467a1ed04f3"},"a":"0x55d07340a4572f2d","b":"0x0a938c37dfb0b6d5","order":"0xa42c14689284d3a7","subgroups":[{"x":"0x8633981c83ed43a2","y":"0x7b5374e9d7997199","order":"0xa42c14689284d3a7","cofactor":"0x1","points":[{"x":"0x8633981c83ed43a2","y":"0x7b5374e9d7997199","order":"0xa42c14689284d3a7"}]}]}]""",
    """[{"field":{"p":"0xea0d9cead19016ab"},"a":"0xcbbfe501c4ef6d92","b":"0x5762de777a6d9178","order":"0xea0d9cea8cd2c857","subgroups":[{"x":"0xe7daa3e061c3111b","y":"0x56ee59a6845c5e93","order":"0xea0d9cea8cd2c857","cofactor":"0x1","points":[{"x":"0xe7daa3e061c3111b","y":"0x56ee59a6845c5e93","order":"0xea0d9cea8cd2c857"}]}]}]""",
    """[{"field":{"p":"0x9c7e7216decb71a7"},"a":"0x324ef48887401a87","b":"0x3ce6f35a00280102","order":"0x9c7e72175ebfe709","subgroups":[{"x":"0x34683229b405418d","y":"0x308c923cae004514","order":"0x9c7e72175ebfe709","cofactor":"0x1","points":[{"x":"0x34683229b405418d","y":"0x308c923cae004514","order":"0x9c7e72175ebfe709"}]}]}]""",
    """[{"field":{"p":"0xeb5779f0bbf1ef5b"},"a":"0x2419e8bbc7b5f8f2","b":"0xe74e5d3064a4f2e3","order":"0xeb5779f21320c2e9","subgroups":[{"x":"0x3b6c269560abeb00","y":"0x29d157628e75e1c0","order":"0xeb5779f21320c2e9","cofactor":"0x1","points":[{"x":"0x3b6c269560abeb00","y":"0x29d157628e75e1c0","order":"0xeb5779f21320c2e9"}]}]}]""",
    """[{"field":{"p":"0x97b6ea097868b95d"},"a":"0x550a41d65e4bcd13","b":"0x47c5e527113b261c","order":"0x97b6ea094947a76b","subgroups":[{"x":"0x1e669fe19c865bd9","y":"0x05a6bb891920440f","order":"0x97b6ea094947a76b","cofactor":"0x1","points":[{"x":"0x1e669fe19c865bd9","y":"0x05a6bb891920440f","order":"0x97b6ea094947a76b"}]}]}]""",
    """[{"field":{"p":"0xa00629e6522032f7"},"a":"0x896f04a7ae302922","b":"0x6bc03365b1f1cb50","order":"0xa00629e5c03cf913","subgroups":[{"x":"0x14b7b48954936d4e","y":"0x670dc776273bf899","order":"0xa00629e5c03cf913","cofactor":"0x1","points":[{"x":"0x14b7b48954936d4e","y":"0x670dc776273bf899","order":"0xa00629e5c03cf913"}]}]}]""",
    """[{"field":{"p":"0xd47ec1d03a62686d"},"a":"0xd00a3ee0f5c86b02","b":"0x457a5b6c47db38d8","order":"0xd47ec1d107db7d6f","subgroups":[{"x":"0x41ebc3b763f3cd1b","y":"0x3d6925f214620e0c","order":"0xd47ec1d107db7d6f","cofactor":"0x1","points":[{"x":"0x41ebc3b763f3cd1b","y":"0x3d6925f214620e0c","order":"0xd47ec1d107db7d6f"}]}]}]""",
    """[{"field":{"p":"0xb1c9115c6f40d755"},"a":"0x79d3ceefafc44ce9","b":"0x8316af84264df42b","order":"0xb1c9115d17f84a45","subgroups":[{"x":"0x8b0a274089b53fe5","y":"0x3508d33c4beba5ad","order":"0xb1c9115d17f84a45","cofactor":"0x1","points":[{"x":"0x8b0a274089b53fe5","y":"0x3508d33c4beba5ad","order":"0xb1c9115d17f84a45"}]}]}]""",
    """[{"field":{"p":"0x8f738fda18cd5dff"},"a":"0x4747f2f9b8628cbf","b":"0x586cdb9378a1389f","order":"0x8f738fd8fc7ebed3","subgroups":[{"x":"0x7ad306c73b64c1b5","y":"0x69e3ca555190da4b","order":"0x8f738fd8fc7ebed3","cofactor":"0x1","points":[{"x":"0x7ad306c73b64c1b5","y":"0x69e3ca555190da4b","order":"0x8f738fd8fc7ebed3"}]}]}]""",
    # a = -1
    """[{"field":{"p":"0xcfef393139c3007f"},"a":"0xcfef393139c3007e","b":"0x950312812acb155f","order":"0xcfef39320179387b","subgroups":[{"x":"0xae2d2f58ca5b5cf7","y":"0xc3a4bf3a1dc10005","order":"0xcfef39320179387b","cofactor":"0x1","points":[{"x":"0xae2d2f58ca5b5cf7","y":"0xc3a4bf3a1dc10005","order":"0xcfef39320179387b"}]}]}]""",
    """[{"field":{"p":"0xb0461c0e4946cbd5"},"a":"0xb0461c0e4946cbd4","b":"0x082c3722016def51","order":"0xb0461c0e07e3e1bf","subgroups":[{"x":"0x5142200263be1fe3","y":"0x14984b7551ed21a9","order":"0xb0461c0e07e3e1bf","cofactor":"0x1","points":[{"x":"0x5142200263be1fe3","y":"0x14984b7551ed21a9","order":"0xb0461c0e07e3e1bf"}]}]}]""",
    """[{"field":{"p":"0xeff607c2dc4f278b"},"a":"0xeff607c2dc4f278a","b":"0x26fd03674f5092d2","order":"0xeff607c30ab8c50d","subgroups":[{"x":"0x004d4a5a9bb849fe","y":"0x80eb7ef89110c149","order":"0xeff607c30ab8c50d","cofactor":"0x1","points":[{"x":"0x004d4a5a9bb849fe","y":"0x80eb7ef89110c149","order":"0xeff607c30ab8c50d"}]}]}]""",
    """[{"field":{"p":"0xedc14fda51686379"},"a":"0xedc14fda51686378","b":"0xb3973a86901e3364","order":"0xedc14fda0cdbc199","subgroups":[{"x":"0xc76f0776feb59336","y":"0x625adaf0fb44ab9f","order":"0xedc14fda0cdbc199","cofactor":"0x1","points":[{"x":"0xc76f0776feb59336","y":"0x625adaf0fb44ab9f","order":"0xedc14fda0cdbc199"}]}]}]""",
    """[{"field":{"p":"0xfc6ee07288f1b78f"},"a":"0xfc6ee07288f1b78e","b":"0xe18821a83ca2ca30","order":"0xfc6ee0713e07f37f","subgroups":[{"x":"0x339d01a4b0db428e","y":"0x68100d42e5ffd979","order":"0xfc6ee0713e07f37f","cofactor":"0x1","points":[{"x":"0x339d01a4b0db428e","y":"0x68100d42e5ffd979","order":"0xfc6ee0713e07f37f"}]}]}]""",
    """[{"field":{"p":"0xa03c03a0072f69b3"},"a":"0xa03c03a0072f69b2","b":"0x3208ecebb633b82c","order":"0xa03c039ff31e37a7","subgroups":[{"x":"0x8134208d53e6f6c0","y":"0x6245db54032630a6","order":"0xa03c039ff31e37a7","cofactor":"0x1","points":[{"x":"0x8134208d53e6f6c0","y":"0x6245db54032630a6","order":"0xa03c039ff31e37a7"}]}]}]""",
    """[{"field":{"p":"0xbc8c6e7ce26746d9"},"a":"0xbc8c6e7ce26746d8","b":"0xb7e2b4fb2d769c4e","order":"0xbc8c6e7ba032dda7","subgroups":[{"x":"0x8e3c9cd771e7ffd8","y":"0x4dd02403ca890c5a","order":"0xbc8c6e7ba032dda7","cofactor":"0x1","points":[{"x":"0x8e3c9cd771e7ffd8","y":"0x4dd02403ca890c5a","order":"0xbc8c6e7ba032dda7"}]}]}]""",
    """[{"field":{"p":"0x9ccda4c062b95787"},"a":"0x9ccda4c062b95786","b":"0x31fcbb278951e3b9","order":"0x9ccda4bfae73e4f5","subgroups":[{"x":"0x303ac583c81644e3","y":"0x76713f6f470e94a0","order":"0x9ccda4bfae73e4f5","cofactor":"0x1","points":[{"x":"0x303ac583c81644e3","y":"0x76713f6f470e94a0","order":"0x9ccda4bfae73e4f5"}]}]}]""",
    """[{"field":{"p":"0xa339e3745518416f"},"a":"0xa339e3745518416e","b":"0x52d39a67f2401673","order":"0xa339e3743950389b","subgroups":[{"x":"0x6b8986f706afac58","y":"0x5c901b1afa0b64da","order":"0xa339e3743950389b","cofactor":"0x1","points":[{"x":"0x6b8986f706afac58","y":"0x5c901b1afa0b64da","order":"0xa339e3743950389b"}]}]}]""",
    """[{"field":{"p":"0x8b7d2baee4e47311"},"a":"0x8b7d2baee4e47310","b":"0x21ab23afb5a9e2ca","order":"0x8b7d2baf201f2bdd","subgroups":[{"x":"0x797c1dec0d73ec64","y":"0x28f90926ea9c6b33","order":"0x8b7d2baf201f2bdd","cofactor":"0x1","points":[{"x":"0x797c1dec0d73ec64","y":"0x28f90926ea9c6b33","order":"0x8b7d2baf201f2bdd"}]}]}]""",
    # a = -3
    """[{"field":{"p":"0x8d79ca36cee026a7"},"a":"0x8d79ca36cee026a4","b":"0x0478c1f80ce2c9c6","order":"0x8d79ca35a428c76f","subgroups":[{"x":"0x2e94a3e38f8b345e","y":"0x83e6c6f0cb8f69c4","order":"0x8d79ca35a428c76f","cofactor":"0x1","points":[{"x":"0x2e94a3e38f8b345e","y":"0x83e6c6f0cb8f69c4","order":"0x8d79ca35a428c76f"}]}]}]""",
    """[{"field":{"p":"0x48e1a125250323a7"},"a":"0x48e1a125250323a4","b":"0x02a4d99f41d23210","order":"0x48e1a124f895db6d","subgroups":[{"x":"0x409e15d65fcae55a","y":"0x207e142056d62d07","order":"0x48e1a124f895db6d","cofactor":"0x1","points":[{"x":"0x409e15d65fcae55a","y":"0x207e142056d62d07","order":"0x48e1a124f895db6d"}]}]}]""",
    """[{"field":{"p":"0xcb5aa8a7a10aa06b"},"a":"0xcb5aa8a7a10aa068","b":"0x31fe9c57c570174f","order":"0xcb5aa8a6cf812191","subgroups":[{"x":"0x84c75d46fc687ff1","y":"0x7424362ac73df187","order":"0xcb5aa8a6cf812191","cofactor":"0x1","points":[{"x":"0x84c75d46fc687ff1","y":"0x7424362ac73df187","order":"0xcb5aa8a6cf812191"}]}]}]""",
    """[{"field":{"p":"0xba965ca9c8aa0a1b"},"a":"0xba965ca9c8aa0a18","b":"0x676535a1eaf5c605","order":"0xba965caae5741b6f","subgroups":[{"x":"0x313d58c47b8ed95f","y":"0x991ba98cbbb0fe9f","order":"0xba965caae5741b6f","cofactor":"0x1","points":[{"x":"0x313d58c47b8ed95f","y":"0x991ba98cbbb0fe9f","order":"0xba965caae5741b6f"}]}]}]""",
    """[{"field":{"p":"0xbfb7747454a17d15"},"a":"0xbfb7747454a17d12","b":"0x611b69881db4ce69","order":"0xbfb7747547fd57d3","subgroups":[{"x":"0x3385044d698640fc","y":"0x50cee623251b559e","order":"0xbfb7747547fd57d3","cofactor":"0x1","points":[{"x":"0x3385044d698640fc","y":"0x50cee623251b559e","order":"0xbfb7747547fd57d3"}]}]}]""",
    """[{"field":{"p":"0x99235f1ed44b3959"},"a":"0x99235f1ed44b3956","b":"0x5d8dda19dbe804d4","order":"0x99235f1d975f376d","subgroups":[{"x":"0x4fed262974c1d800","y":"0x27590c454edd59ca","order":"0x99235f1d975f376d","cofactor":"0x1","points":[{"x":"0x4fed262974c1d800","y":"0x27590c454edd59ca","order":"0x99235f1d975f376d"}]}]}]""",
    """[{"field":{"p":"0xa7ff74a0dc8c161d"},"a":"0xa7ff74a0dc8c161a","b":"0x583b968bb611b284","order":"0xa7ff74a06811ee75","subgroups":[{"x":"0x5f5c76454edf12e7","y":"0x4c73cbfc44f41508","order":"0xa7ff74a06811ee75","cofactor":"0x1","points":[{"x":"0x5f5c76454edf12e7","y":"0x4c73cbfc44f41508","order":"0xa7ff74a06811ee75"}]}]}]""",
    """[{"field":{"p":"0xb52c62ca8703a063"},"a":"0xb52c62ca8703a060","b":"0x0baec43a07b54c21","order":"0xb52c62c963037121","subgroups":[{"x":"0x6fe4a521a29bc1ab","y":"0x3fca7180021f8f0f","order":"0xb52c62c963037121","cofactor":"0x1","points":[{"x":"0x6fe4a521a29bc1ab","y":"0x3fca7180021f8f0f","order":"0xb52c62c963037121"}]}]}]""",
    """[{"field":{"p":"0xb8921f25b6ce5267"},"a":"0xb8921f25b6ce5264","b":"0xa575c9f97563265d","order":"0xb8921f2592b6b39f","subgroups":[{"x":"0x7eb120fada47765c","y":"0x64ef4e51d4159304","order":"0xb8921f2592b6b39f","cofactor":"0x1","points":[{"x":"0x7eb120fada47765c","y":"0x64ef4e51d4159304","order":"0xb8921f2592b6b39f"}]}]}]""",
    """[{"field":{"p":"0xc591b8c4df0afc19"},"a":"0xc591b8c4df0afc16","b":"0x0a1eb46a6e647f0a","order":"0xc591b8c3eb07239f","subgroups":[{"x":"0x1963bfb862cb0bf3","y":"0x30da8bb7fa77277d","order":"0xc591b8c3eb07239f","cofactor":"0x1","points":[{"x":"0x1963bfb862cb0bf3","y":"0x30da8bb7fa77277d","order":"0xc591b8c3eb07239f"}]}]}]""",
    # a = 0
    """[{"field":{"p":"0xceaf446a53f14bc1"},"a":"0x0000000000000000","b":"0x326539376260f173","order":"0xceaf446aae275419","subgroups":[{"x":"0x98fe44948c3f8678","y":"0x3d440ee959a912d7","order":"0xceaf446aae275419","cofactor":"0x1","points":[{"x":"0x98fe44948c3f8678","y":"0x3d440ee959a912d7","order":"0xceaf446aae275419"}]}]}]""",
    """[{"field":{"p":"0xb3c2beca75d66de3"},"a":"0x0000000000000000","b":"0x46069225826b51aa","order":"0xb3c2bec95881b695","subgroups":[{"x":"0x81500c226efa0d5a","y":"0x674e09d296452eee","order":"0xb3c2bec95881b695","cofactor":"0x1","points":[{"x":"0x81500c226efa0d5a","y":"0x674e09d296452eee","order":"0xb3c2bec95881b695"}]}]}]""",
    """[{"field":{"p":"0xd6097c1ce207aae7"},"a":"0x0000000000000000","b":"0x7adaab54e7dfd564","order":"0xd6097c1b407eb413","subgroups":[{"x":"0x151da8fb1f83201e","y":"0x8bfeb90ec1177a91","order":"0xd6097c1b407eb413","cofactor":"0x1","points":[{"x":"0x151da8fb1f83201e","y":"0x8bfeb90ec1177a91","order":"0xd6097c1b407eb413"}]}]}]""",
    """[{"field":{"p":"0x97a3e2d617a2309d"},"a":"0x0000000000000000","b":"0x7f311cba46652247","order":"0x97a3e2d712ffd715","subgroups":[{"x":"0x46d725812af15870","y":"0x727f88365dbd0e80","order":"0x97a3e2d712ffd715","cofactor":"0x1","points":[{"x":"0x46d725812af15870","y":"0x727f88365dbd0e80","order":"0x97a3e2d712ffd715"}]}]}]""",
    """[{"field":{"p":"0xd8d7f39f545a5da7"},"a":"0x0000000000000000","b":"0x09097ddcd4be8d55","order":"0xd8d7f3a0c4115b4b","subgroups":[{"x":"0xc5c771a9827a3251","y":"0x64bf52041ac05b23","order":"0xd8d7f3a0c4115b4b","cofactor":"0x1","points":[{"x":"0xc5c771a9827a3251","y":"0x64bf52041ac05b23","order":"0xd8d7f3a0c4115b4b"}]}]}]""",
    """[{"field":{"p":"0x847de883ab4fbf4d"},"a":"0x0000000000000000","b":"0x09866366b3b45c2d","order":"0x847de8837e6d4477","subgroups":[{"x":"0x62fd7b4bc7c9acb4","y":"0x2d0942774607106b","order":"0x847de8837e6d4477","cofactor":"0x1","points":[{"x":"0x62fd7b4bc7c9acb4","y":"0x2d0942774607106b","order":"0x847de8837e6d4477"}]}]}]""",
    """[{"field":{"p":"0xf0d617c3c47b7c77"},"a":"0x0000000000000000","b":"0xd856b3dcb95764a2","order":"0xf0d617c5512cec85","subgroups":[{"x":"0xeaf9b352a3daac45","y":"0x4e4e557f9fc3febc","order":"0xf0d617c5512cec85","cofactor":"0x1","points":[{"x":"0xeaf9b352a3daac45","y":"0x4e4e557f9fc3febc","order":"0xf0d617c5512cec85"}]}]}]""",
    """[{"field":{"p":"0xce920b656c80b373"},"a":"0x0000000000000000","b":"0xb4a07dfae71ddc62","order":"0xce920b65eee38015","subgroups":[{"x":"0x7895c02b3c5205b5","y":"0x2926be6446b98d62","order":"0xce920b65eee38015","cofactor":"0x1","points":[{"x":"0x7895c02b3c5205b5","y":"0x2926be6446b98d62","order":"0xce920b65eee38015"}]}]}]""",
    """[{"field":{"p":"0xb6d4e25f76cc9df7"},"a":"0x0000000000000000","b":"0x18e95a2283692623","order":"0xb6d4e25ea270ed03","subgroups":[{"x":"0x2da7a97d5d899bc5","y":"0x17d27fd34562e3d9","order":"0xb6d4e25ea270ed03","cofactor":"0x1","points":[{"x":"0x2da7a97d5d899bc5","y":"0x17d27fd34562e3d9","order":"0xb6d4e25ea270ed03"}]}]}]""",
    """[{"field":{"p":"0xe7cd1979ebed69ed"},"a":"0x0000000000000000","b":"0x278e92b83191a7da","order":"0xe7cd197966893365","subgroups":[{"x":"0xc4de44402da5b9a6","y":"0x2b45e7f32e3701ba","order":"0xe7cd197966893365","cofactor":"0x1","points":[{"x":"0xc4de44402da5b9a6","y":"0x2b45e7f32e3701ba","order":"0xe7cd197966893365"}]}]}]"""
    # b = 0 (causes more issues than gain) but the (w12-0) coords actually require it...
    #"""[{"field":{"p":"0x9d9119957f02fe3f"},"a":"0x0106903196d88df9","b":"0x0000000000000000","order":"0x9d9119957f02fe40","subgroups":[{"x":"0x191a36b9cd81de96","y":"0x10f2c6bded391aa9","order":"0x9d9119957f02fe40","cofactor":"0x1","points":[{"x":"0x0000000000000000","y":"0x0000000000000000","order":"0x2"},{"x":"0x95913fae9065da0f","y":"0x5eeddeee7152d6fb","order":"0x276446655fc0bf9"}]}]}]"""
]))

for i, params in enumerate(curves):
    curve = params.curve
    if curve.parameters["a"] == -3:
        params.name = "a=-3"
    elif curve.parameters["a"] == -1:
        params.name = "a=-1"
    elif curve.parameters["a"] == 0:
        params.name = "a=0"
    else:
        params.name = "random"
    params.name += f"[{i}]"

Computing addition chains

First lets fix some scalars, go over the curves and compute the addition chain to obtain information about which multiples of the input point will go into the formulas. The i-th scalar will be used with the i-th curve as defined above. There are 10 unique scalars, so each curve group will share those.

[ ]:
scalars = [0b1111110, 0b1011101, 0b110110, 0b100100, 0b1000110, 0b1001101, 0b101001, 0b1100100, 0b1010110, 0b101010] * 4

chains = []
scalar_map = {}
chain_map = {}
ops = set()
for scalar, params in zip(scalars, curves):
    chain = addition_chain(scalar, params, LTRMultiplier, lambda add,dbl: LTRMultiplier(add, dbl, None, False, AccumulationOrder.PeqRP, True, True))
    chains.append(chain)
    scalar_map[params] = scalar
    chain_map[params] = chain
    ops.update(chain)
print(sorted(list(ops)))

Loading the formulas

Now lets load the formulas, either just those from the EFD or also load the expanded library formulas if they are available. See the formulas notebook.

[ ]:
load_expanded = False

formula_classes = [AdditionFormula, DoublingFormula]
formula_groups = {}
for coord_name, coords in tqdm(model.coordinates.items(), desc=f"Loading {'expanded' if load_expanded else 'EFD'} formulas"):
    groups = []
    for formula_class in formula_classes:
        expanded_path = Path(f"sw_{coord_name}_{formula_class.shortname}s.pickle")
        if load_expanded:
            if not expanded_path.exists():
                raise ValueError(f"Expanded formulas do not exist {expanded_path}.")
            with expanded_path.open("rb") as f:
                expanded = pickle.load(f)
            formula_group = list(expanded)
        else:
            formula_group = list(filter(lambda formula: isinstance(formula, formula_class) and (formula.name.startswith("add") or formula.name.startswith("dbl")), coords.formulas.values()))
        groups.append(formula_group)
    formula_groups[coords] = groups

print(f"Loaded {sum(sum(len(group) for group in pair) for pair in formula_groups.values())} formulas.")

Computing the factor sets

Now, compute the factor sets of the formulas in parallel. There are two options available. xonly_fsets builds the factor sets “x”-only by eliminating y-coords using the curve equation. filter_nonhomo specifies whether to filter out non-homogenous polynomials.

[ ]:
xonly_fsets = False
filter_nonhomo = False

def compute_fsets(formula_group, formula_class, fh, xo):
    from pyecsca.sca.re.zvp import compute_factor_set
    from pyecsca.ec.formula import DoublingFormula
    fsets = []
    for formula in formula_group:
        fset = compute_factor_set(formula, filter_nonhomo=fh, xonly=xo)

        # Fix the factor set polynomials for the case of doubling.
        # TODO: Investigate how this plays with xonly an filter_nonhomo arguments.
        if formula_class == DoublingFormula:
            new_fset = set()
            for poly in fset:
                pl = poly.copy()
                for symbol in poly.free_symbols:
                    original = str(symbol)
                    if original.endswith("1"):
                        new = original.replace("1", "2")
                        pl = pl.subs(original, new)
                new_fset.add(pl)
            fset = new_fset
        fsets.append(fset)
    return fsets

factor_sets = {}
with TaskExecutor(max_workers=22, mp_context=spawn_context) as pool, enable_spawn(compute_fsets) as target:
    for coord_name, coords in model.coordinates.items():
        for formula_group, formula_class in zip(formula_groups[coords], formula_classes):
            pool.submit_task((coords, formula_group, formula_class),
                             target, formula_group, formula_class, filter_nonhomo, xonly_fsets)

    for (coords, formula_group, formula_class), future in tqdm(pool.as_completed(), desc="Computing factor sets", total=len(pool.tasks)):
        if error := future.exception():
            print(coords, formula_class.shortname, error)
            raise error
        else:
            fsets = future.result()
            print(f"Got {coords.name} {formula_class.shortname}s {len(fsets)}.")
            for formula, fset in zip(formula_group, fsets):
                factor_sets[formula] = fset

You can now store the (or load the previously computed) factor sets.

[ ]:
with open("factor_sets.pickle", "wb") as f:
    pickle.dump(factor_sets, f)
    print(f"Stored {len(factor_sets)} factor sets.")
[ ]:
with open("factor_sets.pickle", "rb") as f:
    factor_sets = pickle.load(f)
    print(f"Loaded {len(factor_sets)} factor sets.")

Accumulating polynomials

We now accumulate all of the polynomials for the adds and the dbls. We do so for all the compatible curves from our curves list. We will be looking for ZVP points on all of these curves, to make sure at least some lead to solutions.

[ ]:
polynomials = {}
for coord_name, coords in tqdm(model.coordinates.items(), desc="Accumulating"):
    coord_adds = 0
    coord_dbls = 0
    for chain, affine_params in zip(chains, curves):
        try:
            params = affine_params.to_coords(coords)
        except UnsatisfiedAssumptionError:
            continue
        add_polynomials = set()
        dbl_polynomials = set()
        for formula_group, formula_class in zip(formula_groups[coords], formula_classes):
            for formula in formula_group:
                if formula_class == AdditionFormula:
                    add_polynomials.update(factor_sets.get(formula, []))
                else:
                    dbl_polynomials.update(factor_sets.get(formula, []))
        polynomials[params] = {
            "add": add_polynomials,
            "dbl": dbl_polynomials
        }
        coord_adds += len(add_polynomials)
        coord_dbls += len(dbl_polynomials)
    print(f"Got {coord_adds} add polys and {coord_dbls} dbl polys for {coord_name}.")

Computing ZVP points

Now, lets compute the sets of ZVP points, going over all of the polynomials. Also filter the points such that for each “polynomial, curve category, k” we have only one point, as more are unneccessary.

[ ]:
# bound is the maximal dlog in the hard case of the DCP to be solved
bound = 100
# Note that if you do not have the "pari" extra dependency installed ("cysignals", "cypari2") this bound
# will have to be limited very low and the memory usage will be significant.

all_points = set()
all_points_filtered = {}
dk = set()
with TaskExecutor(max_workers=20, mp_context=spawn_context) as pool:
    for coord_name, coords in tqdm(model.coordinates.items(), desc="Submitting"):
        for chain, affine_params in zip(chains, curves):
            try:
                params = affine_params.to_coords(coords)
            except UnsatisfiedAssumptionError:
                continue
            unique = set()
            for op, ks in chain:
                if len(ks) == 1:
                    k = ks[0]
                else:
                    # zvp_points assumes (P, [k]P)
                    ks_mod = list(map(lambda v: Mod(v, params.order), ks))
                    k = int(ks_mod[1] / ks_mod[0])
                polys = polynomials[params][op]
                for poly in polys:
                    only_1 = all((not str(gen).endswith("2")) for gen in poly.gens)
                    only_2 = all((not str(gen).endswith("1")) for gen in poly.gens)
                    # This is the hard case where a dlog needs to be substituted, so bound it.
                    if not (only_1 or only_2) and k > bound:
                        continue
                    unique.add((poly, k))
            for poly, k in unique:
                pool.submit_task((poly, affine_params, k),
                                 zvp_points, poly, params.curve, k, params.order)
    for (poly, affine_params, k), future in tqdm(pool.as_completed(), desc="Computing", total=len(pool.tasks), smoothing=0):
        params_name_match = re.match("(.+)\[([0-9]+)\]", affine_params.name)
        params_category = params_name_match.group(1)
        if error := future.exception():
            print(error)
        elif result := future.result():
            for point in result:
                all_points.add((point, affine_params))
            all_points_filtered[(poly, params_category, k)] = (result, affine_params)

print(f"Got {len(all_points)} points.")
print(f"Got {len(all_points_filtered)} filtered points")#, but just {len(set(all_points_filtered.values()))} unique.")

You can now store the (or load the previously computed) point sets.

[ ]:
with open("all_points.pickle", "wb") as f:
    pickle.dump((all_points, all_points_filtered), f)
[ ]:
with open("all_points.pickle", "rb") as f:
    all_points, all_points_filtered = pickle.load(f)

print(f"Got {len(all_points)} points.")
print(f"Got {len(all_points_filtered)} filtered points")

Remapping

Our ZVP points might (due to the bounds above thus the incompleteness of our analysis) lead to more zeros than we attribute to them (in more configurations), i.e. “false negatives”. They might also be erroneous and not lead to zeros if the argument filter_nonhomo is False, as non-homogenous intermediate polynomials are not filtered out of the analysis. They introduce “false positives” but also some true positives.

Thus we perform a remapping step where we execute the scalar multiplication with given points and trace whether they introduce the zeros. This gives us a new distinguishing map remapped_hit_point_map, now without “false negatives” or “false positives”.

Note that we also construct two additional maps remapped_count_point_map and remapped_position_point_map which represent the results of remapping using an oracle counting the zeros and an oracle giving the positions of the zeros in the computation, respectively.

[ ]:
def remap(coords, chunk, points, scalar_map):
    import numpy as np
    from pyecsca.ec.mult import LTRMultiplier, AccumulationOrder
    from pyecsca.ec.context import local, DefaultContext
    from pyecsca.ec.error import UnsatisfiedAssumptionError
    from pyecsca.ec.formula import FormulaAction
    lp = len(points)
    lc = len(chunk)
    counts = np.full((lc, lp), -1, dtype=np.int16)
    positions = np.full((lc, lp), None, dtype=object)
    for i, formulas in enumerate(chunk):
        mult = LTRMultiplier(*formulas, None, False, AccumulationOrder.PeqRP, True, True)

        for j, entry in enumerate(points):
            if entry is None:
                continue
            point, params = entry
            mult.init(params, point)
            scalar = scalar_map[params]
            with local(DefaultContext()) as ctx:
                try:
                    mult.multiply(scalar)
                except UnsatisfiedAssumptionError:
                    continue

            zeros = []

            def callback(action):
                if isinstance(action, FormulaAction):
                    for intermediate in action.op_results:
                        zeros.append(int(intermediate.value) == 0)
            ctx.actions.walk(callback)
            count = sum(zeros)
            counts[i, j] = count
            positions[i, j] = tuple(zeros)
    return counts, positions

remapped_hit_point_map = {}
remapped_count_point_map = {}
remapped_position_point_map = {}
#all_points_list = list(all_points)
all_points_list = list(set((list(res[0])[0], res[1]) for res in all_points_filtered.values()))

with TaskExecutor(max_workers=30, mp_context=spawn_context) as pool, enable_spawn(remap) as remap_spawn:
    for coord_name, coords in tqdm(model.coordinates.items()):
        param_map = {}
        points_mapped = []
        scalars_mapped = {}
        mapped = 0
        for point, params in tqdm(all_points_list, desc=f"Map points to {coord_name}", leave=False):
            if params not in param_map:
                try:
                    param_map[params] = params.to_coords(coords)
                except UnsatisfiedAssumptionError:
                    param_map[params] = None
            if param_map[params] is None:
                points_mapped.append(None)
            else:
                mapped += 1
                points_mapped.append((point.to_model(param_map[params].curve.coordinate_model, param_map[params].curve), param_map[params]))
        print(f"{coord_name}: {mapped} are compatible. Remapping...")
        for params, scalar in scalar_map.items():
            if params not in param_map or param_map[params] is None:
                continue
            scalars_mapped[param_map[params]] = scalar

        pairs = list(product(*formula_groups[coords]))
        chunk_size = 10
        chunks = 0
        for i in trange(0, len(pairs), chunk_size, desc=f"Chunking {coord_name}", smoothing=0, leave=False):
            chunk = pairs[i:i+chunk_size]
            chunks += 1
            pool.submit_task(chunk,
                             remap_spawn, coords, chunk, points_mapped, scalars_mapped)

        for chunk, future in tqdm(pool.as_completed(), total=chunks, unit_scale=chunk_size, desc=f"Remapping {coord_name} ({len(pairs)} formula pairs)", leave=False, smoothing=0):
            error = future.exception()
            if error:
                print(j, error)
            else:
                counts, positions = future.result()
                for cfg, counts_row, positions_row in zip(chunk, counts, positions):
                    hit_set = set()
                    hit_map = {}
                    count_map = {}
                    position_map = {}
                    for pp_tuple, count, position in zip(all_points_list, counts_row, positions_row):
                        hit_map[pp_tuple] = None if count == -1 else (count > 0)
                        count_map[pp_tuple] = count
                        position_map[pp_tuple] = position
                    remapped_hit_point_map[cfg] = hit_map
                    remapped_count_point_map[cfg] = count_map
                    remapped_position_point_map[cfg] = position_map
        print(f"{coord_name}: Remapped.")

You can now store the (or load the previously computed) remmapped maps.

[ ]:
with open("remapped.pickle", "wb") as f:
    pickle.dump((remapped_hit_point_map, remapped_count_point_map, remapped_position_point_map), f)
[ ]:
with open("remapped.pickle", "rb") as f:
    remapped_hit_point_map, remapped_count_point_map, remapped_position_point_map = pickle.load(f)

Distinguishing map and distinguishing tree building

Finally, we can build a tree using the remapped distinguishing map. Let’s first wrap the raw mapping with a distinguishing Map object.

[ ]:
cfgs = set()
for coord_name, coords in model.coordinates.items():
    cfgs.update(product(*formula_groups[coords]))

param_categories = {
    "a=-1": ["projective-1"],
    "a=-3": ["projective-3", "jacobian-3", "xyzz-3"],
    "a=0": ["jacobian-0"],
    "generic": ["jacobian", "projective", "modified", "xyzz", "xz"],
    "b=0": ["w12-0"]
}
cfg_categories = {}
for name, coord_names in param_categories.items():
    category_cfgs = set()
    for coord_name in coord_names:
        coords = model.coordinates[coord_name]
        category_cfgs.update(filter(lambda cfg: cfg[0].coordinate_model == coords and cfg[1].coordinate_model == coords, cfgs))
    cfg_categories[name] = category_cfgs
category_map = {cfg: {"category": name} for name, category_cfgs in cfg_categories.items() for cfg in category_cfgs}
dmap_categories = Map.from_io_maps(cfgs, category_map)
[ ]:
dmap_remapped = Map.from_io_maps(cfgs, remapped_hit_point_map)
[ ]:
from copy import deepcopy
dmap_dedup = deepcopy(dmap_remapped)
dmap_dedup.deduplicate()
print(f"Rows before: {len(dmap_remapped.mapping)}, rows after deduplication: {len(dmap_dedup.mapping)}.")

Let’s now watch the tree getting built.

[ ]:
tree_categories = Tree.build(cfgs, dmap_categories)
tree_remapped = tree_categories.expand(dmap_dedup)

The tree is built, we can examine its properties, such as: - the total number of configurations (formula pairs), - its depth, - number of leaves (configurations that cannot be further distinguished), - average leaf size - mean result size (if this tree were to be used for reverse-enginnering).

[ ]:
print(tree_remapped.describe())
[ ]:
print(tree_remapped.render_basic())

We can also investigate the other oracles and the distinguishing trees they can build:

[ ]:
dmap_count = Map.from_io_maps(cfgs, remapped_count_point_map)
dmap_count.deduplicate()
[ ]:
print(dmap_count.describe())
[ ]:
tree_count = tree_categories.expand(dmap_count)
[ ]:
dmap_position = Map.from_io_maps(cfgs, remapped_position_point_map)
dmap_position.deduplicate()
[ ]:
print(dmap_position.describe())
[ ]:
tree_position = tree_categories.expand(dmap_position)
[ ]:
print("Zero hit")
print(tree_remapped.describe())
print("\nZero count")
print(tree_count.describe())
print("\nZero position")
print(tree_position.describe())

Evaluation

[ ]:
correct_rate, precise_rate, amount_rate, query_rate = eval_tree_symmetric1(cfgs, [tree_remapped], num_tries=100, num_cores=30)
[ ]:
store("zvp_re_symmetric.nc", correct_rate, precise_rate, amount_rate, query_rate)
[ ]:
correct_rate, precise_rate, amount_rate, query_rate = load("zvp_re_symmetric.nc")
[ ]:
success_rate_symmetric(correct_rate, None).savefig("zvp_re_success_rate_symmetric.pdf", bbox_inches="tight")
precise_rate_symmetric(precise_rate).savefig("zvp_re_precise_rate_symmetric.pdf", bbox_inches="tight")
query_rate_symmetric(query_rate).savefig("zvp_re_query_rate_symmetric.pdf", bbox_inches="tight")
amount_rate_symmetric(amount_rate).savefig("zvp_re_amount_rate_symmetric.pdf", bbox_inches="tight")
success_rate_vs_query_rate_symmetric(query_rate, correct_rate).savefig("zvp_re_scatter_symmetric.pdf", bbox_inches="tight")
success_rate_vs_majority_symmetric(correct_rate).savefig("zvp_re_plot_symmetric.pdf", bbox_inches="tight")
[ ]:
correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b = eval_tree_asymmetric1(cfgs, [tree_remapped], num_tries=100, num_cores=30)
[ ]:
store("zvp_re_asymmetric.nc", correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b)
[ ]:
correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b = load("zvp_re_asymmetric.nc")
[ ]:
success_rate_asymmetric(correct_rate_b, None).savefig("zvp_re_success_rate_asymmetric.pdf", bbox_inches="tight")
precise_rate_asymmetric(precise_rate_b).savefig("zvp_re_precise_rate_asymmetric.pdf", bbox_inches="tight")
query_rate_asymmetric(query_rate_b).savefig("zvp_re_query_rate_asymmetric.pdf", bbox_inches="tight")
amount_rate_asymmetric(amount_rate_b).savefig("zvp_re_amount_rate_asymmetric.pdf", bbox_inches="tight")
success_rate_vs_majority_asymmetric(correct_rate_b).savefig("zvp_re_plot_asymmetric.pdf", bbox_inches="tight")
[ ]:
correct_rate_c, precise_rate_c, amount_rate_c, query_rate_c = eval_tree_binomial1(cfgs, [tree_count], num_tries=100, num_cores=30)
[ ]:
store("zvp_re_binomial.nc", correct_rate_c, precise_rate_c, amount_rate_c, query_rate_c)
[ ]:
correct_rate_c, precise_rate_c, amount_rate_c, query_rate_c = load("zvp_re_binomial.nc")
[ ]:
success_rate_binomial(correct_rate_c, None).savefig("zvp_re_success_rate_binomial.pdf", bbox_inches="tight")
precise_rate_binomial(precise_rate_c).savefig("zvp_re_precise_rate_binomial.pdf", bbox_inches="tight")
query_rate_binomial(query_rate_c).savefig("zvp_re_query_rate_binomial.pdf", bbox_inches="tight")
amount_rate_binomial(amount_rate_c).savefig("zvp_re_amount_rate_binomial.pdf", bbox_inches="tight")

Factor sets

Now, lets examine the representation of factor sets and the distinguishing trees they can build.

[ ]:
fset_map = {}
fset_nonhomo_map = {}
factor_sets = {}
factor_sets_nonhomo = {}
for coord_name, coords in tqdm(model.coordinates.items()):
    for formula_group in formula_groups[coords]:
        for formula in tqdm(formula_group, leave=False):
            factor_sets[formula] = compute_factor_set(formula)
            factor_sets_nonhomo[formula] = compute_factor_set(formula, filter_nonhomo=False)
    formula_combinations = list(product(*formula_groups[coords]))
    for formulas in tqdm(formula_combinations, leave=False):
        fset = set()
        fset_nonhomo = set()
        for formula in formulas:
            fset.update(factor_sets[formula])
            fset_nonhomo.update(factor_sets_nonhomo[formula])
        fset_map[formulas] = fset
        fset_nonhomo_map[formulas] = fset_nonhomo
[ ]:
with open("factor_sets_extended.pickle", "wb") as f:
    pickle.dump((factor_sets, fset_map), f)
with open("factor_sets_nonhomo_extended.pickle", "wb") as f:
    pickle.dump((factor_sets_nonhomo, fset_nonhomo_map), f)
[ ]:
with open("factor_sets_extended.pickle", "rb") as f:
    factor_sets, fset_map = pickle.load(f)
with open("factor_sets_nonhomo_extended.pickle", "rb") as f:
    factor_sets_nonhomo, fset_nonhomo_map = pickle.load(f)

Build two trees, one with the filtered factor sets and one with unfilitered factor sets (with non-homogenous polynomials present).

[ ]:
dmap_fset = Map.from_sets(set(fset_map.keys()), fset_map, deduplicate=True)
[ ]:
print(dmap_fset.describe())
[ ]:
tree_fset = Tree.build(dmap_fset.cfgs, dmap_fset)
[ ]:
dmap_fset_nonhomo = Map.from_sets(set(fset_nonhomo_map.keys()), fset_nonhomo_map, deduplicate=True)
[ ]:
print(dmap_fset_nonhomo.describe())
[ ]:
tree_fset_nonhomo = Tree.build(dmap_fset_nonhomo.cfgs, dmap_fset_nonhomo)
[ ]:
print("Factor sets")
print(tree_fset.describe())
print("\nFactor sets (unfiltered)")
print(tree_fset_nonhomo.describe())

As we can see, our technique is able to distinguish formulas better than the filtered factor sets, but not better than unfiltered factor sets.

We can further analyze where unused possibilities of distinguishing remain but expanding the tree using the distinguishing map built out of factor sets.

[ ]:
expanded = tree_remapped.expand(dmap_fset)
[ ]:
print("Zero hit + factor sets")
print(expanded.describe())

The tree improves its distinguishing abilities, however, upon closer analysis, the polynomials that it uses to further split the configurations never actually have solutions on the curves with the assumed properties (e.g. a = 0).

Miscellaneous analysis

[ ]:
p = set()
for node in PreOrderIter(expanded.root):
    if isinstance(node.dmap_input, Poly):
        print(node.dmap_input, [len(child.cfgs) for child in node.children])
        print("\t", "\n\t".join(", ".join(f"({cfg[0]} {cfg[1]})" for cfg in child.cfgs) for child in node.children))
        p.add(node.dmap_input)
print("---")
for pp in p:
    print(pp)
    for formula, fset in factor_sets_nonhomo.items():
        if pp in fset:
            print(formula)
[ ]:
rev_point_map = {}
for (poly, params_cat, k), (points, affine_params) in all_points_filtered.items():
    for point in points:
        poly_set = rev_point_map.setdefault((point, affine_params), set())
        poly_set.add((poly, params_cat, k))

for (point, affine_params), poly_set in rev_point_map.items():
    if len(poly_set) > 1:
        print(point, affine_params.curve.parameters, f"p={affine_params.curve.prime}")
        cond = affine_params.name.split("=")[1].split("[")[0] if "=" in affine_params.name else ""
        polys_mapped = set()
        for poly, params_cat, k in poly_set:
            mapd = eliminate_y(poly, affine_params.curve.model)
            polys_mapped.add(mapd)
            print(poly.as_expr(), "|", params_cat, "|", k)
            #for formula, fset in factor_sets.items():
            #    if poly in fset and (cond in formula.coordinate_model.name or "-" not in formula.coordinate_model.name):
            #        print("\t", formula)
        print("->")
        for poly in polys_mapped:
            p = Poly(poly, domain=ZZ)
            print(p.factor_list())
        print("------")
    else:
        print(".", end="") #poly_set, point)
[ ]:
for node in PreOrderIter(tree_remapped.root):
    if node.dmap_input:
        pad = " " * node.depth
        poly_set = rev_point_map[node.dmap_input]
        point, affine_params = node.dmap_input
        print(pad, point, affine_params.name, affine_params.curve.parameters, f"p={affine_params.curve.prime}")
        chain = chain_map[affine_params]
        print(pad, chain)
        cond = affine_params.name.split("=")[1].split("[")[0] if "=" in affine_params.name else ""
        true_formulas = set()
        for poly, params_cat, k in poly_set:
            print(pad, poly.as_expr(), "|", params_cat, "|", k)
            for formula, fset in factor_sets.items():
                if poly in fset and (cond in formula.coordinate_model.name or "-" not in formula.coordinate_model.name):
                    kvar = ("add", (1, k)) if formula.shortname == "add" else ("dbl", (k, ))
                    print(pad, "\t", formula, "*" if kvar in chain else "not")
                    if kvar in chain:
                        true_formulas.add(formula)
        for child in node.children:
            print(pad, child.response)
            for cfg in child.cfgs:
                if child.response == True:
                    print(pad, "ok" if true_formulas.intersection(cfg) else "nok", cfg)
                elif child.response == False:
                    print(pad, "ok" if not true_formulas.intersection(cfg) else "nok", cfg)
                else:
                    print(pad, cfg)
            print("")
        print("------")
[ ]:
for coord_name, coords in model.coordinates.items():
    total = 1
    print(coord_name)
    for formula_group, formula_class in zip(formula_groups[coords], formula_classes):
        lf = len(formula_group)
        print("\t", formula_class.shortname, lf)
        for formula in formula_group:
            print("\t", formula)
        total *= lf
    print("\t", total)
[ ]:
for tree in (tree_remapped, tree_count, tree_position):
    same = 0
    for leaf in tree.leaves:
        cds = set()
        for c in leaf.cfgs:
            cds.add(c[0].coordinate_model)
            cds.add(c[1].coordinate_model)
        same += len(cds) == 1
    print(same / len(tree.leaves))
[ ]:
tree_counter = Counter()
for node in LevelOrderIter(tree_count.root):
    if isinstance(node.response, int):
        tree_counter[node.response] += 1
    if node.children:
        print(list(map(lambda n: n.response, node.children)))
fig, ax = plt.subplots(figsize=(10, 3))
ax.bar(tree_counter.keys(), tree_counter.values(), width=1);
[ ]:
dmap_counter = Counter()
for col in dmap_count.mapping:
    dmap_counter.update(dmap_count.mapping[col])

del dmap_counter[-1]
del dmap_counter[0]
fig, ax = plt.subplots(figsize=(10, 3))
ax.bar(dmap_counter.keys(), dmap_counter.values(), width=1);
[ ]: