Merge branch 'zvp-re'

3 files changed, 878 insertions, 0 deletions

diff --git a/re/epa.ipynb b/re/epa.ipynb
new file mode 100644
index 0000000..67e7973
--- /dev/null
+++ b/re/epa.ipynb
@@ -0,0 +1,262 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "18343664-ebec-4e95-88f2-0231082b6b6e",
+   "metadata": {},
+   "source": [
+    "# EPA-based reverse engineering"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dac3b016-508b-4e63-897a-8ab46cb4f6ed",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import io\n",
+    "import tabulate\n",
+    "from itertools import product\n",
+    "from IPython.display import HTML, display\n",
+    "\n",
+    "from pyecsca.ec.model import ShortWeierstrassModel\n",
+    "from pyecsca.ec.coordinates import AffineCoordinateModel\n",
+    "from pyecsca.ec.curve import EllipticCurve\n",
+    "from pyecsca.ec.params import DomainParameters, load_params_ectester\n",
+    "from pyecsca.ec.mod import Mod, miller_rabin, gcd\n",
+    "from pyecsca.ec.point import Point\n",
+    "from pyecsca.ec.error import NonInvertibleError\n",
+    "from pyecsca.ec.mult import LTRMultiplier\n",
+    "from pyecsca.sca.re.rpa import MultipleContext\n",
+    "from pyecsca.ec.context import local\n",
+    "from pyecsca.ec.error import UnsatisfiedAssumptionError"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9506f889-75cf-4b27-8535-e71cb8d5301e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "curves = [\n",
+    "    # phi(p)/p =\n",
+    "    # 0.8286039438617044\n",
+    "    \"dfb2da5e1b7bd7bb098cb975966293ed,d9c4372806e8131b18d0036e8f832749,bcae41be8e808acdc04bb769dead91e2,0e2f983c0f852bef381f567448f0d488,1599bba77ed1cb8dec41555098958492,10fcabd48fffc71e6300d44acc236157d,0001\",\n",
+    "    # 0.633325621696952\n",
+    "    \"cca6f6718a06cad7094962b2a35f067d,67aa9464eb493fbb7b509d29381b9a9d,cafc69aa517b654a6a608644996cc8d1,4c092beb06cc00751eec39675f680cb8,82800378a47dd6f26ff6a50f69e4c4e6,18a22d20b6de3ff6bdc49329c21163f77,0001\",\n",
+    "    # 0.8508806646440022 \n",
+    "    \"b3755d654bad73114e4191e9f5f36af9,9fe4f88cfbacba71f4b767ace8580c74,4610526fdcfbd69aed453ac2ee6efeef,542d8e0bbafe40dae36f25cbc350add6,68a65f5a5dc304bfd0d8fe963c250206,118a34a1ea295e78b3a3c960b6f680ee1,0001\",\n",
+    "    # 0.9845701775215489    \n",
+    "    \"de1406450d5d7e91d81907956019c0c1,5fbe46b9f1086011e18f5d823c6110ce,a859c36ceeadb39c7a978f7b1b0563ee,1cba89c3f099c29401ecf3fe1806e822,345d7282a0114070be91f95fe3db1faa,0fcd24d24e57a40547814b6766b9ea735,0001\",\n",
+    "    # 0.9547100843537808\n",
+    "    \"f1a8a441b6d0e9600e33ccf16f9b8291,b3f55185bd6a63528e3d560c6a7b729a,c2fee2d65350e870eda0ac5e2b96b810,29b3e793822fad03a3c2ebca3cf62c12,b937d5389b6c5d0212d0f53e26843092,1153442389f9e1da8dd130bc93c6ef42b,0001\",\n",
+    "    # 0.7214369438844093\n",
+    "    \"a4dfa4b6b065c40b45980474266c9fbb,2c3486e725755b44a7c119473c5b9c64,329078ab070fc18edc6ce53047e00a39,9f6209be91b66943d9e8e0b61c4aae4e,05271c9ac628351b9add9e1be69a9fa4,0cef2e52ffe86ebc6dd323912ac7d9a87,0001\",\n",
+    "    # 0.4716485170445178 (160-bit)\n",
+    "    \"db49063db56b7783fa01dd62077c5a88dfa28009,aee572fdd4790bcd4729bb3b612b52a573df46e9,dab9e68366a593ca1df9cb2f20890a578729d6ef,d4a3aaf43bdb25be7c308b69ae54f639e6e32e8c,7b6c82140bb427ac6e2a64507f60775949b2c8ce,34a9fbe62b272f930b2e5027780a32300feb0dd8f,0001\"\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e92cdd8-9040-4c6d-9700-f992d09195ee",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = ShortWeierstrassModel()\n",
+    "affine = AffineCoordinateModel(model)\n",
+    "which = \"projective\"\n",
+    "coords = model.coordinates[which]\n",
+    "\n",
+    "params = load_params_ectester(io.BytesIO(curves[3].encode()), which)\n",
+    "curve = params.curve\n",
+    "g = params.generator\n",
+    "n = params.order"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "703da72a-1c41-41a9-886d-453da1160932",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "adds = list(filter(lambda formula: formula.name.startswith(\"add\"), coords.formulas.values()))\n",
+    "dbls = list(filter(lambda formula: formula.name.startswith(\"dbl\"), coords.formulas.values()))\n",
+    "formula_pairs = list(product(adds, dbls))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ffce5b71-3029-4219-a45a-1c8d78748fee",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "scalars = [int(Mod.random(n)) for _ in range(200)]\n",
+    "\n",
+    "table = [[\"Pair\", \"scalars with trivial gcd\", \"scalars with all multiples with trivial gcds\", \"scalars with invertible final zs\", \"scalars with all multiples's zs invertible\"]]\n",
+    "pair_table = [[None for _ in dbls] for _ in adds]\n",
+    "for pair in formula_pairs:\n",
+    "    mult = LTRMultiplier(*pair, None, False, True, True)\n",
+    "    mult.init(params, g)\n",
+    "    inv = 0\n",
+    "    gcds = 0\n",
+    "    fgcds = 0\n",
+    "    zs = 0\n",
+    "    for scalar in scalars:\n",
+    "        if gcd(scalar, n) == 1:\n",
+    "            gcds += 1\n",
+    "        with local(MultipleContext()) as ctx:\n",
+    "            res = mult.multiply(scalar)\n",
+    "        if all(map(lambda x: gcd(x, n) == 1, ctx.points.values())):\n",
+    "            fgcds += 1\n",
+    "        if all(map(lambda x: gcd(int(x.Z), n) == 1, ctx.points.keys())):\n",
+    "            zs += 1\n",
+    "    \n",
+    "        try:\n",
+    "            res.Z.inverse()\n",
+    "            inv += 1\n",
+    "        except NonInvertibleError as e:\n",
+    "            pass\n",
+    "    #print(\"---\")\n",
+    "    #print(pair[0].name, pair[1].name)\n",
+    "    #print(gcds, \"\\tscalars with trivial gcd\")\n",
+    "    #print(fgcds, \"\\tscalars with all multiples with trivial gcds\")\n",
+    "    #print(inv, \"\\tscalars with invertible final zs\")\n",
+    "    #print(zs, \"\\tscalars with all multiples's zs invertible\")\n",
+    "    pair_table[adds.index(pair[0])][dbls.index(pair[1])] = inv\n",
+    "    table.append([f\"{pair[0].name}, {pair[1].name}\", gcds, fgcds, inv, zs])\n",
+    "for pl, add in zip(pair_table, adds):\n",
+    "    pl.insert(0, add.name)\n",
+    "pair_table.insert(0, [None] + [dbl.name for dbl in dbls])\n",
+    "display(HTML(tabulate.tabulate(table, tablefmt=\"html\", headers=\"firstrow\")))\n",
+    "display(HTML(tabulate.tabulate(pair_table, tablefmt=\"html\", headers=\"firstrow\")))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f47c0743-a6cf-408a-8c6e-dca6732278e2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def simulate_epa_oracle(affine_params, affine_point, scalar):\n",
+    "    real_coords = model.coordinates[\"projective\"]\n",
+    "    real_add = real_coords.formulas[\"add-2007-bl\"]\n",
+    "    real_dbl = real_coords.formulas[\"dbl-2007-bl\"]\n",
+    "    real_mult = LTRMultiplier(real_add, real_dbl, None, False, True, True)\n",
+    "    params = affine_params.to_coords(real_coords)\n",
+    "    point = affine_point.to_model(real_coords, params.curve)\n",
+    "    real_mult.init(params, point)\n",
+    "    res = real_mult.multiply(scalar)\n",
+    "    try:\n",
+    "        res.to_affine()\n",
+    "        return True\n",
+    "    except NonInvertibleError as e:\n",
+    "        return False\n",
+    "\n",
+    "def epa_distinguish(oracle, mult_factory):\n",
+    "    affine_params = load_params_ectester(io.BytesIO(curves[3].encode()), \"affine\")\n",
+    "    model = affine_params.curve.model\n",
+    "    scalars = [int(Mod.random(affine_params.order)) for _ in range(100)]\n",
+    "    responses = [oracle(affine_params,  affine_params.generator, scalar) for scalar in scalars]\n",
+    "    candidates = set()\n",
+    "    print(\"Got responses\")\n",
+    "    total = 0\n",
+    "    for coords in model.coordinates.values():\n",
+    "        adds = list(filter(lambda formula: formula.name.startswith(\"add\"), coords.formulas.values()))\n",
+    "        dbls = list(filter(lambda formula: formula.name.startswith(\"dbl\"), coords.formulas.values()))\n",
+    "        formula_pairs = list(product(adds, dbls))\n",
+    "        total += len(formula_pairs)\n",
+    "        try:\n",
+    "            params = affine_params.to_coords(coords)\n",
+    "        except UnsatisfiedAssumptionError:\n",
+    "            print(f\"Skipping {coords.name}, does not fit\")\n",
+    "            continue\n",
+    "\n",
+    "        for pair in formula_pairs:\n",
+    "            mult = mult_factory(*pair)\n",
+    "            mult.init(params, params.generator)\n",
+    "            abort = False\n",
+    "            print(f\"Trying {coords.name} {pair[0].name} {pair[1].name}\", end=\"\") \n",
+    "            for scalar, target in zip(scalars, responses):\n",
+    "                res = mult.multiply(scalar)\n",
+    "                try:\n",
+    "                    res.to_affine()\n",
+    "                    if not target:\n",
+    "                        # not this one\n",
+    "                        abort = True\n",
+    "                        break\n",
+    "                except NonInvertibleError as e:\n",
+    "                    if target:\n",
+    "                        # not this one\n",
+    "                        abort = True\n",
+    "                        break\n",
+    "            if abort:\n",
+    "                print(\" not\")\n",
+    "                continue\n",
+    "            else:\n",
+    "                print(\" candidate\")\n",
+    "                candidates.add((coords.name, pair[0].name, pair[1].name))\n",
+    "    print(f\"Got {len(candidates)} out of {total} total\")\n",
+    "    return candidates\n",
+    "                "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d2c26c45-a66a-4472-830e-b22010a967e5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c = epa_distinguish(simulate_epa_oracle, lambda add,dbl:LTRMultiplier(add, dbl, None, False, True, True))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dd3e076c-cd52-47e1-beb0-2cb59812dacc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "aa548870-4236-42d8-a4a6-b6b8082b87bb",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/re/rpa.ipynb b/re/rpa.ipynb
new file mode 100644
index 0000000..29bf15e
--- /dev/null
+++ b/re/rpa.ipynb
@@ -0,0 +1,373 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "66c30004-cd2c-4d34-9999-f33f9e6fd5e9",
+   "metadata": {},
+   "source": [
+    "# RPA-based reverse-engineering"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11a5f41d-1471-49c3-ba7c-ac87470a31d0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from collections import Counter\n",
+    "import numpy as np\n",
+    "import holoviews as hv\n",
+    "from scipy.signal import find_peaks\n",
+    "\n",
+    "from pyecsca.ec.model import ShortWeierstrassModel\n",
+    "from pyecsca.ec.coordinates import AffineCoordinateModel\n",
+    "from pyecsca.ec.curve import EllipticCurve\n",
+    "from pyecsca.ec.params import DomainParameters\n",
+    "from pyecsca.ec.formula import FormulaAction\n",
+    "from pyecsca.ec.point import Point\n",
+    "from pyecsca.ec.mod import Mod\n",
+    "from pyecsca.ec.mult import *\n",
+    "from pyecsca.sca.trace.sampling import downsample_average, downsample_max\n",
+    "from pyecsca.sca.trace.process import normalize, rolling_mean\n",
+    "from pyecsca.sca.trace.combine import average, subtract\n",
+    "from pyecsca.sca.trace.test import welch_ttest\n",
+    "from pyecsca.sca.attack.leakage_model import HammingWeight, NormalNoice\n",
+    "from pyecsca.ec.context import DefaultContext, local\n",
+    "from pyecsca.sca.re.rpa import MultipleContext, rpa_distinguish\n",
+    "from pyecsca.sca.trace import Trace\n",
+    "from pyecsca.sca.trace.plot import plot_trace, plot_traces"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3ff1d591-f922-4c94-9e47-ab053fc21cf1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hv.extension(\"bokeh\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "194fff59-1c4b-473a-9ffc-99b256aecc24",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = ShortWeierstrassModel()\n",
+    "coordsaff = AffineCoordinateModel(model)\n",
+    "coords = model.coordinates[\"projective\"]\n",
+    "add = coords.formulas[\"add-2007-bl\"]  # The formulas are irrelevant for this method\n",
+    "dbl = coords.formulas[\"dbl-2007-bl\"]\n",
+    "neg = coords.formulas[\"neg\"]\n",
+    "\n",
+    "# A 64-bit prime order curve for testing things out\n",
+    "p = 0xc50de883f0e7b167\n",
+    "a = Mod(0x4833d7aa73fa6694, p)\n",
+    "b = Mod(0xa6c44a61c5323f6a, p)\n",
+    "gx = Mod(0x5fd1f7d38d4f2333, p)\n",
+    "gy = Mod(0x21f43957d7e20ceb, p)\n",
+    "n = 0xc50de885003b80eb\n",
+    "h = 1\n",
+    "\n",
+    "# A (0, y) RPA point on the above curve, in affine coords.\n",
+    "P0_aff = Point(coordsaff, x=Mod(0, p), y=Mod(0x1742befa24cd8a0d, p))\n",
+    "\n",
+    "infty = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))\n",
+    "g = Point(coords, X=gx, Y=gy, Z=Mod(1, p))\n",
+    "\n",
+    "curve = EllipticCurve(model, coords, p, infty, dict(a=a,b=b))\n",
+    "params = DomainParameters(curve, g, n, h)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "29fb8683-bcad-4a3e-869b-95f0e4b7bde3",
+   "metadata": {},
+   "source": [
+    "## Exploration\n",
+    "First select a bunch of multipliers. We will be trying to distinguish among these."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "febb198f-4370-4abd-8edc-17c5c1da8d0f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "multipliers = []\n",
+    "multipliers.append(LTRMultiplier(add, dbl, None, False, True, True))\n",
+    "multipliers.append(LTRMultiplier(add, dbl, None, True, True, True))\n",
+    "multipliers.append(RTLMultiplier(add, dbl, None, False, True))\n",
+    "multipliers.append(RTLMultiplier(add, dbl, None, True, True))\n",
+    "multipliers.append(SimpleLadderMultiplier(add, dbl, None, True, True))\n",
+    "multipliers.append(BinaryNAFMultiplier(add, dbl, neg, None, True))\n",
+    "multipliers.append(WindowNAFMultiplier(add, dbl, neg, 3, None, True))\n",
+    "multipliers.append(WindowNAFMultiplier(add, dbl, neg, 4, None, True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bea70d56-1359-423b-9273-fae5877f6400",
+   "metadata": {},
+   "source": [
+    "Then select a random scalar and simulate computation using all of the multipliers, track the multiples, print the projective and affine results."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "44e7e2e9-e0ad-4c8d-8605-af68f73d73e2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "scalar = 0b1000000000000000000000000000000000000000000000000\n",
+    "scalar = 0b1111111111111111111111111111111111111111111111111\n",
+    "scalar = 0b1010101010101010101010101010101010101010101010101\n",
+    "scalar = 0b1111111111111111111111110000000000000000000000000\n",
+    "scalar = 123456789123456789\n",
+    "# multiples is a mapping from a multiple (integer) to a set of scalar multipliers that compute said multiple when doing [scalar]\n",
+    "multiples = {}\n",
+    "\n",
+    "for mult in multipliers:\n",
+    "    print(repr(mult))\n",
+    "    with local(MultipleContext()) as ctx:\n",
+    "        mult.init(params, g)\n",
+    "        res = mult.multiply(scalar)\n",
+    "    print(res, res.to_affine())\n",
+    "    for m in ctx.points.values():\n",
+    "        s = multiples.setdefault(m, set())\n",
+    "        s.add(mult)\n",
+    "    print()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ba284641-c29b-4e42-95ec-aa630e305b10",
+   "metadata": {},
+   "source": [
+    "Pick a multiple `k` that is computed by some multiplier for the scalar,\n",
+    "invert it mod n, and do `[k^-1]P0` to obtain a point `P0_target`,\n",
+    "such that, `[k]P0_target = P0` and `P0` has a zero coordinate."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a7fb8a3f-7938-493b-88dc-582ba4d8959d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "k = 7186132\n",
+    "kinv = Mod(k, n).inverse()\n",
+    "P0_target = curve.affine_multiply(P0_aff, int(kinv)).to_model(coords, curve)\n",
+    "\n",
+    "print(\"Original P0\", P0_aff)\n",
+    "print(\"P0_target  \", P0_target.to_affine())\n",
+    "print(\"Verify P0  \", curve.affine_multiply(P0_target.to_affine(), k))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "914a67a5-b6a2-4d02-811f-a423099853c3",
+   "metadata": {},
+   "source": [
+    "Now go over the multipliers with P0_target and the original scalar as input.\n",
+    "Then look whether a zero coordinate point was computed.\n",
+    "Also look at whether the multiple \"k\" was computed. These two should be the same."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8113cb3f-dc06-4cb7-955c-11cedb4fbdd7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for mult in multipliers:\n",
+    "    print(mult.__class__.__name__)\n",
+    "    with local(MultipleContext()) as ctx:\n",
+    "        mult.init(params, P0_target)\n",
+    "        res = mult.multiply(scalar)\n",
+    "    print(\"\\tzero present     \", any(map(lambda P: P.X == 0, ctx.points.keys())))\n",
+    "    print(\"\\tmultiple computed\", k in ctx.points.values())\n",
+    "    print()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0d2b9fe8-5064-4887-b782-dcfe9f42d217",
+   "metadata": {},
+   "source": [
+    "Now lets look at the relation of multiples to multipliers."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "67d7705c-6a41-47d9-ad1e-23ea549aaf00",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for multiple, mults in multiples.items():\n",
+    "    print(multiple, [mult.__class__.__name__ for mult in mults])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e9015f3c-fba6-4614-b722-848b8522d072",
+   "metadata": {},
+   "source": [
+    "## Reverse-engineering\n",
+    "The `simulated_oracle` function simulates an RPA oracle that detect a zero coordinate point in the scalar multiplication.\n",
+    "This can be used by the `rpa_distinguish` function to distinguish the true scalar multiplier."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9bb61ac5-d837-4287-a5de-a9a63c346acf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def simulated_oracle(scalar, affine_point):\n",
+    "    real_mult = BinaryNAFMultiplier(add, dbl, neg, None, True)\n",
+    "    point = affine_point.to_model(params.curve.coordinate_model, params.curve)\n",
+    "    with local(MultipleContext()) as ctx:\n",
+    "        real_mult.init(params, point)\n",
+    "        real_mult.multiply(scalar)\n",
+    "    return any(map(lambda P: P.X == 0 or P.Y == 0, ctx.points.keys()))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "265eb13a-6028-4fde-9c3f-cc23768ba63e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rpa_distinguish(params, multipliers, simulated_oracle)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "62b0b8f2-8149-4abd-9aa7-a056b237ac6e",
+   "metadata": {},
+   "source": [
+    "The `simulate_trace` function simulates a Hamming weight leakage trace of a given multiplier computing a scalar multiple.\n",
+    "This is used by the `simulated_rpa_trace` function that does the RPA attack on simulated traces and returns the differential\n",
+    "trace. This is in turn used to build the `simulated_rpa_oracle` which can be used by the `rpa_distinguish` function to perform\n",
+    "RPA-RE and distinguish the true scalar multiplier."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96bec03e-5397-440b-9e8c-81ba5253921b",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "def simulate_trace(mult, scalar, point):\n",
+    "    with local(DefaultContext()) as ctx:\n",
+    "        mult.init(params, point)\n",
+    "        mult.multiply(scalar)\n",
+    "\n",
+    "    lm = HammingWeight()\n",
+    "    trace = []\n",
+    "\n",
+    "    def callback(action):\n",
+    "        if isinstance(action, FormulaAction):\n",
+    "            for intermediate in action.op_results:\n",
+    "                leak = lm(intermediate.value)\n",
+    "                trace.append(leak)\n",
+    "\n",
+    "    ctx.actions.walk(callback)\n",
+    "    return Trace(np.array(trace))\n",
+    "\n",
+    "def simulated_rpa_trace(mult, scalar, affine_point, noise):\n",
+    "    target_point = affine_point.to_model(params.curve.coordinate_model, params.curve)\n",
+    "    random_point = params.curve.affine_random().to_model(params.curve.coordinate_model, params.curve)\n",
+    "\n",
+    "    random_traces = [noise(simulate_trace(mult, scalar, random_point)) for _ in range(10)]\n",
+    "    target_traces = [noise(simulate_trace(mult, scalar, target_point)) for _ in range(500)]\n",
+    "\n",
+    "    random_avg = average(*random_traces)\n",
+    "    target_avg = average(*target_traces)\n",
+    "\n",
+    "    diff_trace = downsample_max(subtract(random_avg, target_avg), 25)\n",
+    "    return diff_trace\n",
+    "\n",
+    "def simulated_rpa_oracle(scalar, affine_point):\n",
+    "    real_mult = BinaryNAFMultiplier(add, dbl, neg, None, True)\n",
+    "    noise = NormalNoice(0, 1)\n",
+    "    diff_trace = normalize(simulated_rpa_trace(real_mult, scalar, affine_point, noise))\n",
+    "    peaks, props = find_peaks(diff_trace.samples, height=4)\n",
+    "    return len(peaks) != 0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d74a944f-e0a6-434d-8a12-138dfb9d516e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rpa_distinguish(params, multipliers, simulated_rpa_oracle)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e694b3b3-290d-4528-a611-16a183662944",
+   "metadata": {},
+   "source": [
+    "Note that the oracle function above has several parameters, like noise amplitude, amount of traces simulated and peak finding height threshold. The cell below compares the differential RPA trace when the multiple is computed in the simulation vs when it is not."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e47a70d7-7993-4560-9257-dbbc1cf658cf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "diff_real = normalize(simulated_rpa_trace(BinaryNAFMultiplier(add, dbl, neg, None, True), scalar, P0_target.to_affine(), NormalNoice(0, 1)))\n",
+    "diff_nothing = normalize(simulated_rpa_trace(LTRMultiplier(add, dbl, None, False, True, True), scalar, P0_target.to_affine(), NormalNoice(0, 1)))\n",
+    "plot_traces(diff_real, diff_nothing).opts(width=950, height=600)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "51a78d40-0cc1-49e2-9753-ddf49a496ed9",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/re/zvp.ipynb b/re/zvp.ipynb
new file mode 100644
index 0000000..1df4dbb
--- /dev/null
+++ b/re/zvp.ipynb
@@ -0,0 +1,243 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "3ff024c8-7e90-4094-ab07-a051a601d9cd",
+   "metadata": {},
+   "source": [
+    "# ZVP-based reverse-engineering"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "17a9e580-5f1e-45c9-8afd-fc35e833b0b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sympy\n",
+    "from sympy import FF, sympify, symbols\n",
+    "\n",
+    "from pyecsca.ec.model import ShortWeierstrassModel\n",
+    "from pyecsca.ec.coordinates import AffineCoordinateModel\n",
+    "from pyecsca.ec.curve import EllipticCurve\n",
+    "from pyecsca.ec.params import DomainParameters\n",
+    "from pyecsca.ec.point import Point\n",
+    "from pyecsca.ec.mod import Mod, SymbolicMod\n",
+    "from pyecsca.ec.context import DefaultContext, local\n",
+    "from pyecsca.misc.cfg import getconfig"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "d8ea0c4d-86e1-46af-ac20-569e6ef5439d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cfg = getconfig()\n",
+    "cfg.ec.mod_implementation = \"python\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "8de124f1-9498-4e55-8b7f-1bd291ccf3fe",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = ShortWeierstrassModel()\n",
+    "coordsaff = AffineCoordinateModel(model)\n",
+    "coords = model.coordinates[\"projective\"]\n",
+    "add = coords.formulas[\"add-2007-bl\"]\n",
+    "dbl = coords.formulas[\"dbl-2007-bl\"]\n",
+    "neg = coords.formulas[\"neg\"]\n",
+    "\n",
+    "# A 64-bit prime order curve for testing things out\n",
+    "p = 0xc50de883f0e7b167\n",
+    "a = Mod(0x4833d7aa73fa6694, p)\n",
+    "b = Mod(0xa6c44a61c5323f6a, p)\n",
+    "gx = Mod(0x5fd1f7d38d4f2333, p)\n",
+    "gy = Mod(0x21f43957d7e20ceb, p)\n",
+    "n = 0xc50de885003b80eb\n",
+    "h = 1\n",
+    "\n",
+    "field = FF(p)\n",
+    "\n",
+    "infty = Point(coords, X=Mod(0, p), Y=Mod(1, p), Z=Mod(0, p))\n",
+    "g = Point(coords, X=gx, Y=gy, Z=Mod(1, p))\n",
+    "\n",
+    "curve = EllipticCurve(model, coords, p, infty, dict(a=a,b=b))\n",
+    "params = DomainParameters(curve, g, n, h)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "id": "e91d81b6-a3ce-4b1f-8ace-3e55e0d6fadf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# add\n",
+    "# TODO: this is wrooooong, needs to have X1 and X2...\n",
+    "P1 = Point(coords, **{str(var): SymbolicMod(symbols(var + \"1\"), p) for var in coords.variables})\n",
+    "P2 = Point(coords, **{str(var): SymbolicMod(symbols(var + \"2\"), p) for var in coords.variables})\n",
+    "parameters = {str(var): SymbolicMod(var, p) for var in symbols(model.parameter_names)}\n",
+    "with local(DefaultContext()) as ctx:\n",
+    "    res = add(p, P1, P2, **parameters)[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "id": "7f4c7416-49a9-4fe5-a1ff-38f2cb6c5815",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle 2 Z_{1} Z_{2} \\left(Y_{1} Z_{2} + Y_{2} Z_{1}\\right) \\left(Z_{1}^{2} Z_{2}^{2} \\left(Y_{1} Z_{2} + Y_{2} Z_{1}\\right)^{4} + \\left(X_{1} Z_{2} + X_{2} Z_{1}\\right)^{2} - \\left(X_{1} Z_{2} + X_{2} Z_{1} + Z_{1} Z_{2} \\left(Y_{1} Z_{2} + Y_{2} Z_{1}\\right)^{2}\\right)^{2} + 2 \\left(- X_{1} X_{2} Z_{1} Z_{2} + Z_{1}^{2} Z_{2}^{2} a + \\left(X_{1} Z_{2} + X_{2} Z_{1}\\right)^{2}\\right)^{2}\\right)$"
+      ],
+      "text/plain": [
+       "2*Z1*Z2*(Y1*Z2 + Y2*Z1)*(Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 + (X1*Z2 + X2*Z1)**2 - (X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 + 2*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2)"
+      ]
+     },
+     "execution_count": 20,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "res.X.x"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "id": "89a34db8-81d2-4741-ae0a-5268b46d0974",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle - 2 Z_{1}^{2} Z_{2}^{2} \\left(Y_{1} Z_{2} + Y_{2} Z_{1}\\right)^{4} + \\left(- X_{1} X_{2} Z_{1} Z_{2} + Z_{1}^{2} Z_{2}^{2} a + \\left(X_{1} Z_{2} + X_{2} Z_{1}\\right)^{2}\\right) \\left(- 3 Z_{1}^{2} Z_{2}^{2} \\left(Y_{1} Z_{2} + Y_{2} Z_{1}\\right)^{4} - 3 \\left(X_{1} Z_{2} + X_{2} Z_{1}\\right)^{2} + 3 \\left(X_{1} Z_{2} + X_{2} Z_{1} + Z_{1} Z_{2} \\left(Y_{1} Z_{2} + Y_{2} Z_{1}\\right)^{2}\\right)^{2} - 4 \\left(- X_{1} X_{2} Z_{1} Z_{2} + Z_{1}^{2} Z_{2}^{2} a + \\left(X_{1} Z_{2} + X_{2} Z_{1}\\right)^{2}\\right)^{2}\\right)$"
+      ],
+      "text/plain": [
+       "-2*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 + (-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)*(-3*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 - 3*(X1*Z2 + X2*Z1)**2 + 3*(X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 - 4*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2)"
+      ]
+     },
+     "execution_count": 21,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "res.Y.x"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "id": "8d3db4ac-95fb-4483-a572-265744461daf",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle 4 Z_{1}^{3} Z_{2}^{3} \\left(Y_{1} Z_{2} + Y_{2} Z_{1}\\right)^{3}$"
+      ],
+      "text/plain": [
+       "4*Z1**3*Z2**3*(Y1*Z2 + Y2*Z1)**3"
+      ]
+     },
+     "execution_count": 22,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "res.Z.x"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "id": "645ca07a-8d47-446d-93ca-eab4de26dfff",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "U1 = X1*Z2 \t X1*Z2\n",
+      "U2 = X2*Z1 \t X2*Z1\n",
+      "S1 = Y1*Z2 \t Y1*Z2\n",
+      "S2 = Y2*Z1 \t Y2*Z1\n",
+      "ZZ = Z1*Z2 \t Z1*Z2\n",
+      "T = U1+U2 \t X1*Z2 + X2*Z1\n",
+      "TT = T \t (X1*Z2 + X2*Z1)**2\n",
+      "M = S1+S2 \t Y1*Z2 + Y2*Z1\n",
+      "t0 = ZZ \t Z1**2*Z2**2\n",
+      "t1 = a*t0 \t Z1**2*Z2**2*a\n",
+      "t2 = U1*U2 \t X1*X2*Z1*Z2\n",
+      "t3 = TT-t2 \t -X1*X2*Z1*Z2 + (X1*Z2 + X2*Z1)**2\n",
+      "R = t3+t1 \t -X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2\n",
+      "F = ZZ*M \t Z1*Z2*(Y1*Z2 + Y2*Z1)\n",
+      "L = M*F \t Z1*Z2*(Y1*Z2 + Y2*Z1)**2\n",
+      "LL = L \t Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4\n",
+      "t4 = T+L \t X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2\n",
+      "t5 = t4 \t (X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2\n",
+      "t6 = t5-TT \t -(X1*Z2 + X2*Z1)**2 + (X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2\n",
+      "G = t6-LL \t -Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 - (X1*Z2 + X2*Z1)**2 + (X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2\n",
+      "t7 = R \t (-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2\n",
+      "t8 = 2*t7 \t 2*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2\n",
+      "W = t8-G \t Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 + (X1*Z2 + X2*Z1)**2 - (X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 + 2*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2\n",
+      "t9 = F*W \t Z1*Z2*(Y1*Z2 + Y2*Z1)*(Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 + (X1*Z2 + X2*Z1)**2 - (X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 + 2*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2)\n",
+      "X3 = 2*t9 \t 2*Z1*Z2*(Y1*Z2 + Y2*Z1)*(Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 + (X1*Z2 + X2*Z1)**2 - (X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 + 2*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2)\n",
+      "t10 = 2*W \t 2*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 + 2*(X1*Z2 + X2*Z1)**2 - 2*(X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 + 4*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2\n",
+      "t11 = G-t10 \t -3*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 - 3*(X1*Z2 + X2*Z1)**2 + 3*(X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 - 4*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2\n",
+      "t12 = 2*LL \t 2*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4\n",
+      "t13 = R*t11 \t (-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)*(-3*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 - 3*(X1*Z2 + X2*Z1)**2 + 3*(X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 - 4*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2)\n",
+      "Y3 = t13-t12 \t -2*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 + (-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)*(-3*Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**4 - 3*(X1*Z2 + X2*Z1)**2 + 3*(X1*Z2 + X2*Z1 + Z1*Z2*(Y1*Z2 + Y2*Z1)**2)**2 - 4*(-X1*X2*Z1*Z2 + Z1**2*Z2**2*a + (X1*Z2 + X2*Z1)**2)**2)\n",
+      "t14 = F \t Z1**2*Z2**2*(Y1*Z2 + Y2*Z1)**2\n",
+      "t15 = F*t14 \t Z1**3*Z2**3*(Y1*Z2 + Y2*Z1)**3\n",
+      "Z3 = 4*t15 \t 4*Z1**3*Z2**3*(Y1*Z2 + Y2*Z1)**3\n"
+     ]
+    }
+   ],
+   "source": [
+    "for op_result in ctx.actions.get_by_index([0])[0].op_results:\n",
+    "    print(repr(op_result), \"\\t\", op_result.value.x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "08e17876-331a-445a-8fbc-1990a0b5f69d",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}