#!/usr/bin/env python
# coding: utf-8
# # ZVP-based reverse-engineering
# This notebook showcases the ZVP-based reverse-engineering technique for addition formulas.
#
# - [Exploration](#Exploration)
# - [Reverse-engineering](#Reverse-engineering)
# - [Computing addition chains](#Computing-addition-chains)
# - [Loading the formulas](#Loading-the-formulas)
# - [Computing the factor sets](#Computing-the-factor-sets)
# - [Accumulating polynomials](#Accumulating-polynomials)
# - [Computing ZVP points](#Computing-ZVP-points)
# - [Remapping](#Remapping)
# - [Distinguishing map and distinguishing tree building](#Distinguishing-map-and-distinguishing-tree-building)
# - [Evaluation](#Evaluation)
#
# In[ ]:
import io
import holoviews as hv
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, symbols, Poly
from contextlib import contextmanager
from importlib import import_module, invalidate_caches
from pathlib import Path
from itertools import product
from IPython.display import display
from tqdm.notebook import tqdm, trange
from anytree import LevelOrderIter, PreOrderIter
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 AdditionFormula, DoublingFormula
from pyecsca.ec.point import Point
from pyecsca.ec.mod import mod, SymbolicMod
from pyecsca.ec.mult import LTRMultiplier, AccumulationOrder
from pyecsca.ec.error import UnsatisfiedAssumptionError
from pyecsca.sca.re.tree import Map, Tree
from pyecsca.sca.re.zvp import zvp_points, compute_factor_set, 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")
get_ipython().run_line_magic('matplotlib', 'ipympl')
hv.extension("bokeh")
# In[ ]:
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).
# In[ ]:
which = "jacobian"
coords = model.coordinates[which]
# In[ ]:
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.
# In[ ]:
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.
# In[ ]:
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.
# In[ ]:
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](formulas.ipynb) notebook.
# In[ ]:
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.
#
#
#
# This is a resource intensive cell that uses parallelism. Set the `num_workers` variable to something reasonable, like the number of cores of your machine minus two.
#
#
# In[ ]:
num_workers = 30
# In[ ]:
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=num_workers, 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.
# In[ ]:
with open("factor_sets.pickle", "wb") as f:
pickle.dump(factor_sets, f)
print(f"Stored {len(factor_sets)} factor sets.")
# In[ ]:
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.
# In[ ]:
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 the polynomials. Also filter the points such that for each "polynomial, curve category, k" we have only one point, as more are unnecessary.
#
#
#
# This is a resource intensive cell that uses parallelism. Recall the `num_workers` variable.
#
#
# In[ ]:
# 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=num_workers, 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(r"(.+)\[([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.
# In[ ]:
with open("all_points.pickle", "wb") as f:
pickle.dump((all_points, all_points_filtered), f)
# In[ ]:
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.
#
#
#
# This is a resource intensive cell that uses parallelism. Recall the `num_workers` variable.
#
#
# In[ ]:
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[0].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=num_workers, 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(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.
# In[ ]:
with open("remapped.pickle", "wb") as f:
pickle.dump((remapped_hit_point_map, remapped_count_point_map, remapped_position_point_map), f)
# In[ ]:
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.
# In[ ]:
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)
# In[ ]:
dmap_remapped = Map.from_io_maps(cfgs, remapped_hit_point_map)
# In[ ]:
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.
# In[ ]:
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).
# In[ ]:
print(tree_remapped.describe())
# In[ ]:
print(tree_remapped.render_basic())
# We can also investigate the other oracles and the distinguishing trees they can build:
# In[ ]:
dmap_count = Map.from_io_maps(cfgs, remapped_count_point_map)
dmap_count.deduplicate()
# In[ ]:
print(dmap_count.describe())
# In[ ]:
tree_count = tree_categories.expand(dmap_count)
# In[ ]:
dmap_position = Map.from_io_maps(cfgs, remapped_position_point_map)
dmap_position.deduplicate()
# In[ ]:
print(dmap_position.describe())
# In[ ]:
tree_position = tree_categories.expand(dmap_position)
# In[ ]:
print("Zero hit")
print(tree_remapped.describe())
print("\nZero count")
print(tree_count.describe())
print("\nZero position")
print(tree_position.describe())
# ### Evaluation
#
# The following cells evaluate the ZVP-RE method and its behavior under various levels of noise. The cells with the `store` and `load` calls can be used to store the results so that different plots can be printed without re-running the evaluation.
#
#
#
# These are resource intensive cells that use parallelism. Recall the `num_workers` variable.
#
#
# First, lets look at the binary oracle in presence of symmetric noise.
# In[ ]:
correct_rate, precise_rate, amount_rate, query_rate = eval_tree_symmetric1(cfgs, [tree_remapped], num_tries=100, num_cores=num_workers)
# Let's save the results for later and then plot them.
# In[ ]:
store("zvp_re_symmetric.nc", correct_rate, precise_rate, amount_rate, query_rate)
# In[ ]:
correct_rate, precise_rate, amount_rate, query_rate = load("zvp_re_symmetric.nc")
# In[ ]:
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")
# Now, lets look at the binary oracle in presence of asymmetric noise.
# In[ ]:
correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b = eval_tree_asymmetric1(cfgs, [tree_remapped], num_tries=100, num_cores=num_workers)
# Let's save the results for later and then plot them.
# In[ ]:
store("zvp_re_asymmetric.nc", correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b)
# In[ ]:
correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b = load("zvp_re_asymmetric.nc")
# In[ ]:
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")
# Finally, lets analyze the count oracle in the presence of noise.
# In[ ]:
correct_rate_c, precise_rate_c, amount_rate_c, query_rate_c = eval_tree_binomial1(cfgs, [tree_count], num_tries=100, num_cores=num_workers)
# Let's save the results for later and then plot them.
# In[ ]:
store("zvp_re_binomial.nc", correct_rate_c, precise_rate_c, amount_rate_c, query_rate_c)
# In[ ]:
correct_rate_c, precise_rate_c, amount_rate_c, query_rate_c = load("zvp_re_binomial.nc")
# In[ ]:
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.
# In[ ]:
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
# In[ ]:
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)
# In[ ]:
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).
# In[ ]:
dmap_fset = Map.from_sets(set(fset_map.keys()), fset_map, deduplicate=True)
# In[ ]:
print(dmap_fset.describe())
# In[ ]:
tree_fset = Tree.build(dmap_fset.cfgs, dmap_fset)
# In[ ]:
dmap_fset_nonhomo = Map.from_sets(set(fset_nonhomo_map.keys()), fset_nonhomo_map, deduplicate=True)
# In[ ]:
print(dmap_fset_nonhomo.describe())
# In[ ]:
tree_fset_nonhomo = Tree.build(dmap_fset_nonhomo.cfgs, dmap_fset_nonhomo)
# In[ ]:
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.
# In[ ]:
expanded = tree_remapped.expand(dmap_fset)
# In[ ]:
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`).
# In[ ]: