Add more detailed RPA-RE exploration.

1 files changed, 173 insertions, 18 deletions

diff --git a/re/rpa.ipynb b/re/rpa.ipynb
index 995e62d..1dfe668 100644
--- a/re/rpa.ipynb
+++ b/re/rpa.ipynb
@@ -18,11 +18,14 @@
     "from collections import Counter\n",
     "import numpy as np\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",
+    "from concurrent.futures import ProcessPoolExecutor, as_completed\n",
     "\n",
     "from IPython.display import HTML, display\n",
-    "from tqdm.auto import tqdm\n",
+    "from tqdm.auto import tqdm, trange\n",
     "import tabulate\n",
     "\n",
     "from pyecsca.ec.model import ShortWeierstrassModel\n",
@@ -52,6 +55,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "%matplotlib inline\n",
     "hv.extension(\"bokeh\")"
    ]
   },
@@ -113,14 +117,14 @@
     "    BinaryNAFMultiplier(add, dbl, neg, None, ProcessingDirection.LTR, 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, 4, None, AccumulationOrder.PeqPR, False, True),\n",
+    "    #WindowNAFMultiplier(add, dbl, neg, 4, None, AccumulationOrder.PeqPR, False, True),\n",
     "    SlidingWindowMultiplier(add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
     "    SlidingWindowMultiplier(add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
     "    FixedWindowLTRMultiplier(add, dbl, 4, None, AccumulationOrder.PeqPR, True),\n",
     "    FixedWindowLTRMultiplier(add, dbl, 5, 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",
-    "    FullPrecompMultiplier(add, dbl, None, False, ProcessingDirection.RTL, AccumulationOrder.PeqPR, True, True),\n",
+    "    #FullPrecompMultiplier(add, dbl, None, False, ProcessingDirection.RTL, AccumulationOrder.PeqPR, True, True),\n",
     "    BGMWMultiplier(add, dbl, 3, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
     "    BGMWMultiplier(add, dbl, 5, None, ProcessingDirection.LTR, AccumulationOrder.PeqPR, True),\n",
     "    CombMultiplier(add, dbl, 3, None, AccumulationOrder.PeqPR, True),\n",
@@ -140,7 +144,9 @@
    "cell_type": "code",
    "execution_count": null,
    "id": "44e7e2e9-e0ad-4c8d-8605-af68f73d73e2",
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "scalar = 0b1000000000000000000000000000000000000000000000000\n",
@@ -148,6 +154,7 @@
     "scalar = 0b1010101010101010101010101010101010101010101010101\n",
     "scalar = 0b1111111111111111111111110000000000000000000000000\n",
     "scalar = 123456789123456789\n",
+    "scarar = 8750920244948492046\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",
@@ -205,7 +212,9 @@
    "cell_type": "code",
    "execution_count": null,
    "id": "8113cb3f-dc06-4cb7-955c-11cedb4fbdd7",
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "table = [[\"Multiplier\", \"zero present\", \"multiple computed\"]]\n",
@@ -214,9 +223,11 @@
     "    with local(MultipleContext()) as ctx:\n",
     "        mult.init(params, P0_target)\n",
     "        res = mult.multiply(scalar)\n",
-    "    table.append([str(mult), any(map(lambda P: P.X == 0, ctx.points.keys())), k in ctx.points.values()])\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=\"html\", headers=\"firstrow\")))"
+    "display(HTML(tabulate.tabulate(table, tablefmt=\"unsafehtml\", headers=\"firstrow\", colalign=(\"left\", \"center\", \"center\"))))"
    ]
   },
   {
@@ -258,8 +269,10 @@
    "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)."
+    "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) and its error (flip) probability."
    ]
   },
   {
@@ -269,31 +282,171 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def simulated_oracle(scalar, affine_point, simulate_mult_id = 0):\n",
+    "def simulated_oracle(scalar, affine_point, simulate_mult_id=0, flip_proba=0):\n",
     "    real_mult = multipliers[simulate_mult_id]\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()))"
+    "    real_result = any(map(lambda P: P.X == 0 or P.Y == 0, ctx.points.keys()))\n",
+    "    change = bernoulli(flip_proba).rvs()\n",
+    "    return bool(real_result ^ change)"
+   ]
+  },
+  {
+   "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",
+   "execution_count": null,
+   "id": "b6c70d89-1c7d-4d7c-bc65-0cf766b86c0a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res = rpa_distinguish(params, multipliers, simulated_oracle)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "068e1ba5-9884-4d2b-97f6-d54313daddad",
+   "metadata": {},
+   "source": [
+    "Now we can examine how the method performs in the presence of noise and with various majority vote parameters. Note that the code below spawns several processes (`num_cores`) and saturates their CPU fully, so set this to something appropriate."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "265eb13a-6028-4fde-9c3f-cc23768ba63e",
-   "metadata": {
-    "scrolled": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "table = [[\"True multiplier\", \"Reversed\", \"Correct\", \"Remaining\"]]\n",
+    "errs = (0, 0.1, 0.2, 0.3, 0.4, 0.5)\n",
+    "majs = (1, 3, 5, 7)\n",
+    "num_tries = 100\n",
+    "num_cores = 4\n",
+    "\n",
+    "def measure_mult(params, multipliers, simulated_oracle, i, mult, err, majority):\n",
+    "    correct = 0\n",
+    "    precise = 0\n",
+    "    calls = 0\n",
+    "    p = partial(simulated_oracle, simulate_mult_id=i, flip_proba=err)\n",
+    "    def oracle(scalar, affine_point):\n",
+    "        nonlocal calls\n",
+    "        calls += 1\n",
+    "        return p(scalar, affine_point)\n",
+    "    for j in range(num_tries):\n",
+    "        res = rpa_distinguish(params, multipliers, oracle, majority=majority)\n",
+    "        if mult in res:\n",
+    "            correct += 1\n",
+    "            if len(res) == 1:\n",
+    "                precise += 1\n",
+    "    return correct, precise, calls\n",
+    "\n",
     "with TemporaryConfig() as cfg:\n",
     "    cfg.log.enabled = False\n",
-    "    for i, mult in tqdm(enumerate(multipliers)):\n",
-    "        res = rpa_distinguish(params, multipliers, partial(simulated_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\")))"
+    "    with ProcessPoolExecutor(max_workers=num_cores) as pool:\n",
+    "        futures = []\n",
+    "        args = []\n",
+    "        for i, mult in enumerate(multipliers):\n",
+    "            for err in errs:\n",
+    "                for majority in majs:\n",
+    "                    a = (params, multipliers, simulated_oracle, i, mult, err, majority)\n",
+    "                    futures.append(pool.submit(measure_mult, *a))\n",
+    "                    args.append(a)\n",
+    "        results = [None for _ in futures]\n",
+    "        for future in tqdm(as_completed(futures), total=len(futures)):\n",
+    "            j = futures.index(future)\n",
+    "            a = args[j]\n",
+    "            results[j] = future.result()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4d283a04-decc-422f-b804-addd57a8a635",
+   "metadata": {},
+   "source": [
+    "Now we accumulate the results across the error rate and majority vote parameters and plot two heatmaps:\n",
+    " - One for the average number of queries to the oracle.\n",
+    " - One for the success rate of the distinguisher."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6cc0af93-8442-4ccf-958f-3502d2d76a34",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "correct_rate = np.zeros((len(errs), len(majs)))\n",
+    "precise_rate = np.zeros((len(errs), len(majs)))\n",
+    "query_rate = np.zeros((len(errs), len(majs)))\n",
+    "for a, result in zip(args, results):\n",
+    "    i = errs.index(a[5])\n",
+    "    j = len(majs) - majs.index(a[6]) - 1\n",
+    "    correct_rate[i, j] += result[0]\n",
+    "    precise_rate[i, j] += result[1]\n",
+    "    query_rate[i, j] += result[2]\n",
+    "\n",
+    "correct_rate = (correct_rate * 100) / (num_tries * len(multipliers))\n",
+    "precise_rate = (precise_rate * 100) / (num_tries * len(multipliers))\n",
+    "query_rate = query_rate / (num_tries * len(multipliers))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3caf3afd-e27e-4b6f-947b-cd8edb2d568d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "im = ax.imshow(query_rate.T, vmin=0, cmap=\"plasma\")\n",
+    "cbar_ax = fig.add_axes((0.97, 0.15, 0.04, 0.69))\n",
+    "cbar = fig.colorbar(im, cax=cbar_ax)\n",
+    "cbar.ax.set_ylabel(\"Average oracle query rate\", rotation=-90, va=\"bottom\")\n",
+    "\n",
+    "ax.set_xticks(np.arange(len(errs)), labels=errs)\n",
+    "ax.set_yticks(np.arange(len(majs)), labels=reversed(majs))\n",
+    "ax.set_xlabel(\"err\")\n",
+    "ax.set_ylabel(\"majority\")\n",
+    "for i in range(len(errs)):\n",
+    "    for j in range(len(majs)):\n",
+    "        text = ax.text(i, j, f\"{query_rate[i, j]:.2f}\",\n",
+    "                       ha=\"center\", va=\"center\", color=\"w\")\n",
+    "fig.savefig(\"rpa_re_query_rate.png\", bbox_inches=\"tight\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "09c976dd-4940-40be-b7b2-6da73d672f24",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "im = ax.imshow(correct_rate.T, vmin=0, cmap=\"viridis\")\n",
+    "cbar_ax = fig.add_axes((0.97, 0.15, 0.04, 0.69))\n",
+    "cbar = fig.colorbar(im, cax=cbar_ax)\n",
+    "cbar.ax.set_ylabel(\"Success rate\", rotation=-90, va=\"bottom\")\n",
+    "cbar.ax.axhline(100 / len(multipliers), color=\"red\", linestyle=\"--\")\n",
+    "\n",
+    "ax.set_xticks(np.arange(len(errs)), labels=errs)\n",
+    "ax.set_yticks(np.arange(len(majs)), labels=reversed(majs))\n",
+    "ax.set_xlabel(\"err\")\n",
+    "ax.set_ylabel(\"majority\")\n",
+    "for i in range(len(errs)):\n",
+    "    for j in range(len(majs)):\n",
+    "        text = ax.text(i, j, f\"{correct_rate[i, j]:.2f}%\",\n",
+    "                       ha=\"center\", va=\"center\", color=\"w\")\n",
+    "fig.savefig(\"rpa_re_success_rate.png\", bbox_inches=\"tight\")\n",
+    "plt.show()"
    ]
   },
   {
@@ -301,6 +454,8 @@
    "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",