8 files changed, 1906 insertions, 304 deletions

diff --git a/src/cz/crcs/ectester/common/util/ECUtil.java b/src/cz/crcs/ectester/common/util/ECUtil.java
index 9b0949d..d5e8006 100644
--- a/src/cz/crcs/ectester/common/util/ECUtil.java
+++ b/src/cz/crcs/ectester/common/util/ECUtil.java
@@ -362,11 +362,22 @@ public class ECUtil {
     public static BigInteger recoverSignatureNonce(byte[] signature, byte[] data, BigInteger privkey, ECParameterSpec params, String hashType) {
         try {
             int bitSize = params.getOrder().bitLength();
-            MessageDigest md = MessageDigest.getInstance(hashType);
-            byte[] hash = md.digest(data);
+            // Hash the data.
+            byte[] hash;
+            if (hashType.equals("NONE")) {
+                hash = data;
+            } else {
+                MessageDigest md = MessageDigest.getInstance(hashType);
+                hash = md.digest(data);
+            }
+            // Trim bitSize of rightmost bits.
             BigInteger hashInt = new BigInteger(1, hash);
-            hashInt = hashInt.and(BigInteger.ONE.shiftLeft(bitSize + 1).subtract(BigInteger.ONE));
+            int hashBits = hashInt.bitLength();
+            if (hashBits > bitSize) {
+                hashInt = hashInt.shiftRight(hashBits - bitSize);
+            }
 
+            // Parse DERSignature
             BigInteger[] sigPair = fromDERSignature(signature);
             BigInteger r = sigPair[0];
             BigInteger s = sigPair[1];
@@ -420,7 +431,7 @@ public class ECUtil {
             FileInputStream in = new FileInputStream(file);
             result.readCSV(in);
             in.close();
-        } else {
+        } else if (named != null) {
             if (params == EC_Consts.PARAMETER_W) {
                 result = EC_Store.getInstance().getObject(EC_Key.Public.class, named);
             } else if (params == EC_Consts.PARAMETER_S) {
diff --git a/src/cz/crcs/ectester/standalone/ECTesterStandalone.java b/src/cz/crcs/ectester/standalone/ECTesterStandalone.java
index 35db7da..1b4d9b8 100644
--- a/src/cz/crcs/ectester/standalone/ECTesterStandalone.java
+++ b/src/cz/crcs/ectester/standalone/ECTesterStandalone.java
@@ -369,7 +369,7 @@ public class ECTesterStandalone {
         ECPublicKey pubkey = (ECPublicKey) ECUtil.loadKey(EC_Consts.PARAMETER_W, cli.getOptionValue("ecdh.named-public"), cli.getOptionValue("ecdh.public"), (ECParameterSpec) spec);
 
         int amount = Integer.parseInt(cli.getOptionValue("ecdh.amount", "1"));
-        for (int i = 0; i < amount; ++i) {
+        for (int i = 0; i < amount || amount == 0; ++i) {
             if (!cli.hasOption("ecdh.fixed-private") && !cli.hasOption("ecdh.named-private") && !cli.hasOption("ecdh.private")) {
                 one = kpg.genKeyPair();
             }
@@ -506,7 +506,7 @@ public class ECTesterStandalone {
         ECPublicKey pubkey = (ECPublicKey) ECUtil.loadKey(EC_Consts.PARAMETER_W, cli.getOptionValue("ecdsa.named-public"), cli.getOptionValue("ecdsa.public"), spec);
 
         int amount = Integer.parseInt(cli.getOptionValue("ecdsa.amount", "1"));
-        for (int i = 0; i < amount; ++i) {
+        for (int i = 0; i < amount || amount == 0; ++i) {
             if (!cli.hasOption("ecdsa.named-private") || !cli.hasOption("ecdsa.named-public")) {
                 KeyPair one = kpg.genKeyPair();
 
diff --git a/util/plot_dh.ipynb b/util/plot_dh.ipynb
new file mode 100644
index 0000000..4d4edbc
--- /dev/null
+++ b/util/plot_dh.ipynb
@@ -0,0 +1,616 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Analysis of key generation data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:51:29.892989Z",
+     "start_time": "2019-03-17T19:51:29.557783Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    " %matplotlib notebook \n",
+    "import numpy as np\n",
+    "from scipy.stats import describe\n",
+    "from scipy.stats import norm as norm_dist\n",
+    "from scipy.stats.mstats import mquantiles\n",
+    "from math import log, sqrt\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import ticker, colors, gridspec\n",
+    "from copy import deepcopy\n",
+    "from utils import plot_hist, moving_average, hw\n",
+    "from binascii import unhexlify\n",
+    "from IPython.display import display, HTML\n",
+    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
+    "import ipywidgets as widgets\n",
+    "import tabulate"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Settings\n",
+    "Enter your input below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:57:52.012826Z",
+     "start_time": "2019-03-17T19:57:52.008374Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# File name with output from ECTesterReader or ECTesterStandalone ECDH.\n",
+    "fname = \"filename.csv\"\n",
+    "\n",
+    "# The amount of entries skipped from the beginning of the file, as they are usually outliers.\n",
+    "skip_first = 10\n",
+    "\n",
+    "# Whether to plot things in logarithmic scale or not.\n",
+    "log_scale = False\n",
+    "\n",
+    "# Whether to trim the time data outside the 1 - 99 percentile range (adjust below). Quite useful.\n",
+    "trim = True\n",
+    "\n",
+    "# How much to trim? Either a number in [0,1] signifying a quantile, or an absolute value signifying a threshold\n",
+    "trim_low = 0.01\n",
+    "trim_high = 0.99\n",
+    "\n",
+    "# Graphical (matplotlib) style name\n",
+    "style = \"ggplot\"\n",
+    "\n",
+    "# Color map to use, and what color to assign to \"bad\" values (necessary for log_scale)\n",
+    "color_map = plt.cm.plasma\n",
+    "color_map_bad = \"black\"\n",
+    "\n",
+    "# What function to use to calculate number of histogram bins of time\n",
+    "# one of \"sqrt\", \"sturges\", \"rice\", \"scott\" and \"fd\" or a number specifying the number of bins\n",
+    "hist_size = \"rice\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data processing"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:51:36.973070Z",
+     "start_time": "2019-03-17T19:51:36.967369Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Setup plot style\n",
+    "\n",
+    "plt.style.use(style)\n",
+    "\n",
+    "cmap = deepcopy(color_map)\n",
+    "cmap.set_bad(color_map_bad)\n",
+    "\n",
+    "# Normalization, linear or log.\n",
+    "if log_scale:\n",
+    "    norm = colors.LogNorm()\n",
+    "else:\n",
+    "    norm = colors.Normalize()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:51:39.208449Z",
+     "start_time": "2019-03-17T19:51:37.430702Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Read the header line.\n",
+    "\n",
+    "with open(fname, \"r\") as f:\n",
+    "    header = f.readline()\n",
+    "header_names = header.split(\";\")\n",
+    "if len(header_names) != 5:\n",
+    "    print(\"Bad data?\")\n",
+    "    exit(1)\n",
+    "\n",
+    "# Load the data\n",
+    "\n",
+    "hx = lambda x: int(x, 16)\n",
+    "data = np.genfromtxt(fname, delimiter=\";\", skip_header=1, converters={2: unhexlify, 3: hx, 4: hx},\n",
+    "                         dtype=np.dtype([(\"index\", \"u4\"), (\"time\", \"u4\"), (\"pub\", \"O\"), (\"priv\", \"O\"), (\"secret\", \"O\")]))\n",
+    "\n",
+    "time_unit = \"ms\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:57:56.363502Z",
+     "start_time": "2019-03-17T19:57:56.331005Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Setup the data\n",
+    "\n",
+    "# Skip first (outliers?)\n",
+    "\n",
+    "data = data[skip_first:]\n",
+    "\n",
+    "# If in nanoseconds, scale to microseconds\n",
+    "if header_names[1].endswith(\"[nano]\") and time_unit == \"ms\":\n",
+    "    time_unit = r\"$\\mu s$\"\n",
+    "    np.floor_divide(data[\"time\"], 1000, out=data[\"time\"])\n",
+    "\n",
+    "\n",
+    "# Trim times\n",
+    "quant_low_bound = trim_low if 0 <= trim_low <= 1 else 0.01\n",
+    "quant_high_bound = trim_high if 0 <= trim_high <= 1 else 0.95\n",
+    "quantiles = mquantiles(data[\"time\"], prob=(quant_low_bound, 0.25, 0.5, 0.75, quant_high_bound))\n",
+    "if trim:\n",
+    "    low_bound = quantiles[0] if 0 <= trim_low <= 1 else trim_low\n",
+    "    high_bound = quantiles[4] if 0 <= trim_high <= 1 else trim_high\n",
+    "    data_trimmed = data[np.logical_and(data[\"time\"] >= low_bound,\n",
+    "                                       data[\"time\"] <= high_bound)]\n",
+    "    quantiles_trim = mquantiles(data_trimmed[\"time\"], prob=(quant_low_bound, 0.25, 0.5, 0.75, quant_high_bound))\n",
+    "else:\n",
+    "    low_bound = None\n",
+    "    high_bound = None\n",
+    "    data_trimmed = data\n",
+    "    quantiles_trim = quantiles_gen\n",
+    "\n",
+    "description = describe(data[\"time\"])\n",
+    "description_trim = describe(data_trimmed[\"time\"])\n",
+    "\n",
+    "max_time = description.minmax[1]\n",
+    "min_time = description.minmax[0]\n",
+    "bit_size = len(bin(max(data[\"priv\"]))) - 2\n",
+    "byte_size = (bit_size + 7) // 8\n",
+    "\n",
+    "if hist_size == \"sqrt\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(sqrt(n)) + 1\n",
+    "elif hist_size == \"sturges\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(log(n, 2)) + 1\n",
+    "elif hist_size == \"rice\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(2 * n**(1/3))\n",
+    "elif hist_size == \"scott\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: (xmax - xmin) // int((3.5 * sqrt(var)) / (n**(1/3)))\n",
+    "elif hist_size == \"fd\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: (xmax - xmin) // int(2 * (xupper - xlower) / (n**(1/3)))\n",
+    "else:\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: hist_size\n",
+    "\n",
+    "hist_size_time = hist_size_func(description.nobs, min_time, max_time, description.variance, quantiles[1], quantiles[3])\n",
+    "hist_size_time_trim = hist_size_func(description_trim.nobs, description_trim.minmax[0], description_trim.minmax[1], description_trim.variance, quantiles_trim[1], quantiles_trim[3])\n",
+    "\n",
+    "if hist_size_time < 30:\n",
+    "    hist_size_time = max_time - min_time\n",
+    "if hist_size_time_trim < 30:\n",
+    "    hist_size_time_trim = description_trim.minmax[1] - description_trim.minmax[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Analysis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Summary"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:57:59.534102Z",
+     "start_time": "2019-03-17T19:57:59.507172Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "display(\"Raw\")\n",
+    "desc = [(\"N\", \"min, max\", \"mean\", \"variance\", \"skewness\", \"kurtosis\"),\n",
+    "        description]\n",
+    "display(HTML(tabulate.tabulate(desc, tablefmt=\"html\")))\n",
+    "display(\"Trimmed\")\n",
+    "desc = [(\"N\", \"min, max\", \"mean\", \"variance\", \"skewness\", \"kurtosis\"),\n",
+    "        description_trim]\n",
+    "display(HTML(tabulate.tabulate(desc, tablefmt=\"html\")))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Selected quantiles"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:58:00.833677Z",
+     "start_time": "2019-03-17T19:58:00.827736Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "tbl = [(quant_low_bound, \"0.25\", \"0.5\", \"0.75\", quant_high_bound),\n",
+    "       list(map(lambda x: \"{} {}\".format(x, time_unit), quantiles))]\n",
+    "display(HTML(tabulate.tabulate(tbl, tablefmt=\"html\")))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Info"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:58:01.954382Z",
+     "start_time": "2019-03-17T19:58:01.947339Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "display(\"Bitsize: {}\".format(bit_size))\n",
+    "display(\"Histogram time bins: {}\".format(hist_size_time))\n",
+    "display(\"Histogram time bins(trimmed): {}\".format(hist_size_time_trim))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plots"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key MSB vs time heatmap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:58:03.641387Z",
+     "start_time": "2019-03-17T19:58:03.572612Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_private = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "axe_private = fig_private.add_subplot(1, 1, 1)\n",
+    "priv_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data_trimmed[\"priv\"])), dtype=np.dtype(\"u1\"))\n",
+    "max_msb = max(priv_msb)\n",
+    "min_msb = min(priv_msb)\n",
+    "heatmap, xedges, yedges = np.histogram2d(priv_msb, data_trimmed[\"time\"],\n",
+    "                                         bins=[max_msb - min_msb + 1, hist_size_time_trim])\n",
+    "extent = [min_msb, max_msb, yedges[0], yedges[-1]]\n",
+    "im = axe_private.imshow(heatmap.T, extent=extent, aspect=\"auto\", cmap=cmap, origin=\"low\",\n",
+    "                   interpolation=\"nearest\", norm=norm)\n",
+    "axe_private.set_xlabel(\"private key MSB value\")\n",
+    "axe_private.set_ylabel(\"key agreement time ({})\".format(time_unit))\n",
+    "fig_private.colorbar(im, ax=axe_private)\n",
+    "\n",
+    "del priv_msb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key Hamming Weight vs time heatmap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:58:07.768683Z",
+     "start_time": "2019-03-17T19:58:06.938237Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_priv_hist = plt.figure(figsize=(10.5, 12), dpi=90)\n",
+    "gs = gridspec.GridSpec(2, 1, height_ratios=[2.5, 1])\n",
+    "axe_priv_hist = fig_priv_hist.add_subplot(gs[0])\n",
+    "axe_priv_hist_hw = fig_priv_hist.add_subplot(gs[1], sharex = axe_priv_hist)\n",
+    "priv_hw = np.array(list(map(hw, data_trimmed[\"priv\"])), dtype=np.dtype(\"u2\"))\n",
+    "h, xe, ye = np.histogram2d(priv_hw, data_trimmed[\"time\"], bins=[max(priv_hw) - min(priv_hw), hist_size_time_trim])\n",
+    "im = axe_priv_hist.imshow(h.T, origin=\"low\", cmap=cmap, aspect=\"auto\", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=norm)\n",
+    "axe_priv_hist.axvline(x=bit_size//2, alpha=0.7, linestyle=\"dotted\", color=\"white\", label=str(bit_size//2) + \" bits\")\n",
+    "axe_priv_hist.set_xlabel(\"private key Hamming weight\")\n",
+    "axe_priv_hist.set_ylabel(\"key agreement time ({})\".format(time_unit))\n",
+    "axe_priv_hist.legend(loc=\"best\")\n",
+    "\n",
+    "plot_hist(axe_priv_hist_hw, priv_hw, \"private key Hamming weight\", log_scale, None)\n",
+    "\n",
+    "param = norm_dist.fit(priv_hw)\n",
+    "pdf_range = np.arange(min(priv_hw), max(priv_hw))\n",
+    "norm_pdf = norm_dist.pdf(pdf_range, *param[:-2], loc=param[-2], scale=param[-1]) * description_trim.nobs\n",
+    "axe_priv_hist_hw.plot(pdf_range, norm_pdf, label=\"fitted normal distribution\")\n",
+    "axe_priv_hist_hw.legend(loc=\"best\")\n",
+    "fig_priv_hist.colorbar(im, ax=[axe_priv_hist, axe_priv_hist_hw])\n",
+    "\n",
+    "display(HTML(\"<b>Private key Hamming weight fitted with normal distribution:</b>\"))\n",
+    "display(HTML(tabulate.tabulate([(\"Mean\", \"Variance\"), param], tablefmt=\"html\")))\n",
+    "\n",
+    "del priv_hw"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Key agreement time histogram"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:58:17.986917Z",
+     "start_time": "2019-03-17T19:58:11.101449Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_ka_hist = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "axe_hist_full = fig_ka_hist.add_subplot(2, 1, 1)\n",
+    "axe_hist_trim = fig_ka_hist.add_subplot(2, 1, 2)\n",
+    "plot_hist(axe_hist_full, data[\"time\"], \"key agreement time ({})\".format(time_unit), log_scale, hist_size_time);\n",
+    "plot_hist(axe_hist_trim, data_trimmed[\"time\"], \"key agreement time ({})\".format(time_unit), log_scale, hist_size_time_trim);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Moving averages of key agreement time"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:51:57.934476Z",
+     "start_time": "2019-03-17T19:51:57.877729Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_avg = plt.figure(figsize=(10.5, 7), dpi=90)\n",
+    "axe_avg = fig_avg.add_subplot(1, 1, 1)\n",
+    "avg_100 = moving_average(data[\"time\"], 100)\n",
+    "avg_1000 = moving_average(data[\"time\"], 1000)\n",
+    "axe_avg.plot(avg_100, label=\"window = 100\")\n",
+    "axe_avg.plot(avg_1000, label=\"window = 1000\")\n",
+    "if low_bound is not None:\n",
+    "    axe_avg.axhline(y=low_bound, alpha=0.7, linestyle=\"dotted\", color=\"green\", label=\"Low trim bound = {}\".format(low_bound))\n",
+    "if high_bound is not None:\n",
+    "    axe_avg.axhline(y=high_bound, alpha=0.7, linestyle=\"dotted\", color=\"orange\", label=\"Hight trim bound = {}\".format(high_bound))\n",
+    "axe_avg.set_ylabel(\"key agreement time ({})\".format(time_unit))\n",
+    "axe_avg.set_xlabel(\"index\")\n",
+    "axe_avg.legend(loc=\"best\")\n",
+    "\n",
+    "del avg_100, avg_1000"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key MSB and LSB histograms"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:51:58.466578Z",
+     "start_time": "2019-03-17T19:51:57.937797Z"
+    },
+    "hide_input": false
+   },
+   "outputs": [],
+   "source": [
+    "fig_priv_hists = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "priv_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data[\"priv\"])), dtype=np.dtype(\"u1\"))\n",
+    "priv_lsb = np.array(list(map(lambda x: x & 0xff, data[\"priv\"])), dtype=np.dtype(\"u1\"))\n",
+    "axe_msb_s_hist = fig_priv_hists.add_subplot(2, 1, 1)\n",
+    "axe_lsb_s_hist = fig_priv_hists.add_subplot(2, 1, 2)\n",
+    "msb_h = plot_hist(axe_msb_s_hist, priv_msb, \"private key MSB\", log_scale, False, False)\n",
+    "lsb_h = plot_hist(axe_lsb_s_hist, priv_lsb, \"private key LSB\", log_scale, False, False)\n",
+    "\n",
+    "del priv_msb, priv_lsb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Public key coordinate MSB and LSB histograms"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:52:21.184705Z",
+     "start_time": "2019-03-17T19:52:20.589707Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def _split(xy):\n",
+    "    x = int.from_bytes(xy[1:byte_size + 1], byteorder=\"big\")\n",
+    "    y = int.from_bytes(xy[1 + byte_size:], byteorder=\"big\")\n",
+    "    return (x, y)\n",
+    "\n",
+    "pub_coords = np.array(list(map(_split, data[\"pub\"])), dtype=np.dtype(\"O\"))\n",
+    "xs = pub_coords[...,0]\n",
+    "ys = pub_coords[...,1]\n",
+    "fig_pub_hists = plt.figure(figsize=(10.5, 14), dpi=90)\n",
+    "\n",
+    "def _plot_coord(data, name, offset):\n",
+    "    axe_msb_pub_hist = fig_pub_hists.add_subplot(4, 1, offset)\n",
+    "    axe_lsb_pub_hist = fig_pub_hists.add_subplot(4, 1, offset + 1)\n",
+    "    pub_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data)))\n",
+    "    pub_lsb = np.array(list(map(lambda x: x & 0xff, data)))\n",
+    "    plot_hist(axe_msb_pub_hist, pub_msb, \"{} coordinate MSB\".format(name), log_scale)\n",
+    "    plot_hist(axe_lsb_pub_hist, pub_lsb, \"{} coordinate LSB\".format(name), log_scale)\n",
+    "    del pub_msb, pub_lsb\n",
+    "\n",
+    "_plot_coord(xs, \"X\", 1)\n",
+    "_plot_coord(ys, \"Y\", 3)\n",
+    "\n",
+    "del pub_coords, xs, ys"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key bit length histogram"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:52:07.657216Z",
+     "start_time": "2019-03-17T19:52:07.549731Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_bl = plt.figure(figsize=(10.5, 12), dpi=90)\n",
+    "gs = gridspec.GridSpec(2, 1, height_ratios=[2.5, 1])\n",
+    "axe_bl_heat = fig_bl.add_subplot(gs[0])\n",
+    "axe_bl_hist = fig_bl.add_subplot(gs[1], sharex=axe_bl_heat)\n",
+    "bl_data = np.array(list(map(lambda x: x.bit_length(), data_trimmed[\"priv\"])), dtype=np.dtype(\"u2\"))\n",
+    "\n",
+    "h, xe, ye = np.histogram2d(bl_data, data_trimmed[\"time\"], bins=[max(bl_data) - min(bl_data), hist_size_time_trim])\n",
+    "im = axe_bl_heat.imshow(h.T, origin=\"low\", cmap=cmap, aspect=\"auto\", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=norm)\n",
+    "axe_bl_heat.set_xlabel(\"private key bit length\")\n",
+    "axe_bl_heat.set_ylabel(\"key agreement time ({})\".format(time_unit))\n",
+    "\n",
+    "plot_hist(axe_bl_hist, bl_data, \"Private key bit length\", log_scale, align=\"right\")\n",
+    "fig_bl.colorbar(im, ax=[axe_bl_heat, axe_bl_hist])\n",
+    "\n",
+    "del bl_data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "@webio": {
+   "lastCommId": "954c1f99782e402895d668a42553e22f",
+   "lastKernelId": "0b8e59f0-d640-4f72-ae7f-1b327e75910b"
+  },
+  "hide_input": false,
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.7.2"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/util/plot_dh.py b/util/plot_dh.py
deleted file mode 100755
index 60e20ae..0000000
--- a/util/plot_dh.py
+++ /dev/null
@@ -1,107 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: UTF-8 -*-
-#
-# Script for plotting ECTester ECDH results.
-#
-# Example usage:
-#
-#     > java -jar ECTesterReader.jar -dh 10000 -b 192 -fp -o dh.csv
-#     ...
-#     > ./plot_dh.py dh.csv
-#     ...
-#
-
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib import ticker, colors
-import argparse
-from copy import deepcopy
-from operator import itemgetter
-
-from utils import hw, moving_average, plot_hist
-
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser(description="Plot ECTester ECDH timing.")
-    parser.add_argument("-o", "--output", dest="output", type=argparse.FileType("wb"), help="Write image to [file], do not display.", metavar="file")
-    parser.add_argument("--priv", dest="priv", action="store_true", help="Show private key MSB heatmap plot.")
-    parser.add_argument("--hist", dest="hist", action="store_true", help="Show time histogram.")
-    parser.add_argument("--hw-hist", dest="hw_hist", action="store_true", help="Show Hamming weight heatmap (private key Hamming weight and time).")
-    parser.add_argument("--avg", dest="avg", action="store_true", help="Show moving average of time.")
-    parser.add_argument("--log", dest="log", action="store_true", help="Use logarithmic scale.")
-    parser.add_argument("--skip-first", dest="skip_first", nargs="?", const=1, type=int, help="Skip first entry, as it's usually a large outlier.")
-    parser.add_argument("-t", "--title", dest="title", nargs="?", default="", type=str, help="What title to give the figure.")
-    parser.add_argument("file", type=str, help="The file to plot(csv).")
-
-    opts = parser.parse_args()
-
-    with open(opts.file, "r") as f:
-        header = f.readline()
-    header_names = header.split(";")
-
-    hx = lambda x: int(x, 16)
-    data = np.genfromtxt(opts.file, delimiter=";", skip_header=1, converters={2: hx, 3: hx, 4: hx}, dtype=np.dtype([("index","u4"), ("time","u4"), ("pub", "O"), ("priv", "O"), ("secret","O")]))
-    if opts.skip_first:
-        data = data[opts.skip_first:]
-
-    time_data = data["time"]
-    if "nano" in header_names[1]:
-        unit = r"$\mu s$"
-        time_data = np.array(list(map(lambda x: x//1000, time_data)))
-    else:
-        unit = r"ms"
-    priv_data = data["priv"]
-    pub_data = data["pub"]
-    secret_data = data["secret"]
-
-    plt.style.use("ggplot")
-    fig = plt.figure()
-    layout_kwargs = {}
-    if opts.title is None:
-        fig.suptitle(opts.file)
-        layout_kwargs["rect"] = [0, 0.02, 1, 0.98]
-    elif opts.title:
-        fig.suptitle(opts.title)
-        layout_kwargs["rect"] = [0, 0.02, 1, 0.98]
-    fig.tight_layout(**layout_kwargs)
-
-    time_max = max(time_data)
-    time_min = min(time_data)
-    bit_size = len(bin(max(priv_data))) - 2
-
-    cmap = deepcopy(plt.cm.plasma)
-    cmap.set_bad("black")
-
-    norm = colors.Normalize()
-    if opts.log:
-        norm = colors.LogNorm()
-
-    axe_private = fig.add_subplot(3,1,1)
-    priv_msb = np.array(list(map(lambda x: x >> (bit_size - 8), priv_data)), dtype=np.dtype("u1"))
-    heatmap, xedges, yedges = np.histogram2d(priv_msb, time_data, bins=[128, time_max - time_min])
-    extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
-    axe_private.imshow(heatmap.T, extent=extent, aspect="auto", cmap=cmap, origin="low", interpolation="nearest", norm=norm)
-    axe_private.set_xlabel("private key MSB value")
-    axe_private.set_ylabel("ECDH time ({})".format(unit))
-
-    axe_hist = fig.add_subplot(3,1,2)
-    plot_hist(axe_hist, time_data, "ECDH time ({})".format(unit), opts.log)
-    axe_hist.legend(loc="best")
-
-    axe_priv_hist = fig.add_subplot(3,1,3)
-    priv_hw = np.array(list(map(hw, priv_data)), dtype=np.dtype("u2"))
-    h, xe, ye = np.histogram2d(priv_hw, time_data, bins=[max(priv_hw) - min(priv_hw), time_max - time_min])
-    im = axe_priv_hist.imshow(h.T, origin="low", cmap=cmap, aspect="auto", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=colors.LogNorm())
-    axe_priv_hist.axvline(x=bit_size//2, alpha=0.7, linestyle="dotted", color="white", label=str(bit_size//2) + " bits")    
-    axe_priv_hist.set_xlabel("private key Hamming weight")
-    axe_priv_hist.set_ylabel("time ({})".format(unit))
-    axe_priv_hist.legend(loc="best")
-    fig.colorbar(im, ax=axe_priv_hist)
-
-    fig.text(0.01, 0.02, "Data size: {}".format(len(time_data)), size="small")
-
-    if opts.output is None:
-        plt.show()
-    else:
-        fig.set_size_inches(12, 10)
-        ext = opts.output.name.split(".")[-1]
-        plt.savefig(opts.output, format=ext, dpi=400, bbox_inches='tight')
diff --git a/util/plot_dsa.ipynb b/util/plot_dsa.ipynb
new file mode 100644
index 0000000..dbfb38b
--- /dev/null
+++ b/util/plot_dsa.ipynb
@@ -0,0 +1,613 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Analysis of key generation data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:00:25.518989Z",
+     "start_time": "2019-03-17T23:00:24.501601Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    " %matplotlib notebook \n",
+    "import numpy as np\n",
+    "from scipy.stats import describe\n",
+    "from scipy.stats import norm as norm_dist\n",
+    "from scipy.stats.mstats import mquantiles\n",
+    "from math import log, sqrt\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import ticker, colors, gridspec\n",
+    "from copy import deepcopy\n",
+    "from utils import plot_hist, moving_average, hw\n",
+    "from binascii import unhexlify\n",
+    "from IPython.display import display, HTML\n",
+    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
+    "import ipywidgets as widgets\n",
+    "import tabulate"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Settings\n",
+    "Enter your input below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:06:29.704432Z",
+     "start_time": "2019-03-17T23:06:29.694540Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# File name with output from ECTesterReader or ECTesterStandalone signatures.\n",
+    "fname = \"filename.csv\"\n",
+    "\n",
+    "# The amount of entries skipped from the beginning of the file, as they are usually outliers.\n",
+    "skip_first = 10\n",
+    "\n",
+    "# Whether to plot things in logarithmic scale or not.\n",
+    "log_scale = False\n",
+    "\n",
+    "# Whether to trim the time data outside the 1 - 99 percentile range (adjust below). Quite useful.\n",
+    "trim = True\n",
+    "\n",
+    "# How much to trim? Either a number in [0,1] signifying a quantile, or an absolute value signifying a threshold\n",
+    "trim_low = 0.01\n",
+    "trim_high = 0.99\n",
+    "\n",
+    "# Graphical (matplotlib) style name\n",
+    "style = \"ggplot\"\n",
+    "\n",
+    "# Color map to use, and what color to assign to \"bad\" values (necessary for log_scale)\n",
+    "color_map = plt.cm.plasma\n",
+    "color_map_bad = \"black\"\n",
+    "\n",
+    "# What function to use to calculate number of histogram bins of time\n",
+    "# one of \"sqrt\", \"sturges\", \"rice\", \"scott\" and \"fd\" or a number specifying the number of bins\n",
+    "hist_size = \"sturges\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data processing"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:06:30.551732Z",
+     "start_time": "2019-03-17T23:06:30.545202Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Setup plot style\n",
+    "\n",
+    "plt.style.use(style)\n",
+    "\n",
+    "cmap = deepcopy(color_map)\n",
+    "cmap.set_bad(color_map_bad)\n",
+    "\n",
+    "# Normalization, linear or log.\n",
+    "if log_scale:\n",
+    "    norm = colors.LogNorm()\n",
+    "else:\n",
+    "    norm = colors.Normalize()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:00:38.023486Z",
+     "start_time": "2019-03-17T23:00:27.178465Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Read the header line.\n",
+    "\n",
+    "with open(fname, \"r\") as f:\n",
+    "    header = f.readline()\n",
+    "header_names = header.split(\";\")\n",
+    "if len(header_names) != 9:\n",
+    "    print(\"Bad data?\")\n",
+    "    exit(1)\n",
+    "\n",
+    "# Load the data\n",
+    "\n",
+    "hx = lambda x: int(x, 16)\n",
+    "data = np.genfromtxt(fname, delimiter=\";\", skip_header=1, converters={3: unhexlify, 4: unhexlify,\n",
+    "                                                                      5: hx, 6: unhexlify, 7: hx,\n",
+    "                                                                      8: lambda b: bool(int(b))},\n",
+    "                     dtype=np.dtype([(\"index\", \"u4\"), (\"sign_time\", \"u4\"), (\"verify_time\", \"u4\"),\n",
+    "                                     (\"data\", \"O\"), (\"pub\", \"O\"), (\"priv\", \"O\"), (\"signature\", \"O\"),\n",
+    "                                    (\"nonce\", \"O\"), (\"valid\", \"b\")]))\n",
+    "\n",
+    "    \n",
+    "sign_unit = \"ms\"\n",
+    "verify_unit = \"ms\"\n",
+    "# Setup the datatrimmed = False"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:00:38.465677Z",
+     "start_time": "2019-03-17T23:00:38.025692Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Setup the data\n",
+    "\n",
+    "# Skip first (outliers?)\n",
+    "\n",
+    "data = data[skip_first:]\n",
+    "\n",
+    "# If in nanoseconds, scale to microseconds\n",
+    "if header_names[1].endswith(\"[nano]\") and sign_unit == \"ms\":\n",
+    "    sign_unit = r\"$\\mu s$\"\n",
+    "    np.floor_divide(data[\"sign_time\"], 1000, out=data[\"sign_time\"])\n",
+    "\n",
+    "if header_names[2].endswith(\"[nano]\") and verify_unit == \"ms\":\n",
+    "    verify_unit = r\"$\\mu s$\"\n",
+    "    np.floor_divide(data[\"verify_time\"], 1000, out=data[\"verify_time\"])\n",
+    "\n",
+    "# Trim times\n",
+    "quant_low_bound = trim_low if 0 <= trim_low <= 1 else 0.01\n",
+    "quant_high_bound = trim_high if 0 <= trim_high <= 1 else 0.95\n",
+    "quantiles_sign = mquantiles(data[\"sign_time\"], prob=(quant_low_bound, 0.25, 0.5, 0.75, quant_high_bound))\n",
+    "if trim:\n",
+    "    low_bound = quantiles_sign[0] if 0 <= trim_low <= 1 else trim_low\n",
+    "    high_bound = quantiles_sign[4] if 0 <= trim_high <= 1 else trim_high\n",
+    "    data_trimmed = data[np.logical_and(data[\"sign_time\"] >= low_bound,\n",
+    "                                       data[\"sign_time\"] <= high_bound)]\n",
+    "    quantiles_sign_trim = mquantiles(data_trimmed[\"sign_time\"], prob=(quant_low_bound, 0.25, 0.5, 0.75, quant_high_bound))\n",
+    "else:\n",
+    "    low_bound = None\n",
+    "    high_bound = None\n",
+    "    data_trimmed = data\n",
+    "    quantiles_sign_trim = quantiles_sign\n",
+    "\n",
+    "description_sign = describe(data[\"sign_time\"])\n",
+    "description_sign_trim = describe(data_trimmed[\"sign_time\"])\n",
+    "\n",
+    "max_sign_time = description_sign.minmax[1]\n",
+    "min_sign_time = description_sign.minmax[0]\n",
+    "bit_size = len(bin(max(data[\"priv\"]))) - 2\n",
+    "byte_size = (bit_size + 7) // 8\n",
+    "bit_size = byte_size * 8\n",
+    "\n",
+    "if hist_size == \"sqrt\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(sqrt(n)) + 1\n",
+    "elif hist_size == \"sturges\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(log(n, 2)) + 1\n",
+    "elif hist_size == \"rice\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(2 * n**(1/3))\n",
+    "elif hist_size == \"scott\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: (xmax - xmin) // int((3.5 * sqrt(var)) / (n**(1/3)))\n",
+    "elif hist_size == \"fd\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: (xmax - xmin) // int(2 * (xupper - xlower) / (n**(1/3)))\n",
+    "else:\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: hist_size\n",
+    "\n",
+    "hist_size_sign_time = hist_size_func(description_sign.nobs, min_sign_time, max_sign_time, description_sign.variance, quantiles_sign[1], quantiles_sign[3])\n",
+    "hist_size_sign_time_trim = hist_size_func(description_sign_trim.nobs, description_sign_trim.minmax[0], description_sign_trim.minmax[1], description_sign_trim.variance, quantiles_sign_trim[1], quantiles_sign_trim[3])\n",
+    "\n",
+    "if hist_size_sign_time < 30:\n",
+    "    hist_size_sign_time = max_sign_time - min_sign_time\n",
+    "if hist_size_sign_time_trim < 30:\n",
+    "    hist_size_sign_time_trim = description_sign_trim.minmax[1] - description_sign_trim.minmax[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Analysis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Summary"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:00:39.540701Z",
+     "start_time": "2019-03-17T23:00:39.511019Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "display(\"Raw\")\n",
+    "desc = [(\"N\", \"min, max\", \"mean\", \"variance\", \"skewness\", \"kurtosis\"),\n",
+    "        description_sign]\n",
+    "display(HTML(tabulate.tabulate(desc, tablefmt=\"html\")))\n",
+    "display(\"Trimmed\")\n",
+    "desc = [(\"N\", \"min, max\", \"mean\", \"variance\", \"skewness\", \"kurtosis\"),\n",
+    "        description_sign_trim]\n",
+    "display(HTML(tabulate.tabulate(desc, tablefmt=\"html\")))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Selected quantiles"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:00:40.974497Z",
+     "start_time": "2019-03-17T23:00:40.953755Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "tbl = [(quant_low_bound, \"0.25\", \"0.5\", \"0.75\", quant_high_bound),\n",
+    "       list(map(lambda x: \"{} {}\".format(x, sign_unit), quantiles_sign))]\n",
+    "display(HTML(tabulate.tabulate(tbl, tablefmt=\"html\")))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Info"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:00:41.961541Z",
+     "start_time": "2019-03-17T23:00:41.949385Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "display(\"Bitsize:\", bit_size)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plots"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Nonce MSB vs signature time heatmap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:06:34.030472Z",
+     "start_time": "2019-03-17T23:06:33.761991Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_nonce = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "axe_nonce = fig_nonce.add_subplot(1, 1, 1)\n",
+    "nonce_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data_trimmed[\"nonce\"])), dtype=np.dtype(\"u1\"))\n",
+    "max_msb = max(nonce_msb)\n",
+    "min_msb = min(nonce_msb)\n",
+    "heatmap, xedges, yedges = np.histogram2d(nonce_msb, data_trimmed[\"sign_time\"],\n",
+    "                                         bins=[max_msb - min_msb + 1, hist_size_sign_time_trim])\n",
+    "extent = [min_msb, max_msb, yedges[0], yedges[-1]]\n",
+    "im = axe_nonce.imshow(heatmap.T, extent=extent, aspect=\"auto\", cmap=cmap, origin=\"low\",\n",
+    "                   interpolation=\"nearest\", norm=norm)\n",
+    "axe_nonce.set_xlabel(\"nonce key MSB value\")\n",
+    "axe_nonce.set_ylabel(\"signature time ({})\".format(sign_unit))\n",
+    "fig_nonce.colorbar(im, ax=axe_nonce)\n",
+    "\n",
+    "del nonce_msb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Nonce Hamming Weight vs signature time heatmap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:05:26.804859Z",
+     "start_time": "2019-03-17T23:05:18.214110Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_nonce_hist = plt.figure(figsize=(10.5, 12), dpi=90)\n",
+    "gs = gridspec.GridSpec(2, 1, height_ratios=[2.5, 1])\n",
+    "axe_nonce_hist = fig_nonce_hist.add_subplot(gs[0])\n",
+    "axe_nonce_hist_hw = fig_nonce_hist.add_subplot(gs[1], sharex = axe_nonce_hist)\n",
+    "nonce_hw = np.array(list(map(hw, data_trimmed[\"nonce\"])), dtype=np.dtype(\"u2\"))\n",
+    "h, xe, ye = np.histogram2d(nonce_hw, data_trimmed[\"sign_time\"], bins=[max(nonce_hw) - min(nonce_hw), hist_size_sign_time_trim])\n",
+    "im = axe_nonce_hist.imshow(h.T, origin=\"low\", cmap=cmap, aspect=\"auto\", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=norm)\n",
+    "axe_nonce_hist.axvline(x=bit_size//2, alpha=0.7, linestyle=\"dotted\", color=\"white\", label=str(bit_size//2) + \" bits\")\n",
+    "axe_nonce_hist.set_xlabel(\"nonce Hamming weight\")\n",
+    "axe_nonce_hist.set_ylabel(\"signature time ({})\".format(sign_unit))\n",
+    "axe_nonce_hist.legend(loc=\"best\")\n",
+    "\n",
+    "plot_hist(axe_nonce_hist_hw, nonce_hw, \"nonce Hamming weight\", log_scale, True, True)\n",
+    "\n",
+    "param = norm_dist.fit(nonce_hw)\n",
+    "pdf_range = np.arange(min(nonce_hw), max(nonce_hw))\n",
+    "norm_pdf = norm_dist.pdf(pdf_range, *param[:-2], loc=param[-2], scale=param[-1]) * description_sign_trim.nobs\n",
+    "axe_nonce_hist_hw.plot(pdf_range, norm_pdf, label=\"fitted normal distribution\")\n",
+    "axe_nonce_hist_hw.legend(loc=\"best\")\n",
+    "fig_nonce_hist.colorbar(im, ax=[axe_nonce_hist, axe_nonce_hist_hw])\n",
+    "\n",
+    "display(HTML(\"<b>Nonce Hamming weight fitted with normal distribution:</b>\"))\n",
+    "display(HTML(tabulate.tabulate([(\"Mean\", \"Variance\"), param], tablefmt=\"html\")))\n",
+    "\n",
+    "del nonce_hw"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Signature time histogram"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:05:32.395983Z",
+     "start_time": "2019-03-17T23:05:32.068823Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_sig_hist = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "axe_hist_full = fig_sig_hist.add_subplot(2, 1, 1)\n",
+    "axe_hist_trim = fig_sig_hist.add_subplot(2, 1, 2)\n",
+    "plot_hist(axe_hist_full, data[\"sign_time\"], \"signature time ({})\".format(sign_unit), log_scale, hist_size_sign_time);\n",
+    "plot_hist(axe_hist_trim, data_trimmed[\"sign_time\"], \"signature time ({})\".format(sign_unit), log_scale, hist_size_sign_time_trim);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Verification time histogram"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:05:33.358613Z",
+     "start_time": "2019-03-17T23:05:32.963791Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_vrfy_hist = plt.figure(figsize=(10.5, 5), dpi=90)\n",
+    "axe_hist_full = fig_vrfy_hist.add_subplot(1, 1, 1)\n",
+    "plot_hist(axe_hist_full, data[\"verify_time\"], \"verification time ({})\".format(verify_unit), log_scale, hist_size_sign_time);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Moving averages of signature and verification times"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:05:33.971385Z",
+     "start_time": "2019-03-17T23:05:33.732857Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_avg = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "axe_sign_avg = fig_avg.add_subplot(2, 1, 1)\n",
+    "axe_vrfy_avg = fig_avg.add_subplot(2, 1, 2, sharex=axe_sign_avg)\n",
+    "avg_sign_100 = moving_average(data[\"sign_time\"], 100)\n",
+    "avg_sign_1000 = moving_average(data[\"sign_time\"], 1000)\n",
+    "axe_sign_avg.plot(avg_sign_100, label=\"window = 100\")\n",
+    "axe_sign_avg.plot(avg_sign_1000, label=\"window = 1000\")\n",
+    "if low_bound is not None:\n",
+    "    axe_sign_avg.axhline(y=low_bound, alpha=0.7, linestyle=\"dotted\", color=\"green\", label=\"Low trim bound = {}\".format(low_bound))\n",
+    "if high_bound is not None:\n",
+    "    axe_sign_avg.axhline(y=high_bound, alpha=0.7, linestyle=\"dotted\", color=\"orange\", label=\"Hight trim bound = {}\".format(high_bound))\n",
+    "axe_sign_avg.set_ylabel(\"signature time ({})\".format(sign_unit))\n",
+    "axe_sign_avg.set_xlabel(\"index\")\n",
+    "axe_sign_avg.legend(loc=\"best\")\n",
+    "\n",
+    "avg_vrfy_100 = moving_average(data[\"verify_time\"], 100)\n",
+    "avg_vrfy_1000 = moving_average(data[\"verify_time\"], 1000)\n",
+    "axe_vrfy_avg.plot(avg_vrfy_100, label=\"window = 100\")\n",
+    "axe_vrfy_avg.plot(avg_vrfy_1000, label=\"window = 1000\")\n",
+    "axe_vrfy_avg.set_ylabel(\"verification time ({})\".format(verify_unit))\n",
+    "axe_vrfy_avg.set_xlabel(\"index\")\n",
+    "axe_vrfy_avg.legend(loc=\"best\")\n",
+    "\n",
+    "del avg_sign_100, avg_sign_1000, avg_vrfy_100, avg_vrfy_1000"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Nonce MSB and LSB histograms"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:05:36.256032Z",
+     "start_time": "2019-03-17T23:05:35.302194Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_nonce_hists = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "nonce_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data[\"nonce\"])), dtype=np.dtype(\"u1\"))\n",
+    "nonce_lsb = np.array(list(map(lambda x: x & 0xff, data[\"nonce\"])), dtype=np.dtype(\"u1\"))\n",
+    "axe_msb_n_hist = fig_nonce_hists.add_subplot(2, 1, 1)\n",
+    "axe_lsb_n_hist = fig_nonce_hists.add_subplot(2, 1, 2)\n",
+    "plot_hist(axe_msb_n_hist, nonce_msb, \"nonce MSB\", log_scale, False, False)\n",
+    "plot_hist(axe_lsb_n_hist, nonce_lsb, \"nonce LSB\", log_scale, False, False)\n",
+    "\n",
+    "del nonce_msb, nonce_lsb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Nonce bit length histogram"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T23:05:45.320760Z",
+     "start_time": "2019-03-17T23:05:44.951189Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_bl = plt.figure(figsize=(10.5, 12), dpi=90)\n",
+    "gs = gridspec.GridSpec(2, 1, height_ratios=[2.5, 1])\n",
+    "axe_bl_heat = fig_bl.add_subplot(gs[0])\n",
+    "axe_bl_hist = fig_bl.add_subplot(gs[1], sharex=axe_bl_heat)\n",
+    "bl_data = np.array(list(map(lambda x: x.bit_length(), data_trimmed[\"nonce\"])), dtype=np.dtype(\"u2\"))\n",
+    "\n",
+    "h, xe, ye = np.histogram2d(bl_data, data_trimmed[\"sign_time\"], bins=[max(bl_data) - min(bl_data), hist_size_sign_time_trim])\n",
+    "im = axe_bl_heat.imshow(h.T, origin=\"low\", cmap=cmap, aspect=\"auto\", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=norm)\n",
+    "axe_bl_heat.set_xlabel(\"nonce bit length\")\n",
+    "axe_bl_heat.set_ylabel(\"signature time ({})\".format(sign_unit))\n",
+    "\n",
+    "plot_hist(axe_bl_hist, bl_data, \"nonce bit length\", log_scale, align=\"right\")\n",
+    "fig_bl.colorbar(im, ax=[axe_bl_heat, axe_bl_hist])\n",
+    "\n",
+    "del bl_data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "@webio": {
+   "lastCommId": "2218f80b2f784436bce6ffed5d971ea3",
+   "lastKernelId": "1c06331a-17dd-4743-9e31-832000e597c8"
+  },
+  "hide_input": false,
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.7.2"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/util/plot_gen.ipynb b/util/plot_gen.ipynb
new file mode 100644
index 0000000..5f02a00
--- /dev/null
+++ b/util/plot_gen.ipynb
@@ -0,0 +1,649 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Analysis of key generation data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:38.893311Z",
+     "start_time": "2019-03-17T19:16:37.845017Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    " %matplotlib notebook \n",
+    "import numpy as np\n",
+    "from scipy.stats import describe\n",
+    "from scipy.stats import norm as norm_dist\n",
+    "from scipy.stats.mstats import mquantiles\n",
+    "from math import log, sqrt\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import ticker, colors, gridspec\n",
+    "from copy import deepcopy\n",
+    "from utils import plot_hist, moving_average, hw\n",
+    "from binascii import unhexlify\n",
+    "from IPython.display import display, HTML\n",
+    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
+    "import ipywidgets as widgets\n",
+    "import tabulate"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Settings\n",
+    "Enter your input below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:38.911566Z",
+     "start_time": "2019-03-17T19:16:38.900168Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# File name with output from ECTesterReader or ECTesterStandalone key generation.\n",
+    "fname = \"filename.csv\"\n",
+    "\n",
+    "# The amount of entries skipped from the beginning of the file, as they are usually outliers.\n",
+    "skip_first = 10\n",
+    "\n",
+    "# Whether to plot things in logarithmic scale or not.\n",
+    "log_scale = False\n",
+    "\n",
+    "# Whether to trim the time data outside the 1 - 99 percentile range (adjust below). Quite useful.\n",
+    "trim = True\n",
+    "\n",
+    "# How much to trim? Either a number in [0,1] signifying a quantile, or an absolute value signifying a threshold\n",
+    "trim_low = 0.01\n",
+    "trim_high = 0.99\n",
+    "\n",
+    "# Graphical (matplotlib) style name\n",
+    "style = \"ggplot\"\n",
+    "\n",
+    "# Color map to use, and what color to assign to \"bad\" values (necessary for log_scale)\n",
+    "color_map = plt.cm.plasma\n",
+    "color_map_bad = \"black\"\n",
+    "\n",
+    "# What function to use to calculate number of histogram bins of time\n",
+    "# one of \"sqrt\", \"sturges\", \"rice\", \"scott\" and \"fd\" or a number specifying the number of bins\n",
+    "hist_size = \"sturges\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data processing"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:39.733575Z",
+     "start_time": "2019-03-17T19:16:39.728385Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Setup plot style\n",
+    "\n",
+    "plt.style.use(style)\n",
+    "\n",
+    "cmap = deepcopy(color_map)\n",
+    "cmap.set_bad(color_map_bad)\n",
+    "\n",
+    "# Normalization, linear or log.\n",
+    "if log_scale:\n",
+    "    norm = colors.LogNorm()\n",
+    "else:\n",
+    "    norm = colors.Normalize()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:42.300146Z",
+     "start_time": "2019-03-17T19:16:40.259135Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Read the header line.\n",
+    "\n",
+    "with open(fname, \"r\") as f:\n",
+    "    header = f.readline()\n",
+    "header_names = header.split(\";\")\n",
+    "if len(header_names) not in (4, 5):\n",
+    "    print(\"Bad data?\")\n",
+    "    exit(1)\n",
+    "\n",
+    "# Load the data\n",
+    "\n",
+    "hx = lambda x: int(x, 16)\n",
+    "if len(header_names) == 4:\n",
+    "    data = np.genfromtxt(fname, delimiter=\";\", skip_header=1, converters={2: unhexlify, 3: hx},\n",
+    "                         dtype=np.dtype([(\"index\", \"u4\"), (\"gen_time\", \"u4\"), (\"pub\", \"O\"), (\"priv\", \"O\")]))\n",
+    "else:\n",
+    "    data = np.genfromtxt(fname, delimiter=\";\", skip_header=1, converters={3: unhexlify, 4: hx},\n",
+    "                         dtype=np.dtype([(\"index\", \"u4\"), (\"gen_time\", \"u4\"), (\"export_time\", \"u4\"),\n",
+    "                                         (\"pub\", \"O\"), (\"priv\", \"O\")]))\n",
+    "    \n",
+    "gen_unit = \"ms\"\n",
+    "export_unit = \"ms\"\n",
+    "# Setup the datatrimmed = False"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:42.417415Z",
+     "start_time": "2019-03-17T19:16:42.302353Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Setup the data\n",
+    "\n",
+    "# Skip first (outliers?)\n",
+    "\n",
+    "data = data[skip_first:]\n",
+    "\n",
+    "# If in nanoseconds, scale to microseconds\n",
+    "if header_names[1].endswith(\"[nano]\") and gen_unit == \"ms\":\n",
+    "    gen_unit = r\"$\\mu s$\"\n",
+    "    np.floor_divide(data[\"gen_time\"], 1000, out=data[\"gen_time\"])\n",
+    "\n",
+    "if len(header_names) == 5 and header_names[2].endswith(\"[nano]\") and export_unit == \"ms\":\n",
+    "    export_unit = r\"$\\mu s$\"\n",
+    "    np.floor_divide(data[\"export_time\"], 1000, out=data[\"export_time\"])\n",
+    "\n",
+    "# Trim gen times\n",
+    "quant_low_bound = trim_low if 0 <= trim_low <= 1 else 0.01\n",
+    "quant_high_bound = trim_high if 0 <= trim_high <= 1 else 0.95\n",
+    "quantiles_gen = mquantiles(data[\"gen_time\"], prob=(quant_low_bound, 0.25, 0.5, 0.75, quant_high_bound))\n",
+    "if trim:\n",
+    "    low_bound = quantiles_gen[0] if 0 <= trim_low <= 1 else trim_low\n",
+    "    high_bound = quantiles_gen[4] if 0 <= trim_high <= 1 else trim_high\n",
+    "    data_trimmed = data[np.logical_and(data[\"gen_time\"] >= low_bound,\n",
+    "                                       data[\"gen_time\"] <= high_bound)]\n",
+    "    quantiles_gen_trim = mquantiles(data_trimmed[\"gen_time\"], prob=(quant_low_bound, 0.25, 0.5, 0.75, quant_high_bound))\n",
+    "else:\n",
+    "    low_bound = None\n",
+    "    high_bound = None\n",
+    "    data_trimmed = data\n",
+    "    quantiles_gen_trim = quantiles_gen\n",
+    "\n",
+    "description_gen = describe(data[\"gen_time\"])\n",
+    "description_gen_trim = describe(data_trimmed[\"gen_time\"])\n",
+    "\n",
+    "max_gen_time = description_gen.minmax[1]\n",
+    "min_gen_time = description_gen.minmax[0]\n",
+    "bit_size = len(bin(max(data[\"priv\"]))) - 2\n",
+    "byte_size = (bit_size + 7) // 8\n",
+    "\n",
+    "if hist_size == \"sqrt\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(sqrt(n)) + 1\n",
+    "elif hist_size == \"sturges\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(log(n, 2)) + 1\n",
+    "elif hist_size == \"rice\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: int(2 * n**(1/3))\n",
+    "elif hist_size == \"scott\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: (xmax - xmin) // int((3.5 * sqrt(var)) / (n**(1/3)))\n",
+    "elif hist_size == \"fd\":\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: (xmax - xmin) // int(2 * (xupper - xlower) / (n**(1/3)))\n",
+    "else:\n",
+    "    hist_size_func = lambda n, xmin, xmax, var, xlower, xupper: hist_size\n",
+    "\n",
+    "hist_size_gen_time = hist_size_func(description_gen.nobs, min_gen_time, max_gen_time, description_gen.variance, quantiles_gen[1], quantiles_gen[3])\n",
+    "hist_size_gen_time_trim = hist_size_func(description_gen_trim.nobs, description_gen_trim.minmax[0], description_gen_trim.minmax[1], description_gen_trim.variance, quantiles_gen_trim[1], quantiles_gen_trim[3])\n",
+    "\n",
+    "if hist_size_gen_time < 30:\n",
+    "    hist_size_gen_time = max_gen_time - min_gen_time\n",
+    "if hist_size_gen_time_trim < 30:\n",
+    "    hist_size_gen_time_trim = description_gen_trim.minmax[1] - description_gen_trim.minmax[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Analysis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Summary"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:43.343937Z",
+     "start_time": "2019-03-17T19:16:43.329900Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "display(\"Raw\")\n",
+    "desc = [(\"N\", \"min, max\", \"mean\", \"variance\", \"skewness\", \"kurtosis\"),\n",
+    "        description_gen]\n",
+    "display(HTML(tabulate.tabulate(desc, tablefmt=\"html\")))\n",
+    "display(\"Trimmed\")\n",
+    "desc = [(\"N\", \"min, max\", \"mean\", \"variance\", \"skewness\", \"kurtosis\"),\n",
+    "        description_gen_trim]\n",
+    "display(HTML(tabulate.tabulate(desc, tablefmt=\"html\")))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Selected quantiles"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:44.058425Z",
+     "start_time": "2019-03-17T19:16:44.043877Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "tbl = [(quant_low_bound, \"0.25\", \"0.5\", \"0.75\", quant_high_bound),\n",
+    "       list(map(lambda x: \"{} {}\".format(x, gen_unit), quantiles_gen))]\n",
+    "display(HTML(tabulate.tabulate(tbl, tablefmt=\"html\")))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Info"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:44.688872Z",
+     "start_time": "2019-03-17T19:16:44.684485Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "display(\"Bitsize:\", bit_size)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plots"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key MSB vs time heatmap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:45.995145Z",
+     "start_time": "2019-03-17T19:16:45.802741Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_private = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "axe_private = fig_private.add_subplot(1, 1, 1)\n",
+    "priv_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data_trimmed[\"priv\"])), dtype=np.dtype(\"u1\"))\n",
+    "max_msb = max(priv_msb)\n",
+    "min_msb = min(priv_msb)\n",
+    "heatmap, xedges, yedges = np.histogram2d(priv_msb, data_trimmed[\"gen_time\"],\n",
+    "                                         bins=[max_msb - min_msb + 1, hist_size_gen_time_trim])\n",
+    "extent = [min_msb, max_msb, yedges[0], yedges[-1]]\n",
+    "im = axe_private.imshow(heatmap.T, extent=extent, aspect=\"auto\", cmap=cmap, origin=\"low\",\n",
+    "                   interpolation=\"nearest\", norm=norm)\n",
+    "axe_private.set_xlabel(\"private key MSB value\")\n",
+    "axe_private.set_ylabel(\"keygen time ({})\".format(gen_unit))\n",
+    "fig_private.colorbar(im, ax=axe_private)\n",
+    "\n",
+    "del priv_msb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key Hamming Weight vs time heatmap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:49.890330Z",
+     "start_time": "2019-03-17T19:16:47.357225Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_priv_hist = plt.figure(figsize=(10.5, 12), dpi=90)\n",
+    "gs = gridspec.GridSpec(2, 1, height_ratios=[2.5, 1])\n",
+    "axe_priv_hist = fig_priv_hist.add_subplot(gs[0])\n",
+    "axe_priv_hist_hw = fig_priv_hist.add_subplot(gs[1], sharex = axe_priv_hist)\n",
+    "priv_hw = np.array(list(map(hw, data_trimmed[\"priv\"])), dtype=np.dtype(\"u2\"))\n",
+    "h, xe, ye = np.histogram2d(priv_hw, data_trimmed[\"gen_time\"], bins=[max(priv_hw) - min(priv_hw), hist_size_gen_time_trim])\n",
+    "im = axe_priv_hist.imshow(h.T, origin=\"low\", cmap=cmap, aspect=\"auto\", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=norm)\n",
+    "axe_priv_hist.axvline(x=bit_size//2, alpha=0.7, linestyle=\"dotted\", color=\"white\", label=str(bit_size//2) + \" bits\")\n",
+    "axe_priv_hist.set_xlabel(\"private key Hamming weight\")\n",
+    "axe_priv_hist.set_ylabel(\"keygen time ({})\".format(gen_unit))\n",
+    "axe_priv_hist.legend(loc=\"best\")\n",
+    "\n",
+    "plot_hist(axe_priv_hist_hw, priv_hw, \"private key Hamming weight\", log_scale, None)\n",
+    "\n",
+    "param = norm_dist.fit(priv_hw)\n",
+    "pdf_range = np.arange(min(priv_hw), max(priv_hw))\n",
+    "norm_pdf = norm_dist.pdf(pdf_range, *param[:-2], loc=param[-2], scale=param[-1]) * description_gen_trim.nobs\n",
+    "axe_priv_hist_hw.plot(pdf_range, norm_pdf, label=\"fitted normal distribution\")\n",
+    "axe_priv_hist_hw.legend(loc=\"best\")\n",
+    "fig_priv_hist.colorbar(im, ax=[axe_priv_hist, axe_priv_hist_hw])\n",
+    "\n",
+    "display(HTML(\"<b>Private key Hamming weight fitted with normal distribution:</b>\"))\n",
+    "display(HTML(tabulate.tabulate([(\"Mean\", \"Variance\"), param], tablefmt=\"html\")))\n",
+    "\n",
+    "del priv_hw"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Key generation time histogram"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:52.605277Z",
+     "start_time": "2019-03-17T19:16:50.114281Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_kg_hist = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "axe_hist_full = fig_kg_hist.add_subplot(2, 1, 1)\n",
+    "axe_hist_trim = fig_kg_hist.add_subplot(2, 1, 2)\n",
+    "plot_hist(axe_hist_full, data[\"gen_time\"], \"keygen time ({})\".format(gen_unit), log_scale, hist_size_gen_time);\n",
+    "plot_hist(axe_hist_trim, data_trimmed[\"gen_time\"], \"keygen time ({})\".format(gen_unit), log_scale, hist_size_gen_time_trim);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Key export time histogram\n",
+    "*Available only for ECTesterReader and keys generated on cards.*"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:52.610858Z",
+     "start_time": "2019-03-17T19:16:52.607191Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "if \"export_time\" in data.dtype.names:\n",
+    "    fig_exp_hist = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "    axe_hist_full = fig_exp_hist.add_subplot(2, 1, 1)\n",
+    "    axe_hist_trim = fig_exp_hist.add_subplot(2, 1, 2)\n",
+    "    plot_hist(axe_hist_full, data[\"export_time\"], \"export time ({})\".format(export_unit), log_scale, hist_size_gen_time);\n",
+    "    plot_hist(axe_hist_trim, data_trimmed[\"export_time\"], \"export time ({})\".format(export_unit), log_scale, hist_size_gen_time_trim);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Moving averages of key generation time"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:54.504830Z",
+     "start_time": "2019-03-17T19:16:54.409189Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_avg = plt.figure(figsize=(10.5, 7), dpi=90)\n",
+    "axe_avg = fig_avg.add_subplot(1, 1, 1)\n",
+    "avg_100 = moving_average(data[\"gen_time\"], 100)\n",
+    "avg_1000 = moving_average(data[\"gen_time\"], 1000)\n",
+    "axe_avg.plot(avg_100, label=\"window = 100\")\n",
+    "axe_avg.plot(avg_1000, label=\"window = 1000\")\n",
+    "if low_bound is not None:\n",
+    "    axe_avg.axhline(y=low_bound, alpha=0.7, linestyle=\"dotted\", color=\"green\", label=\"Low trim bound = {}\".format(low_bound))\n",
+    "if high_bound is not None:\n",
+    "    axe_avg.axhline(y=high_bound, alpha=0.7, linestyle=\"dotted\", color=\"orange\", label=\"Hight trim bound = {}\".format(high_bound))\n",
+    "axe_avg.set_ylabel(\"keygen time ({})\".format(gen_unit))\n",
+    "axe_avg.set_xlabel(\"index\")\n",
+    "axe_avg.legend(loc=\"best\")\n",
+    "del avg_100, avg_1000"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key MSB and LSB histograms"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:16:55.155285Z",
+     "start_time": "2019-03-17T19:16:54.508407Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_priv_hists = plt.figure(figsize=(10.5, 8), dpi=90)\n",
+    "priv_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data[\"priv\"])), dtype=np.dtype(\"u1\"))\n",
+    "priv_lsb = np.array(list(map(lambda x: x & 0xff, data[\"priv\"])), dtype=np.dtype(\"u1\"))\n",
+    "axe_msb_s_hist = fig_priv_hists.add_subplot(2, 1, 1)\n",
+    "axe_lsb_s_hist = fig_priv_hists.add_subplot(2, 1, 2)\n",
+    "plot_hist(axe_msb_s_hist, priv_msb, \"private key MSB\", log_scale)\n",
+    "plot_hist(axe_lsb_s_hist, priv_lsb, \"private key LSB\", log_scale)\n",
+    "del priv_msb, priv_lsb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Public key coordinate MSB and LSB histograms"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:17:06.443596Z",
+     "start_time": "2019-03-17T19:17:05.516616Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def _split(xy):\n",
+    "    x = int.from_bytes(xy[1:byte_size + 1], byteorder=\"big\")\n",
+    "    y = int.from_bytes(xy[1 + byte_size:], byteorder=\"big\")\n",
+    "    return (x, y)\n",
+    "\n",
+    "pub_coords = np.array(list(map(_split, data[\"pub\"])), dtype=np.dtype(\"O\"))\n",
+    "xs = pub_coords[...,0]\n",
+    "ys = pub_coords[...,1]\n",
+    "fig_pub_hists = plt.figure(figsize=(10.5, 14), dpi=90)\n",
+    "\n",
+    "def _plot_coord(data, name, offset):\n",
+    "    axe_msb_pub_hist = fig_pub_hists.add_subplot(4, 1, offset)\n",
+    "    axe_lsb_pub_hist = fig_pub_hists.add_subplot(4, 1, offset + 1)\n",
+    "    pub_msb = np.array(list(map(lambda x: x >> (bit_size - 8), data)))\n",
+    "    pub_lsb = np.array(list(map(lambda x: x & 0xff, data)))\n",
+    "    plot_hist(axe_msb_pub_hist, pub_msb, \"{} coordinate MSB\".format(name), log_scale)\n",
+    "    plot_hist(axe_lsb_pub_hist, pub_lsb, \"{} coordinate LSB\".format(name), log_scale)\n",
+    "    del pub_msb, pub_lsb\n",
+    "\n",
+    "_plot_coord(xs, \"X\", 1)\n",
+    "_plot_coord(ys, \"Y\", 3)\n",
+    "\n",
+    "del pub_coords, xs, ys"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Private key bit length vs time heatmap"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2019-03-17T19:25:51.126642Z",
+     "start_time": "2019-03-17T19:25:50.929170Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig_bl = plt.figure(figsize=(10.5, 12), dpi=90)\n",
+    "gs = gridspec.GridSpec(2, 1, height_ratios=[2.5, 1])\n",
+    "axe_bl_heat = fig_bl.add_subplot(gs[0])\n",
+    "axe_bl_hist = fig_bl.add_subplot(gs[1], sharex=axe_bl_heat)\n",
+    "\n",
+    "bl_data = np.array(list(map(lambda x: x.bit_length(), data_trimmed[\"priv\"])), dtype=np.dtype(\"u2\"))\n",
+    "\n",
+    "h, xe, ye = np.histogram2d(bl_data, data_trimmed[\"gen_time\"], bins=[max(bl_data) - min(bl_data), hist_size_gen_time_trim])\n",
+    "im = axe_bl_heat.imshow(h.T, origin=\"low\", cmap=cmap, aspect=\"auto\", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=norm)\n",
+    "axe_bl_heat.set_xlabel(\"private key bit length\")\n",
+    "axe_bl_heat.set_ylabel(\"keygen time ({})\".format(gen_unit))\n",
+    "\n",
+    "plot_hist(axe_bl_hist, bl_data, \"Private key bit length\", log_scale, align=\"right\")\n",
+    "fig_priv_hist.colorbar(im, ax=[axe_bl_heat, axe_bl_hist])\n",
+    "\n",
+    "del bl_data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "@webio": {
+   "lastCommId": "eaf134b2342a4b71afa1ac9334e37e07",
+   "lastKernelId": "95d0d16c-1dde-451e-94cb-8bd0fefb8378"
+  },
+  "hide_input": false,
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.7.2"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/util/plot_gen.py b/util/plot_gen.py
deleted file mode 100755
index 4ee1ddc..0000000
--- a/util/plot_gen.py
+++ /dev/null
@@ -1,186 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: UTF-8 -*-
-#
-# Script for plotting ECTester key generation results.
-#
-# Example usage:
-#
-#     > java -jar ECTesterReader.jar -g 10000 -b 192 -fp -o gen.csv
-#     ...
-#     > ./plot_gen.py gen.csv
-#     ...
-#
-
-import numpy as np
-from scipy.stats import entropy
-import matplotlib.pyplot as plt
-from matplotlib import ticker, colors
-from copy import deepcopy
-import argparse
-
-from utils import hw, moving_average, plot_hist, miller_correction
-
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser(description="Plot results of ECTester key generation timing.")
-    parser.add_argument("-o", "--output", dest="output", type=argparse.FileType("wb"), help="Write image to [file], do not display.", metavar="file")
-    parser.add_argument("--priv", dest="priv", action="store_true", help="Show private key MSB heatmap plot.")
-    parser.add_argument("--entropy", dest="entropy", action="store_true", help="Show estimated entropy of private key MSB conditioned on time of generation.")
-    parser.add_argument("--hist", dest="hist", action="store_true", help="Show keygen time histogram.")
-    parser.add_argument("--export-hist", dest="export_hist", action="store_true", help="Show export time histogram.")
-    parser.add_argument("--avg", dest="avg", action="store_true", help="Show moving average of keygen time.")
-    parser.add_argument("--hw-hist", dest="hw_hist", action="store_true", help="Show Hamming weight heatmap (private key Hamming weight and keygen time).")
-    parser.add_argument("--log", dest="log", action="store_true", help="Use logarithmic scale.")
-    parser.add_argument("--skip-first", dest="skip_first", nargs="?", const=1, type=int, help="Skip first entry, as it's usually a large outlier.")
-    parser.add_argument("-t", "--title", dest="title", type=str, nargs="?", default="", help="What title to give the figure.")
-    parser.add_argument("file", type=str, help="The file to plot(csv).")
-
-    opts = parser.parse_args()
-
-    with open(opts.file, "r") as f:
-        header = f.readline()
-    header_names = header.split(";")
-    if len(header_names) not in (4, 5):
-        print("Bad data?")
-        exit(1)
-
-    plots = [opts.priv, opts.hist, opts.export_hist, opts.avg, opts.hw_hist]
-    n_plots = sum(plots)
-    if n_plots == 0:
-        plots = [True for _ in range(5)]
-        if len(header_names) == 4:
-            n_plots = 4
-            plots[2] = False
-        else:
-            n_plots = 5
-
-
-    if plots[2] and len(header_names) != 5:
-        n_plots = n_plots - 1
-        if n_plots == 0:
-            print("Nothing to plot.")
-            exit(1)
-        plots[2] = False
-
-    hx = lambda x: int(x, 16)
-    if len(header_names) == 4:
-        data = np.genfromtxt(opts.file, delimiter=";", skip_header=1, converters={2: hx, 3: hx}, dtype=np.dtype([("index", "u4"), ("gen_time", "u4"), ("pub", "O"), ("priv", "O")]))
-    else:
-        data = np.genfromtxt(opts.file, delimiter=";", skip_header=1, converters={3: hx, 4: hx}, dtype=np.dtype([("index", "u4"), ("gen_time", "u4"), ("export_time", "u4"), ("pub", "O"), ("priv", "O")]))
-
-    if opts.skip_first:
-        data = data[opts.skip_first:]
-
-
-    gen_time_data = data["gen_time"]
-    export_time_data = None
-    if "export_time" in data.dtype.names:
-        export_time_data = data["export_time"]
-    pub_data = data["pub"]
-    priv_data = data["priv"]
-
-    gen_unit = "ms"
-    if header_names[1].endswith("[nano]"):
-        gen_unit = r"$\mu s$"
-        np.floor_divide(gen_time_data, 1000, out=gen_time_data)
-    export_unit = "ms"
-    if len(header_names) == 5 and header_names[2].endswith("[nano]"):
-        export_unit = r"$\mu s$"
-        np.floor_divide(export_time_data, 1000, out=export_time_data)
-
-    plt.style.use("ggplot")
-    fig = plt.figure()
-    layout_kwargs = {}
-    if opts.title is None:
-        fig.suptitle(opts.file)
-        layout_kwargs["rect"] = [0, 0.02, 1, 0.98]
-    elif opts.title:
-        fig.suptitle(opts.title)
-        layout_kwargs["rect"] = [0, 0.02, 1, 0.98]
-    fig.tight_layout(**layout_kwargs)
-
-    max_gen_time = max(gen_time_data)
-    min_gen_time = min(gen_time_data)
-    bit_size = len(bin(max(priv_data))) - 2
-
-    sorted_data = np.sort(data, order="gen_time")
-
-    cmap = deepcopy(plt.cm.plasma)
-    cmap.set_bad("black")
-
-    norm = colors.Normalize()
-    if opts.log:
-        norm = colors.LogNorm()
-
-    plot_i = 1
-    if plots[0]:
-        axe_private = fig.add_subplot(n_plots, 1, plot_i)
-        priv_msb = np.array(list(map(lambda x: x >> (bit_size - 8), priv_data)), dtype=np.dtype("u1"))
-        max_msb = max(priv_msb)
-        min_msb = min(priv_msb)
-        heatmap, xedges, yedges = np.histogram2d(priv_msb, gen_time_data, bins=[max_msb - min_msb, max_gen_time - min_gen_time])
-        extent = [min_msb, max_msb, yedges[0], yedges[-1]]
-        axe_private.imshow(heatmap.T, extent=extent, aspect="auto", cmap=cmap, origin="low", interpolation="nearest", norm=norm)
-        axe_private.set_xlabel("private key MSB value")
-        axe_private.set_ylabel("keygen time ({})".format(gen_unit))
-        plot_i += 1
-
-    if plots[1]:
-        axe_hist = fig.add_subplot(n_plots, 1, plot_i)
-        plot_hist(axe_hist, gen_time_data, "keygen time ({})".format(gen_unit), opts.log)
-        plot_i += 1
-
-    if plots[2]:
-        axe_hist = fig.add_subplot(n_plots, 1, plot_i)
-        plot_hist(axe_hist, export_time_data, "export time ({})".format(export_unit), opts.log)
-        plot_i += 1
-
-    if plots[3]:
-        axe_avg = fig.add_subplot(n_plots, 1, plot_i)
-        #if len(header_names) == 5:
-        #   axe_other = axe_avg.twinx()
-        #   axe_other.plot(moving_average(export_time_data, 100), color="green", alpha=0.6, label="export, window = 100")
-        #   axe_other.plot(moving_average(export_time_data, 1000), color="yellow", alpha=0.6, label="export, window = 1000")
-        #   axe_other.legend(loc="lower right")
-        axe_avg.plot(moving_average(gen_time_data, 100), label="window = 100")
-        axe_avg.plot(moving_average(gen_time_data, 1000), label="window = 1000")
-        axe_avg.set_ylabel("keygen time ({})".format(gen_unit))
-        axe_avg.set_xlabel("index")
-        axe_avg.legend(loc="best")
-        plot_i += 1
-
-    if plots[4]:
-        axe_priv_hist = fig.add_subplot(n_plots, 1, plot_i)
-        priv_hw = np.array(list(map(hw, priv_data)), dtype=np.dtype("u2"))
-        h, xe, ye = np.histogram2d(priv_hw, gen_time_data, bins=[max(priv_hw) - min(priv_hw), max_gen_time - min_gen_time])
-        im = axe_priv_hist.imshow(h.T, origin="low", cmap=cmap, aspect="auto", extent=[xe[0], xe[-1], ye[0], ye[-1]], norm=norm)
-        axe_priv_hist.axvline(x=bit_size//2, alpha=0.7, linestyle="dotted", color="white", label=str(bit_size//2) + " bits")
-        axe_priv_hist.set_xlabel("private key Hamming weight")
-        axe_priv_hist.set_ylabel("keygen time ({})".format(gen_unit))
-        axe_priv_hist.legend(loc="best")
-        fig.colorbar(im, ax=axe_priv_hist)
-
-    fig.text(0.01, 0.02, "Data size: {}".format(len(gen_time_data)), size="small")
-
-    if opts.entropy:
-        i = 0
-        entropies = {}
-        while i < len(data):
-            time_val = sorted_data["gen_time"][i]
-            j = i
-            msbs = [0 for _ in range(256)]
-            while j < len(data) and sorted_data["gen_time"][j] == time_val:
-                msbs[(sorted_data["priv"][j] >> (bit_size - 8)) & 0xff] += 1
-                j += 1
-            if j - 100 > i:
-                entropies[time_val] = miller_correction(entropy(msbs, base=2), j - i, 256)
-            i = j
-    
-        entropy = sum(entropies.values())/len(entropies)
-        fig.text(0.01, 0.04, "Entropy of privkey MSB(estimated): {:.2f} b".format(entropy), size="small")
-
-    if opts.output is None:
-        plt.show()
-    else:
-        fig.set_size_inches(12, 10)
-        ext = opts.output.name.split(".")[-1]
-        plt.savefig(opts.output, format=ext, dpi=400, bbox_inches='tight')
diff --git a/util/utils.py b/util/utils.py
index 664373f..d6b9aed 100644
--- a/util/utils.py
+++ b/util/utils.py
@@ -16,18 +16,24 @@ def moving_average(a, n) :
     return ret[n - 1:] / n
 
 
-def plot_hist(axes, data, xlabel=None, log=False):
+def plot_hist(axes, data, xlabel=None, log=False, avg=True, median=True, bins=None, **kwargs):
     time_max = max(data)
     time_min = min(data)
     time_avg = np.average(data)
     time_median = np.median(data)
-    axes.hist(data, bins=time_max - time_min + 1, log=log, align="mid")
-    axes.axvline(x=time_avg, alpha=0.7, linestyle="dotted", color="blue", label="avg = {}".format(time_avg))
-    axes.axvline(x=time_median, alpha=0.7, linestyle="dotted", color="green", label="median = {}".format(time_median))
+    if bins is None:
+        bins = time_max - time_min + 1
+    hist = axes.hist(data, bins=bins, log=log, **kwargs)
+    if avg:
+        axes.axvline(x=time_avg, alpha=0.7, linestyle="dotted", color="blue", label="avg = {}".format(time_avg))
+    if median:
+        axes.axvline(x=time_median, alpha=0.7, linestyle="dotted", color="green", label="median = {}".format(time_median))
     axes.set_ylabel("count" + ("\n(log)" if log else ""))
     axes.set_xlabel("time" if xlabel is None else xlabel)
     axes.xaxis.set_major_locator(ticker.MaxNLocator())
-    axes.legend(loc="best")
+    if avg or median:
+        axes.legend(loc="best")
+    return hist
 
 
 def miller_correction(entropy, samples, bins):