path: root/re/rpa.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "66c30004-cd2c-4d34-9999-f33f9e6fd5e9",
   "metadata": {},
   "source": [
    "# RPA-based reverse-engineering\n",
    "This notebook showcases the RPA-based reverse-engineering technique for scalar multipliers.\n",
    "\n",
    " - [Exploration](#Exploration)\n",
    " - [Reverse-engineering](#Reverse-engineering)\n",
    "   - [Oracle simulation](#Oracle-simulation)\n",
    "     - [Symmetric noise](#What-about-(symmetric)-noise?)\n",
    "     - [Asymmetric noise](#What-about-(asymmetric)-noise?)\n",
    "   - [Method simulation](#Method-simulation)"
   ]
  },
  {
   "cell_type": "code",
   "id": "11a5f41d-1471-49c3-ba7c-ac87470a31d0",
   "metadata": {},
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import holoviews as hv\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.signal import find_peaks\n",
    "from functools import partial\n",
    "from scipy.stats import bernoulli\n",
    "\n",
    "from IPython.display import HTML, display\n",
    "from tqdm.auto import tqdm\n",
    "import tabulate\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, get_params\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.misc.utils import silent, TaskExecutor\n",
    "from pyecsca.sca.trace.process import normalize\n",
    "from pyecsca.sca.trace.combine import average, subtract\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, RPA\n",
    "from pyecsca.sca.trace import Trace\n",
    "\n",
    "from eval import (eval_tree_symmetric, eval_tree_asymmetric,\n",
    "                    success_rate_symmetric, success_rate_asymmetric,\n",
    "                    query_rate_symmetric, query_rate_asymmetric,\n",
    "                    precise_rate_symmetric, precise_rate_asymmetric,\n",
    "                    amount_rate_symmetric, amount_rate_asymmetric,\n",
    "                    success_rate_vs_majority_symmetric, success_rate_vs_majority_asymmetric,\n",
    "                    success_rate_vs_query_rate_symmetric, load, store)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "3ff1d591-f922-4c94-9e47-ab053fc21cf1",
   "metadata": {},
   "source": [
    "%matplotlib ipympl\n",
    "hv.extension(\"bokeh\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "194fff59-1c4b-473a-9ffc-99b256aecc24",
   "metadata": {},
   "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)\n",
    "\n",
    "# And P-256 for eval\n",
    "p256 = get_params(\"secg\", \"secp256r1\", \"projective\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "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",
   "id": "febb198f-4370-4abd-8edc-17c5c1da8d0f",
   "metadata": {},
   "source": [
    "multipliers = [\n",
    "    LTRMultiplier(add, dbl, None, False, AccumulationOrder.PeqPR, True, True),\n",
    "    LTRMultiplier(add, dbl, None, True, AccumulationOrder.PeqPR, True, True),\n",
    "    RTLMultiplier(add, dbl, None, False, AccumulationOrder.PeqPR, True),\n",
    "    RTLMultiplier(add, dbl, None, True, AccumulationOrder.PeqPR, False),\n",
    "    SimpleLadderMultiplier(add, dbl, None, True, True),\n",
    "    BinaryNAFMultiplier(add, dbl, neg, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    BinaryNAFMultiplier(add, dbl, neg, None, ProcessingDirection.RTL, AccumulationOrder.PeqPR, True),\n",
    "    WindowNAFMultiplier(add, dbl, neg, 3, None, AccumulationOrder.PeqPR, True, True),\n",
    "    WindowNAFMultiplier(add, dbl, neg, 4, None, AccumulationOrder.PeqPR, True, True),\n",
    "    WindowNAFMultiplier(add, dbl, neg, 5, None, AccumulationOrder.PeqPR, True, True),\n",
    "    WindowBoothMultiplier(add, dbl, neg, 3, None, AccumulationOrder.PeqPR, True, True),\n",
    "    WindowBoothMultiplier(add, dbl, neg, 4, None, AccumulationOrder.PeqPR, True, True),\n",
    "    WindowBoothMultiplier(add, dbl, neg, 5, None, AccumulationOrder.PeqPR, True, True),\n",
    "    SlidingWindowMultiplier(add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    SlidingWindowMultiplier(add, dbl, 4, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    SlidingWindowMultiplier(add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    SlidingWindowMultiplier(add, dbl, 3, None, ProcessingDirection.RTL, AccumulationOrder.PeqPR, True),\n",
    "    SlidingWindowMultiplier(add, dbl, 4, None, ProcessingDirection.RTL, AccumulationOrder.PeqPR, True),\n",
    "    SlidingWindowMultiplier(add, dbl, 5, None, ProcessingDirection.RTL, AccumulationOrder.PeqPR, True),\n",
    "    FixedWindowLTRMultiplier(add, dbl, 3, None, AccumulationOrder.PeqPR, True),\n",
    "    FixedWindowLTRMultiplier(add, dbl, 4, None, AccumulationOrder.PeqPR, True),\n",
    "    FixedWindowLTRMultiplier(add, dbl, 5, None, AccumulationOrder.PeqPR, True),\n",
    "    FixedWindowLTRMultiplier(add, dbl, 8, None, AccumulationOrder.PeqPR, True),\n",
    "    FixedWindowLTRMultiplier(add, dbl, 16, None, AccumulationOrder.PeqPR, True),\n",
    "    FullPrecompMultiplier(add, dbl, None, True, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True, True),\n",
    "    FullPrecompMultiplier(add, dbl, None, False, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True, True),\n",
    "    BGMWMultiplier(add, dbl, 2, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    BGMWMultiplier(add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    BGMWMultiplier(add, dbl, 4, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    BGMWMultiplier(add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
    "    CombMultiplier(add, dbl, 2, None, AccumulationOrder.PeqPR, True),\n",
    "    CombMultiplier(add, dbl, 3, None, AccumulationOrder.PeqPR, True),\n",
    "    CombMultiplier(add, dbl, 4, None, AccumulationOrder.PeqPR, True),\n",
    "    CombMultiplier(add, dbl, 5, None, AccumulationOrder.PeqPR, True)\n",
    "]\n",
    "print(len(multipliers))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "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",
   "id": "44e7e2e9-e0ad-4c8d-8605-af68f73d73e2",
   "metadata": {},
   "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]P\n",
    "multiples = {}\n",
    "\n",
    "table = [[\"Multiplier\", \"multiples\"]]\n",
    "\n",
    "for mult in multipliers:\n",
    "    with local(MultipleContext()) as ctx:\n",
    "        mult.init(params, g)\n",
    "        res = mult.multiply(scalar)\n",
    "    for m in ctx.points.values():\n",
    "        s = multiples.setdefault(m, set())\n",
    "        s.add(mult)\n",
    "    table.append([str(mult), str(list(ctx.points.values()))])\n",
    "\n",
    "display(HTML(tabulate.tabulate(table, tablefmt=\"html\", headers=\"firstrow\")))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "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",
   "id": "a7fb8a3f-7938-493b-88dc-582ba4d8959d",
   "metadata": {},
   "source": [
    "k = 108\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))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "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",
   "id": "8113cb3f-dc06-4cb7-955c-11cedb4fbdd7",
   "metadata": {},
   "source": [
    "table = [[\"Multiplier\", \"zero present\", \"multiple computed\"]]\n",
    "\n",
    "for mult in multipliers:\n",
    "    with local(MultipleContext()) as ctx:\n",
    "        mult.init(params, P0_target)\n",
    "        res = mult.multiply(scalar)\n",
    "    zero = any(map(lambda P: P.X == 0 or P.Y == 0, ctx.points.keys()))\n",
    "    multiple = k in ctx.points.values()\n",
    "    table.append([str(mult), f\"<b>{zero}</b>\" if zero else zero, f\"<b>{multiple}</b>\" if multiple else multiple])\n",
    "\n",
    "display(HTML(tabulate.tabulate(table, tablefmt=\"unsafehtml\", headers=\"firstrow\", colalign=(\"left\", \"center\", \"center\"))))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "0d2b9fe8-5064-4887-b782-dcfe9f42d217",
   "metadata": {},
   "source": [
    "Now lets look at the relation of multiples to multipliers."
   ]
  },
  {
   "cell_type": "code",
   "id": "67d7705c-6a41-47d9-ad1e-23ea549aaf00",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "table = [[\"Multiple\", \"Multipliers\"]]\n",
    "for multiple, mults in multiples.items():\n",
    "    table.append([bin(multiple), [mult.__class__.__name__ for mult in mults]])\n",
    "\n",
    "display(HTML(tabulate.tabulate(table, tablefmt=\"html\", headers=\"firstrow\")))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "9b8e4338-a2a8-468e-8873-c18c77260cfc",
   "metadata": {},
   "source": [
    "Note that all of the exploration so far was in a context of a fixed scalar. Even though for a given scalar some multipliers might be indistinguishable from the perspective of the multiples they compute, there may be other scalars that distinguish them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9015f3c-fba6-4614-b722-848b8522d072",
   "metadata": {},
   "source": [
    "## Reverse-engineering\n",
    "\n",
    "### Oracle simulation\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. The oracle is parametrized with the simulated multiplier index in the table of multipliers (it simulates this \"real\" multiplier). Furthermore, lets also examine a `noisy_oracle` (with a flip probability) and a `biased_oracle` (with asymmetric flip probability).\n",
    "\n",
    "Note that the oracle has two additional parameters `measure_init` and `measure_multiply` which determine whether the oracle considers the zero coordinate point in scalar multiplier initialization (precomputation) and in scalar multiplier multiplication, respectively. This is important for scalar multipliers with precomputation as there one might be able to separate the precomputation and multiplication stages and obtain oracle answers on both separately."
   ]
  },
  {
   "cell_type": "code",
   "id": "9bb61ac5-d837-4287-a5de-a9a63c346acf",
   "metadata": {},
   "source": [
    "def simulated_oracle(scalar, affine_point, simulate_mult_id=0, measure_init=True, measure_multiply=True, randomize=False):\n",
    "    real_mult = multipliers[simulate_mult_id]\n",
    "    point = affine_point.to_model(params.curve.coordinate_model, params.curve, randomized=randomize)\n",
    "    \n",
    "    # Simulate the multiplier init\n",
    "    with local(MultipleContext()) as ctx:\n",
    "        real_mult.init(params, point)\n",
    "    init_points = set(ctx.parents.keys())\n",
    "    init_parents = set(sum((ctx.parents[point] for point in init_points), []))\n",
    "    # Did zero happen in some input point during the init?\n",
    "    init_zero = any(map(lambda P: P.X == 0 or P.Y == 0, init_parents))\n",
    "    \n",
    "    # Simulate the multiplier multiply\n",
    "    with local(ctx) as ctx:\n",
    "        real_mult.multiply(scalar)\n",
    "    all_points = set(ctx.parents.keys())\n",
    "    multiply_parents = set(sum((ctx.parents[point] for point in all_points - init_points), []))\n",
    "    # Did zero happen in some input point during the multiply?\n",
    "    multiply_zero = any(map(lambda P: P.X == 0 or P.Y == 0, multiply_parents))\n",
    "    real_result = (init_zero and measure_init) or (multiply_zero and measure_multiply)\n",
    "    return real_result\n",
    "\n",
    "def noisy_oracle(oracle, flip_proba=0):\n",
    "    def noisy(*args, **kwargs):\n",
    "        real_result = oracle(*args, **kwargs)\n",
    "        change = bernoulli(flip_proba).rvs()\n",
    "        return bool(real_result ^ change)\n",
    "    return noisy\n",
    "\n",
    "def biased_oracle(oracle, flip_0=0, flip_1=0):\n",
    "    def biased(*args, **kwargs):\n",
    "        real_result = oracle(*args, **kwargs)\n",
    "        change = bernoulli(flip_1).rvs() if real_result else bernoulli(flip_0).rvs()\n",
    "        return bool(real_result ^ change)\n",
    "    return biased"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c1ad6235-52bf-4bab-90ff-55b67165390e",
   "metadata": {},
   "source": [
    "We can see how the RPA-RE method distinguishes a given multiplier:"
   ]
  },
  {
   "cell_type": "code",
   "id": "b6c70d89-1c7d-4d7c-bc65-0cf766b86c0a",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "res = rpa_distinguish(params, multipliers, simulated_oracle)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "6361c477-5ddf-46ff-8918-864a453b676b",
   "metadata": {},
   "source": [
    "Let's see if the result is correct. You can replace the `simulated_oracle` above with `noisy_oracle(simulated_oracle, flip_proba=0.2)` or with `biased_oracle(simulated_oracle, flip_0=0.2, flip_1=0.1)` to see how the process and result changes with noise."
   ]
  },
  {
   "cell_type": "code",
   "id": "1556e604-d0ab-403c-bd44-d53be8e18283",
   "metadata": {},
   "source": [
    "print(multipliers[0] in res)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "4cecc7dc-f609-4073-b0db-ac9631ac3edf",
   "metadata": {},
   "source": [
    "We can also have a look at the distinguishing tree that the method builds for this set of multipliers."
   ]
  },
  {
   "cell_type": "code",
   "id": "838ae83b-771d-4e4c-8977-ba997aa4fbb6",
   "metadata": {},
   "source": [
    "re = RPA(set(multipliers))\n",
    "with silent():\n",
    "    re.build_tree(p256, tries=10)\n",
    "print(re.tree.describe())"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "9a99df60-f6ec-40e1-a43a-d95eb64beb8c",
   "metadata": {},
   "source": [
    "We can also look at the rough tree structure."
   ]
  },
  {
   "cell_type": "code",
   "id": "136000fe-4a4a-4261-bacc-a49102080749",
   "metadata": {},
   "source": [
    "print(re.tree.render_basic())"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "068e1ba5-9884-4d2b-97f6-d54313daddad",
   "metadata": {},
   "source": [
    "#### What about (symmetric) noise?\n",
    "Now we can examine how the method performs in the presence of noise and with various majority vote parameters. 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.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "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.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "num_workers = 30",
   "id": "74b83f35d681ab9b",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "04bd8191-1af8-4183-ac95-31b1e1c48e1b",
   "metadata": {},
   "source": [
    "def build_tree(cfgs):\n",
    "    with silent():\n",
    "        re = RPA(set(cfgs))\n",
    "        re.build_tree(p256, tries=10)\n",
    "        return re.tree\n",
    "\n",
    "correct_rate, precise_rate, amount_rate, query_rate = eval_tree_symmetric(set(multipliers), build_tree, num_trees=100, num_tries=100, num_cores=num_workers)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "9bda1baa-359a-4f9b-889d-f64e055deff6",
   "metadata": {},
   "source": [
    "We can plot several heatmaps:\n",
    " - One for the average number of queries to the oracle.\n",
    " - One for the success rate of the reverse-engineering.\n",
    " - One for the precision of the reverse-engineering."
   ]
  },
  {
   "cell_type": "code",
   "id": "7be15799-d042-43ae-b052-ef7b103e1cca",
   "metadata": {},
   "source": [
    "success_rate_symmetric(correct_rate, 100 / len(multipliers)).savefig(\"rpa_re_success_rate_symmetric.pdf\", bbox_inches=\"tight\")\n",
    "query_rate_symmetric(query_rate).savefig(\"rpa_re_query_rate_symmetric.pdf\", bbox_inches=\"tight\")\n",
    "precise_rate_symmetric(precise_rate).savefig(\"rpa_re_precise_rate_symmetric.pdf\", bbox_inches=\"tight\")\n",
    "amount_rate_symmetric(amount_rate).savefig(\"rpa_re_amount_rate_symmetric.pdf\", bbox_inches=\"tight\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b8eeb49e-7734-4fda-bf28-ac39f6c8a626",
   "metadata": {},
   "source": [
    "Another way to look at these metrics is a scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "id": "d5176863-d6b6-4205-8a3a-4453b8177305",
   "metadata": {},
   "source": [
    "success_rate_vs_query_rate_symmetric(query_rate, correct_rate).savefig(\"rpa_re_scatter_symmetric.pdf\", bbox_inches=\"tight\")\n",
    "success_rate_vs_majority_symmetric(correct_rate).savefig(\"rpa_re_plot_symmetric.pdf\", bbox_inches=\"tight\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "8c747434-84bb-4acd-a994-45a4b4859a6e",
   "metadata": {},
   "source": [
    "And save the results for later."
   ]
  },
  {
   "cell_type": "code",
   "id": "857f39ad-f6ba-4006-8da9-92f4318819e2",
   "metadata": {},
   "source": [
    "store(\"rpa_re_symmetric.nc\", correct_rate, precise_rate, amount_rate, query_rate)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "efb26038-87ea-4a5d-8682-f855c2bec13f",
   "metadata": {},
   "source": [
    "correct_rate, precise_rate, amount_rate, query_rate = load(\"rpa_re_symmetric.nc\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "e3981f44-ed3b-43f5-b4bd-e2d4b2ff95e8",
   "metadata": {},
   "source": [
    "#### What about (asymmetric) noise?\n",
    "The oracle may not only be noisy, but biased, this computation evaluates that case. Beware, for the same parameters this is about 6x slower because of the other dimension (two error probabilities instead of one).\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "This is a resource intensive cell that uses parallelism. Recall the `num_workers` variable.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "id": "06268c58-63c9-4565-ba52-4414a9939581",
   "metadata": {},
   "source": [
    "def build_tree(cfgs):\n",
    "    with silent():\n",
    "        re = RPA(set(cfgs))\n",
    "        re.build_tree(p256, tries=10)\n",
    "        return re.tree\n",
    "\n",
    "correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b = eval_tree_asymmetric(set(multipliers), build_tree, num_trees=100, num_tries=100, num_cores=num_workers)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "43a0055a-8d62-4f0a-be9f-e6ccb47b1d11",
   "metadata": {},
   "source": [
    "success_rate_asymmetric(correct_rate_b, 100 / len(multipliers)).savefig(\"rpa_re_success_rate_asymmetric.pdf\", bbox_inches=\"tight\")\n",
    "query_rate_asymmetric(query_rate_b).savefig(\"rpa_re_query_rate_asymmetric.pdf\", bbox_inches=\"tight\")\n",
    "precise_rate_asymmetric(precise_rate_b).savefig(\"rpa_re_precise_rate_asymmetric.pdf\", bbox_inches=\"tight\")\n",
    "amount_rate_asymmetric(amount_rate_b).savefig(\"rpa_re_amount_rate_asymmetric.pdf\", bbox_inches=\"tight\")\n",
    "success_rate_vs_majority_asymmetric(correct_rate_b).savefig(\"rpa_re_plot_asymmetric.pdf\", bbox_inches=\"tight\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "67a6bba8-84e8-48f1-9b3b-b0b8715c71a6",
   "metadata": {},
   "source": [
    "And save the results for later."
   ]
  },
  {
   "cell_type": "code",
   "id": "85933a5e-f526-4002-9203-0e4ae4f731d0",
   "metadata": {},
   "source": [
    "store(\"rpa_re_asymmetric.nc\", correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "ddde34f8-7cde-449d-a2e3-289e1aefba72",
   "metadata": {},
   "source": [
    "correct_rate_b, precise_rate_b, amount_rate_b, query_rate_b = load(\"rpa_re_asymmetric.nc\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "62b0b8f2-8149-4abd-9aa7-a056b237ac6e",
   "metadata": {},
   "source": [
    "### Method simulation\n",
    "\n",
    "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. The oracle is parametrized with the simulated multiplier index in the table of multipliers (it simulates this \"real\" multiplier)."
   ]
  },
  {
   "cell_type": "code",
   "id": "96bec03e-5397-440b-9e8c-81ba5253921b",
   "metadata": {
    "scrolled": true
   },
   "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[0].walk(callback)\n",
    "    return Trace(np.array(trace))\n",
    "\n",
    "def simulated_rpa_trace(mult, scalar, affine_point, noise, num_target=10, num_random=10):\n",
    "    random_traces = [noise(normalize(simulate_trace(mult, scalar, params.curve.affine_random().to_model(params.curve.coordinate_model, params.curve, randomized=True)))) for _ in range(num_random)]\n",
    "    target_traces = [noise(normalize(simulate_trace(mult, scalar, affine_point.to_model(params.curve.coordinate_model, params.curve, randomized=True)))) for _ in range(num_target)]\n",
    "\n",
    "    random_avg = average(*random_traces)\n",
    "    target_avg = average(*target_traces)\n",
    "\n",
    "    diff_trace = subtract(random_avg, target_avg)\n",
    "    return diff_trace\n",
    "\n",
    "def simulated_rpa_oracle(scalar, affine_point, simulate_mult_id = 0, variance=1):\n",
    "    real_mult = multipliers[simulate_mult_id]\n",
    "    noise = NormalNoice(0, variance)\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"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "d74a944f-e0a6-434d-8a12-138dfb9d516e",
   "metadata": {},
   "source": [
    "table = [[\"True multiplier\", \"Reversed\", \"Correct\", \"Remaining\"]]\n",
    "with silent():\n",
    "    for i, mult in tqdm(enumerate(multipliers), total=len(multipliers)):\n",
    "        res = rpa_distinguish(params, multipliers, partial(simulated_rpa_oracle, simulate_mult_id = i))\n",
    "        table.append([mult, res, mult in res, len(res)])\n",
    "display(HTML(tabulate.tabulate(table, tablefmt=\"html\", headers=\"firstrow\")))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "e694b3b3-290d-4528-a611-16a183662944",
   "metadata": {},
   "source": [
    "Note that the oracle function above has several parameters, like noise standard deviation, amount of traces simulated, and peak finding height threshold. Below we analyze how these parameters influence the resulting error probabilities. *If you changed any of the `scalar`, `k`, `P0`, `P0_target` or `multipliers` variables above, the cell below may not work correctly.*\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "This is a resource intensive cell that uses parallelism. Recall the `num_workers` variable.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "id": "e47a70d7-7993-4560-9257-dbbc1cf658cf",
   "metadata": {},
   "source": [
    "def eval(threshold, sdev, tries, num_traces):\n",
    "    # This requires that scalar and P0_target variables are consistent, i.e. that:\n",
    "    # - P0 is computed by scalarmult multipliers[0] when computing [scalar]P0_target\n",
    "    # - but P0 is not computed by scalarmult multipliers[7] when computing [scalar]P0_target\n",
    "    noise = NormalNoice(0, sdev)\n",
    "    aff = P0_target.to_affine()\n",
    "    true_pos = 0\n",
    "    false_pos = 0\n",
    "    for _ in range(tries):\n",
    "        diff_real = normalize(simulated_rpa_trace(multipliers[0], scalar, aff, noise, num_random=num_traces, num_target=num_traces))\n",
    "        true_pos += len(find_peaks(diff_real.samples, height=threshold)[0]) > 0\n",
    "        diff_nothing = normalize(simulated_rpa_trace(multipliers[7], scalar, aff, noise, num_random=num_traces, num_target=num_traces))\n",
    "        false_pos += len(find_peaks(diff_nothing.samples, height=threshold)[0]) > 0\n",
    "    false_neg = tries - true_pos\n",
    "    true_neg = tries - false_pos\n",
    "    return true_pos / tries, true_neg / tries, false_pos / tries, false_neg / tries\n",
    "\n",
    "threshold_range = [4]\n",
    "sdev_range = list(range(0, 11))\n",
    "traces_range = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]\n",
    "e0 = xr.DataArray(np.zeros((len(threshold_range), len(sdev_range), len(traces_range))),\n",
    "                  dims=(\"threshold\", \"sdev\", \"traces\"),\n",
    "                  coords={\"threshold\": threshold_range, \"sdev\": sdev_range, \"traces\": traces_range}, name=\"e0\")\n",
    "e1 = xr.DataArray(np.zeros((len(threshold_range), len(sdev_range), len(traces_range))),\n",
    "                  dims=(\"threshold\", \"sdev\", \"traces\"),\n",
    "                  coords={\"threshold\": threshold_range, \"sdev\": sdev_range, \"traces\": traces_range}, name=\"e1\")\n",
    "tries = 200\n",
    "with TaskExecutor(max_workers=num_workers) as pool:\n",
    "    for threshold in threshold_range:\n",
    "        for sdev in sdev_range:\n",
    "            for num_traces in traces_range:\n",
    "                pool.submit_task((threshold, sdev, num_traces),\n",
    "                                 eval, threshold, sdev, tries, num_traces)\n",
    "    for (threshold, sdev, num_traces), future in tqdm(pool.as_completed(), total=len(pool.tasks), smoothing=0):\n",
    "        true_pos, true_neg, false_pos, false_neg = future.result()\n",
    "        e0.loc[threshold, sdev, num_traces] = false_pos\n",
    "        e1.loc[threshold, sdev, num_traces] = false_neg"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "0b915148-4b0d-4c8d-a904-6b6ae85fafd8",
   "metadata": {},
   "source": [
    "fig, axs = plt.subplots(ncols=2, figsize=(10, 4), sharey=\"row\")\n",
    "for i, threshold in enumerate(threshold_range):\n",
    "    res0 = e0.sel(threshold=threshold).plot(ax=axs[0], vmin=0, vmax=1, cmap=\"plasma\", add_colorbar=False)\n",
    "    for j, sdev in enumerate(sdev_range):\n",
    "        for k, traces in enumerate(traces_range):\n",
    "            val = e0.sel(threshold=threshold, sdev=sdev, traces=traces)\n",
    "            sval = f\"{val:.2f}\"\n",
    "            color = \"white\" if val < 0.5 else \"black\"\n",
    "            if sval == \"0.00\":\n",
    "                color = \"grey\"\n",
    "            axs[0].text(traces, sdev, sval.lstrip(\"0\"), ha=\"center\", va=\"center\", color=color)\n",
    "    axs[0].set_title(\"$e_0$\")\n",
    "    axs[0].set_ylabel(\"noise $\\sigma$\")\n",
    "    axs[0].set_xlabel(\"traces per group\")\n",
    "    res1 = e1.sel(threshold=threshold).plot(ax=axs[1], vmin=0, vmax=1, cmap=\"plasma\", add_colorbar=False)\n",
    "    for j, sdev in enumerate(sdev_range):\n",
    "        for k, traces in enumerate(traces_range):\n",
    "            val = e1.sel(threshold=threshold, sdev=sdev, traces=traces)\n",
    "            sval = f\"{val:.2f}\"\n",
    "            color = \"white\" if val < 0.5 else \"black\"\n",
    "            if sval == \"0.00\":\n",
    "                color = \"grey\"\n",
    "            axs[1].text(traces, sdev, sval.lstrip(\"0\"), ha=\"center\", va=\"center\", color=color)\n",
    "    axs[1].set_title(\"$e_1$\")\n",
    "    axs[1].set_ylabel(\"noise $\\sigma$\")\n",
    "    axs[1].set_xlabel(\"traces per group\")\n",
    "    fig.tight_layout(h_pad=1.5, rect=(0, 0, 0.9, 1))\n",
    "    cbar_ymin, cbar_ymax = axs[0].get_position().ymin, axs[0].get_position().ymax\n",
    "    cbar_ax = fig.add_axes((0.92, 0.145, 0.02, 0.77))\n",
    "    cbar = fig.colorbar(res0, cax=cbar_ax, label=\"error probability\")\n",
    "    cbar.ax.yaxis.set_label_coords(2.8, 0.5)\n",
    "    cbar.ax.set_ylabel(\"error probability\", rotation=-90, va=\"bottom\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "989bcd51-888a-4a7a-bbfe-2eea4fee7adc",
   "metadata": {},
   "source": [
    "fig.savefig(\"rpa_re_errors.pdf\", bbox_inches=\"tight\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "1f809503-d572-4051-9731-a8a7d16ede33",
   "metadata": {},
   "source": [],
   "outputs": [],
   "execution_count": null
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}