RCTF 2020 - infantECC

0writeup.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RCTF 2020 - infantECC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the challenge code:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```py\n",
    "from Crypto.Util.number import getStrongPrime, bytes_to_long, long_to_bytes\n",
    "from hashlib import sha256\n",
    "\n",
    "flag = open(\"flag.txt\",\"rb\").read()\n",
    "\n",
    "p=getStrongPrime(512)\n",
    "q=getStrongPrime(512)\n",
    "R=Zmod(p*q)\n",
    "\n",
    "Mx=R.random_element()\n",
    "My=R.random_element()\n",
    "b=My^2-Mx^3\n",
    "E=EllipticCurve(R, [0,b])\n",
    "Ep=EllipticCurve(GF(p), [0,b])\n",
    "Eq=EllipticCurve(GF(q), [0,b])\n",
    "Ecard=Ep.cardinality()*Eq.cardinality()\n",
    "r=random_prime((p^^q)>>100)\n",
    "s=inverse_mod(r, Ecard)\n",
    "\n",
    "print((s,b))\n",
    "print(s*E(Mx,My))\n",
    "print(randint(0,Ecard)*E(Mx,My))\n",
    "print(r^^(bytes_to_long(sha256(long_to_bytes(Mx)).digest())^^bytes_to_long(flag))<<256)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are given a curve over $\\mathbb{Z}_n$ with $n=pq$ an RSA-like modulus. First, note the special form of the curve: $y^2 \\equiv x^3 + b \\pmod{n}$. Whenever $p \\equiv 2$, the curve $y^2 \\equiv x^3 + b \\pmod{p}$ has $p+1$ points for any $b$, such curves are called *supersingular*. Otherwise, the group order varies in $p\\pm \\mathcal{O}(\\sqrt{p})$, which is hard to work with. Thus, we aim to get $p \\equiv q \\equiv 2 \\pmod{3}$, which happens quite often.\n",
    "\n",
    "**Remark:** we are not given $n$, but from $b$ and two points on the curve it is easy to recover."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can immediately dismiss $n\\equiv 2 \\pmod{3}$ since this implies $p\\equiv 1\\pmod{3}$ or $q\\equiv 1\\pmod{3}$. However, we can not easily distinguish the bad case $p\\equiv q \\equiv 1\\pmod{3}$ and the good case $p\\equiv q\\equiv 2\\pmod{3}$, so the attack will only work with probability 1/2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-01T12:43:28.301884Z",
     "iopub.status.busy": "2020-06-01T12:43:28.301287Z",
     "iopub.status.idle": "2020-06-01T12:43:28.307671Z",
     "shell.execute_reply": "2020-06-01T12:43:28.306668Z",
     "shell.execute_reply.started": "2020-06-01T12:43:28.301784Z"
    }
   },
   "source": [
    "In the fortunate case, the cardinality of the curve over $\\mathbb{Z}_n$ is equal to $(p+1)(q+1)=n+p+q+1$. So, we get an equation on the private/public exponent:\n",
    "$$\n",
    "rs \\equiv 1 \\pmod{n+p+q-1},\n",
    "$$\n",
    "$$\n",
    "\\Rightarrow r\\cdot s - k(n+p+q-1)=1, ~~~ k < r < 2^{412}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This equation is very similar to the one in the Boneh-Durfee attack on small secret RSA exponent. However, that attack requires $k < r < n^{0.292}$, while our $r$ is much larger: $r \\approx n^{0.402}$. Do we have any other information?\n",
    "\n",
    "Yes, it is somewhat hidden, but the least 256 bits of $r$ are given out in the last printed value! Can we use it to improve the Boneh-Durfee attack? Unfortunately, not directly, since the attack considers the equation mod $s$, and thus, the value of $r$ does not matter (it does matter indirectly by giving a bound on $k$, but learning bits of $r$ does not change the bound).\n",
    "\n",
    "A possibility is to consider the equation modulo $2^{256}s$. Then we get an equation of the form $c+k(a+z) \\equiv 0\\pmod{2^{256}s}$, where $k<2^{412}, z=p+q<2^{513}$, and the modulus is close to $2^{1280}$. This could be comparable to the Boneh-Durfee bounds: $k<2^{0.282\\cdot1280}=2^{360},z<2^{641}$. However, I did not manage to get it work because of $c\\ne 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative Approach\n",
    "\n",
    "Instead, let's try to tackle with the equation by hands (and LLL).\n",
    "\n",
    "Let's generate a sample instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:43:47.239931Z",
     "iopub.status.busy": "2020-06-02T14:43:47.239428Z",
     "iopub.status.idle": "2020-06-02T14:43:47.466283Z",
     "shell.execute_reply": "2020-06-02T14:43:47.465620Z",
     "shell.execute_reply.started": "2020-06-02T14:43:47.239857Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   s 1020.453\n",
      "   r 404.274\n",
      "rmax 410.041\n"
     ]
    }
   ],
   "source": [
    "from sage.all import log\n",
    "from tqdm import tqdm\n",
    "\n",
    "log2 = lambda val: f\"{float(RR(log(abs(val),2))):.3f}\"\n",
    "proof.arithmetic(False)\n",
    "\n",
    "from random import randint, seed\n",
    "seed(101)\n",
    "\n",
    "p = q = 1\n",
    "while p % 3 != 2:\n",
    "    p = next_prime(randint(1,2**512))\n",
    "while q % 3 != 2:\n",
    "    q = next_prime(randint(1,2**512))\n",
    "\n",
    "n = p*q\n",
    "R = Zmod(p*q)\n",
    "Mx = R.random_element()\n",
    "My = R.random_element()\n",
    "b = My**2 - Mx**3\n",
    "\n",
    "Ep = EllipticCurve(GF(p), [0,b])\n",
    "Eq = EllipticCurve(GF(q), [0,b])\n",
    "E = EllipticCurve(R, [0,b])\n",
    "\n",
    "Ecard = Ep.cardinality()*Eq.cardinality()\n",
    "\n",
    "#r = random_prime((p^^q)>>100)\n",
    "r = next_prime(randint(0, (p^^q)>>100))\n",
    "s = inverse_mod(r, Ecard)\n",
    "\n",
    "print(\"   s\", log2(s))\n",
    "print(\"   r\", log2(r))\n",
    "print(\"rmax\", log2((p^^q)>>100))\n",
    "\n",
    "spt = s*E(Mx, My)\n",
    "randpt = randint(0, Ecard)*E(Mx, My)\n",
    "\n",
    "from hashlib import sha256\n",
    "from Crypto.Util.number import bytes_to_long, long_to_bytes\n",
    "flag = b\"RCTF{Copper5mith_M3thod_f0r_ECC}\"\n",
    "v = r^^(bytes_to_long(sha256(long_to_bytes(Mx)).digest())^^bytes_to_long(flag))<<256"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:43:47.467304Z",
     "iopub.status.busy": "2020-06-02T14:43:47.466942Z",
     "iopub.status.idle": "2020-06-02T14:43:47.473862Z",
     "shell.execute_reply": "2020-06-02T14:43:47.473116Z",
     "shell.execute_reply.started": "2020-06-02T14:43:47.467237Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k 401.430\n"
     ]
    }
   ],
   "source": [
    "k = (r*s-1)//(n+p+q+1)\n",
    "print(\"k\", log2(k))\n",
    "assert r*s - k*(n+p+q+1) == 1\n",
    "assert (1 + k*(n+1+p+q)) % s == 0\n",
    "\n",
    "r0 = v % 2**256  # given\n",
    "r1 = r >> 256\n",
    "assert (1-r0*s + k*(n+1+p+q)) % 2**256 == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $r = 2^{256}r_1 + r_0$, where $r_0 < 2^{256}$. We can rewrite the main equation modulo $2^{256}s$:\n",
    "$$\n",
    "k(n+1)+k(p+q) \\equiv sr_0-1 \\pmod{2^{256}s}.\n",
    "$$\n",
    "Let's try to find a multiple of this equation that will make the left part less than the modulus, or overflow it by a small amount."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:43:47.474940Z",
     "iopub.status.busy": "2020-06-02T14:43:47.474601Z",
     "iopub.status.idle": "2020-06-02T14:43:47.482299Z",
     "shell.execute_reply": "2020-06-02T14:43:47.481629Z",
     "shell.execute_reply.started": "2020-06-02T14:43:47.474897Z"
    }
   },
   "outputs": [],
   "source": [
    "mod = s*2**256\n",
    "assert k*(n+1+p+q) % mod == (s*r0-1) % mod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:43:47.483681Z",
     "iopub.status.busy": "2020-06-02T14:43:47.483314Z",
     "iopub.status.idle": "2020-06-02T14:43:47.556856Z",
     "shell.execute_reply": "2020-06-02T14:43:47.556295Z",
     "shell.execute_reply.started": "2020-06-02T14:43:47.483616Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t 381.838\n",
      "tn 895.010\n",
      "overflow 20.756 1771398\n"
     ]
    }
   ],
   "source": [
    "weight = 2**512\n",
    "m = matrix(ZZ, [\n",
    "    [n+1, weight*1],\n",
    "    [mod, 0]\n",
    "])\n",
    "\n",
    "for tn, t in m.LLL():\n",
    "    if t * tn < 0: continue    \n",
    "    t //= weight\n",
    "    if t < 0:\n",
    "        tn = -tn\n",
    "        t = -t\n",
    "    break\n",
    "else:\n",
    "    raise Exception(\"negative case, let's skip\")\n",
    "    \n",
    "print(\"t\", log2(t))\n",
    "print(\"tn\", log2(tn))\n",
    "\n",
    "w = k * ((n + 1 + p + q) * t % mod) // mod\n",
    "print(\"overflow\", log2(w), w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The overflow has rather large variance over problem instances, it is about 20-32 bits.\n",
    "\n",
    "**... a part of writeup with intermediate results is missing ...**\n",
    "\n",
    "Define\n",
    "$$\n",
    "\\begin{align}\n",
    "h &= \\left\\lfloor \\frac{r_0t}{2^{256}} \\right\\rfloor,\\\\\n",
    "u &= \\left\\lfloor \\frac{t(n+1)}{2^{256}s} \\right\\rfloor,\\\\\n",
    "w &= \\left\\lfloor \\frac{(~t(n + 1 + p + q) \\mod{2^{256}s})k}{2^{256}s} \\right\\rfloor.\n",
    "\\end{align}\n",
    "$$\n",
    "Note that $u$ and $h$ are known, and $w$ is unknown but rather small (experimentally, around 20-32 bits). Then with overwhelming probability\n",
    "$$\n",
    "ku + w - r_1t = h.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:43:47.557971Z",
     "iopub.status.busy": "2020-06-02T14:43:47.557608Z",
     "iopub.status.idle": "2020-06-02T14:43:47.562981Z",
     "shell.execute_reply": "2020-06-02T14:43:47.562314Z",
     "shell.execute_reply.started": "2020-06-02T14:43:47.557913Z"
    }
   },
   "outputs": [],
   "source": [
    "u = t * (n+1) // mod\n",
    "h = r0 * t // 2**256\n",
    "w = k * ((n + 1 + p + q) * t % mod) // mod\n",
    "assert h == k*u + w - r1*t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once $w$ is guessed, we can recover $r_1$ modulo $u$:\n",
    "$$\n",
    "r_1 \\equiv \\frac{w-h}{t} \\pmod{u}\n",
    "$$\n",
    "Let $r_1 = r_{11}u + r_{10}$, where $0 \\le r_{10} < u$. Note that $r_{11}$ is of size 20-32 bits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:43:47.564376Z",
     "iopub.status.busy": "2020-06-02T14:43:47.563968Z",
     "iopub.status.idle": "2020-06-02T14:43:47.573136Z",
     "shell.execute_reply": "2020-06-02T14:43:47.571697Z",
     "shell.execute_reply.started": "2020-06-02T14:43:47.564319Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "to guess (bit size):\n",
      "r11 19.593 19.593\n",
      "  w 20.756\n"
     ]
    }
   ],
   "source": [
    "r10 = r1 % u\n",
    "r11 = r1 // u\n",
    "assert (-h + w) * inverse_mod(t, u) % u == r1 % u == r10\n",
    "print(\"to guess (bit size):\")\n",
    "print(\"r11\", log2(r11), log2(k//t))\n",
    "print(\"  w\", log2(w))\n",
    "assert r == r0 + 2**256*u*r11 + 2**256*r10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T12:24:03.859510Z",
     "iopub.status.busy": "2020-06-02T12:24:03.859118Z",
     "iopub.status.idle": "2020-06-02T12:24:03.864164Z",
     "shell.execute_reply": "2020-06-02T12:24:03.863115Z",
     "shell.execute_reply.started": "2020-06-02T12:24:03.859462Z"
    }
   },
   "source": [
    "Guessing both $r_{11}$ and $w$ is too much. However, we can exploit the curve group and mount a BSGS-like attack. We have\n",
    "$$\n",
    "r_1 = r_{11}u + r_{10} = r_{11}u + \\left(\\frac{w-h}{t}\\mod{u}\\right).\n",
    "$$\n",
    "Since $r_0$ is known, we can express full secret $r$:\n",
    "$$\n",
    "r = r_0 + 2^{256}\\left(\\frac{w-h}{t}\\mod{u}\\right) + 2^{256}ur_{11}.\n",
    "$$\n",
    "Let $G$ be an arbitrary point on the curve $E$ (e.g. one of those given). We know that $[rs-1]G=O$, therefore:\n",
    "$$\n",
    "[(r_0+2^{256}ur_{11})s-1]G = -[2^{256}\\left(\\frac{w-h}{t}\\mod{u}\\right)s]G.\n",
    "$$\n",
    "We can precompute the left part for all possible values of $r_{11}$ and store in the table. Then we guess $w$ and compute the right part and check against the table. It's possible to optimize both steps so that 1-2 curve point additions per guess are needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:43:47.575014Z",
     "iopub.status.busy": "2020-06-02T14:43:47.574545Z",
     "iopub.status.idle": "2020-06-02T14:47:02.772083Z",
     "shell.execute_reply": "2020-06-02T14:47:02.771514Z",
     "shell.execute_reply.started": "2020-06-02T14:43:47.574957Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 4194304/4194304 [03:14<00:00, 21537.50it/s]\n"
     ]
    }
   ],
   "source": [
    "G = spt\n",
    "\n",
    "assert (r0 + 2**256*u*r11) * s * G - G == -(2**256*r10) * s * G\n",
    "\n",
    "# reasonable bounds:\n",
    "# high enough probability to solve,\n",
    "# low enough memory req.\n",
    "if 1:  # full search\n",
    "    start_r = 0\n",
    "    start_w = 0\n",
    "    total_r = 2**22\n",
    "    total_w = 2**23\n",
    "else:  # fast check on known secrets\n",
    "    start_r = r11 - 5\n",
    "    start_w = w - 5\n",
    "    total_r = 10\n",
    "    total_w = 10\n",
    "\n",
    "mask = 2**100-1  # hashing points\n",
    "\n",
    "tab = {}\n",
    "acc = (r0*s-1)*G\n",
    "step = (2**256*u*s)*G\n",
    "acc += start_r * step\n",
    "\n",
    "perc = 0\n",
    "for i in tqdm(range(total_r)):\n",
    "    test = int(acc.xy()[0]) & mask\n",
    "    tab[test] = start_r + i \n",
    "    acc += step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:47:02.773384Z",
     "iopub.status.busy": "2020-06-02T14:47:02.772965Z",
     "iopub.status.idle": "2020-06-02T14:49:32.082398Z",
     "shell.execute_reply": "2020-06-02T14:49:32.081888Z",
     "shell.execute_reply.started": "2020-06-02T14:47:02.773336Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 21%|██        | 1771398/8388608 [02:29<09:17, 11875.62it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution: 1771398 790644\n",
      "recovered r 49943357289587115406335857308667372798949001275969321697728163810970501325410259988956553512194719612541025798846577063031\n",
      "b'RCTF{Copper5mith_M3thod_f0r_ECC}'\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "H = -2**256*s*G\n",
    "base = (-h) * inverse_mod(t, u) % u\n",
    "step_e = inverse_mod(t, u) % u\n",
    "cur_e = base\n",
    "\n",
    "cur = cur_e * H\n",
    "step = step_e * H\n",
    "red = u * H\n",
    "\n",
    "cur_e += start_w * step_e\n",
    "cur_e %= u\n",
    "cur = cur_e * H\n",
    "\n",
    "for w in tqdm(range(total_w)):\n",
    "    test = int(cur.xy()[0]) & mask\n",
    "    if test in tab:\n",
    "        r11 = tab[test]\n",
    "        w += start_w\n",
    "        print(\"Solution:\", w, r11)\n",
    "        r10 = (-h + w) * step_e % u\n",
    "        r = r0 + 2**256 * r10 + 2**256*u*r11\n",
    "        Mx = (r * spt).xy()[0]\n",
    "        print(\"recovered r\", r)\n",
    "        m = (v ^^ r) >> 256\n",
    "        m = m ^^ bytes_to_long(sha256(long_to_bytes( int(Mx) )).digest())\n",
    "        print(long_to_bytes(m))\n",
    "        break\n",
    "\n",
    "    # optimized reduction mod u\n",
    "    # track the exponent and reduce if needed\n",
    "    cur_e += step_e\n",
    "    cur += step\n",
    "    if cur_e >= u:\n",
    "        cur_e -= u\n",
    "        cur -= red"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On my laptop the whole attack with the chosen bounds runs in 6 minutes.\n",
    "Note that the bounds hold only with some probability, so multiple (typically around 10) instances have to be attempted (recall that $p,q$ can lead to non-supersingular curves and thus the probabilty is halved further)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "The method I used is rather unusual and I still have many gaps in its understanding. It is not exactly clear how to derive good bounds for $r_{11}$ and $w$ and why they have so large variance in size. A possible explanation is to look at the relations e.g. $k < r < ((p\\oplus q)\\gg 100)$. While $k$ can be close to the maximum, it has higher chances to be much lower: if $p$ and $q$ agree in several most significant bits, the bound is decreased; then if $r$ is small by chance, the bound is decreased further.\n",
    "\n",
    "Also, accidentally I noticed that $\\lfloor k/t\\rfloor = r_{11} = \\lfloor r_1/u \\rfloor = \\lfloor r/{2^{256}u} \\rfloor$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-06-02T14:50:15.452051Z",
     "iopub.status.busy": "2020-06-02T14:50:15.451622Z",
     "iopub.status.idle": "2020-06-02T14:50:15.456224Z",
     "shell.execute_reply": "2020-06-02T14:50:15.455627Z",
     "shell.execute_reply.started": "2020-06-02T14:50:15.452003Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "790644 790644\n"
     ]
    }
   ],
   "source": [
    "print(k//t, r1//u)\n",
    "assert k // t == r11 == r1 // u == r // (2**256*u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is rather surprising, as $r$ an $k$ are unknowns which are related by a quantity dependent on $p+q$, and we can find such $t,u$ to make this relation to hold without factoring $n$. Also, $t$ was chosen rather arbitrarily!\n",
    "\n",
    "Finally, it feels like the problem should be very easy once $w$ is guessed, but I didn't find a good method avoiding use of the curve.\n",
    "\n",
    "PS: The intended solution seems to be using Coppersmith methods, but I haven't seen it yet."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.0",
   "language": "sage",
   "name": "sagemath3"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}

challenge.sage

from Crypto.Util.number import getStrongPrime, bytes_to_long, long_to_bytes
from hashlib import sha256

flag = open("flag.txt","rb").read()

p=getStrongPrime(512)
q=getStrongPrime(512)
R=Zmod(p*q)

Mx=R.random_element()
My=R.random_element()
b=My^2-Mx^3
E=EllipticCurve(R, [0,b])
Ep=EllipticCurve(GF(p), [0,b])
Eq=EllipticCurve(GF(q), [0,b])
Ecard=Ep.cardinality()*Eq.cardinality()
r=random_prime((p^^q)>>100)
s=inverse_mod(r, Ecard)

print((s,b))
print(s*E(Mx,My))
print(randint(0,Ecard)*E(Mx,My))
print(r^^(bytes_to_long(sha256(long_to_bytes(Mx)).digest())^^bytes_to_long(flag))<<256)

A lot of weird stuff going on!

https://nbviewer.jupyter.org/gist/hellman/e6352b9f12c9759ef23dcd87b64ec8c0/0writeup.ipynb

is nicer display