{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sage is not required"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def legendre(a, p):\n",
    "    return pow(a, (p - 1) // 2, p)\n",
    " \n",
    "def tonelli(n, p):\n",
    "    assert legendre(n, p) == 1, \"not a square (mod p)\"\n",
    "    q = p - 1\n",
    "    s = 0\n",
    "    while q % 2 == 0:\n",
    "        q //= 2\n",
    "        s += 1\n",
    "    if s == 1:\n",
    "        return pow(n, (p + 1) // 4, p)\n",
    "    for z in range(2, p):\n",
    "        if p - 1 == legendre(z, p):\n",
    "            break\n",
    "    c = pow(z, q, p)\n",
    "    r = pow(n, (q + 1) // 2, p)\n",
    "    t = pow(n, q, p)\n",
    "    m = s\n",
    "    t2 = 0\n",
    "    while (t - 1) % p != 0:\n",
    "        t2 = (t * t) % p\n",
    "        for i in range(1, m):\n",
    "            if (t2 - 1) % p == 0:\n",
    "                break\n",
    "            t2 = (t2 * t2) % p\n",
    "        b = pow(c, 1 << (m - i - 1), p)\n",
    "        r = (r * b) % p\n",
    "        c = (b * b) % p\n",
    "        t = (t * c) % p\n",
    "        m = i\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prerequisites\n",
    "- Elliptic curves\n",
    "\n",
    "**Reminder**  \n",
    "Weierstrass curve\n",
    "$$y^2 = x^3 + ax + b$$\n",
    "with $4a^3 + 27b^2 \\neq 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Costello and Smith paper](https://eprint.iacr.org/2017/212.pdf) - must read, I followed from here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Definition**\n",
    "> A Montgomery curve over a field $K$ is defined by the equation\n",
    "> $$\\boxed{M_{A,B}:By^{2}=x^{3}+Ax^{2}+x}$$\n",
    ">for 2 parameters $A, B \\in K$ and with $B(A^2 − 4) \\neq 0$. \n",
    "\n",
    "Or in projective coordinates\n",
    "$$BY^2Z = X^3 + AX^2Z + XZ^2 \\subseteq \\mathbb P^2$$\n",
    "With $\\mathcal O = (0: 1: 0)$ being the point at infinity\n",
    "\n",
    "**J-invariant**\n",
    "$$J(M_{A, B}) = \\dfrac {256(A^2 - 3)^3} {A^2 - 4}$$\n",
    "- Notice that the $\\overline K$ isomorphism class is determined only by $A$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_on(P, C):\n",
    "    \"\"\"\n",
    "    Checks if a point P is on the montgomery curve C\n",
    "    \"\"\"\n",
    "    A, B, p = C\n",
    "    x, y = P\n",
    "    return (B* y**2) % p == (x**3 + A * x**2 + x) % p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lift(x, C):\n",
    "    \"\"\"\n",
    "    Lifts a coordinate `x` to a point on the curve `C` if possible\n",
    "    \"\"\"\n",
    "    A, B, p = C\n",
    "    y = tonelli((x**3 + A * x**2 + x) % p, p)\n",
    "    return (x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The group structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Identity**\n",
    ": $\\mathcal{O}$ acts as identity\n",
    "\n",
    "**Negation**\n",
    ": $(x:y:1) \\mapsto (x:-y:1)$\n",
    "\n",
    "**Addition**  \n",
    "Let $P = (x_P, y_P), \\ Q = (x_Q, y_Q)$ be a two points on our curve. Then for $R = P + Q = (x_R, y_R)$ we have\n",
    "- $x_R = B\\lambda ^2 - (x_P + x_Q) - A$\n",
    "- $y_R = (2x_P + x_Q + A)\\lambda - B\\lambda ^3 - y_P\n",
    "= \\lambda(x_P - x_R) - yP$\n",
    "\n",
    "For \n",
    "$\\lambda = \\begin{cases}\n",
    "\\cfrac {y_Q - y_P} {x_Q - x_P} \\text{ if } P \\neq Q, \\ P \\neq -Q \\\\\n",
    "\\cfrac {3x^2_P + 2Ax_P + 1)} {2By_P} \\text{ if } P = Q \n",
    "\\end{cases}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Crypto.Util.number import inverse\n",
    "\n",
    "def montgomery_addition(P1, P2, C):\n",
    "    \"\"\"\n",
    "    Adds 2 points P1, P2 on the curve C\n",
    "    \"\"\"\n",
    "    \n",
    "    assert(is_on(P1, C)), f\"The point {P1} is not on curve {C}\"\n",
    "    assert(is_on(P2, C)), f\"The point {P2} is not on curve {C}\"\n",
    "    \n",
    "    A, B, p = C\n",
    "    x1, y1 = P1\n",
    "    x2, y2 = P2\n",
    "    \n",
    "    #assert(x1 != x2 and y1 != y2), \"P must not be equal to Q\"\n",
    "    lam = ((y2 - y1) * inverse(x2 - x1, p)) % p\n",
    "    x3 = (B * lam ** 2 - A - x1 - x2) % p\n",
    "    y3 = (lam * (x1 - x3) - y1) % p\n",
    "    \n",
    "    return (x3, y3)\n",
    "\n",
    "def montgomery_doubling(P, C):\n",
    "    \"\"\"\n",
    "    Double a point P on the curve C\n",
    "    \"\"\"\n",
    "    \n",
    "    assert(is_on(P, C)), f\"The point {P} is not on curve {C}\"\n",
    "    \n",
    "    A, B, p = C\n",
    "    x, y = P\n",
    "    \n",
    "    lam = ((3*x**2 + 2*A*x + 1) * inverse(2 * B * y, p)) % p\n",
    "    x3 = (B * lam ** 2 - A - 2 * x) % p\n",
    "    y3 = (lam * (x - x3) - y) % p\n",
    "    \n",
    "    return (x3, y3)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fast arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The paper presents **pseudo-addition** and **pseudo-doubling** which perform x-addition and x-doubling using $x(P), x(Q), x(P-Q)$ with the `xADD` and `xDBL` algorithms  \n",
    "- $+$ efficient x-line arithmetic\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Montgomery ladder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is `algorithm 3` from the paper and it is used to compute $P \\mapsto [k]P$ on a curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def montgomery_binary(P, k, C):\n",
    "    \"\"\"\n",
    "    Return k*P on C given P, k, C\n",
    "    \"\"\"\n",
    "    assert(is_on(P, C)), f\"The point {P} is not on curve {C}\"\n",
    "\n",
    "    R0 = P\n",
    "    R1 = montgomery_doubling(P, C)\n",
    "    l = k.bit_length()\n",
    "    assert((k >> (l-1)) & 1 == 1), \"the l-bit k must have the l-1 bit equal to 1\"\n",
    "    for i in range(l-2, -1, -1):\n",
    "        if (k>>i) & 1 == 0:\n",
    "            R0, R1= montgomery_doubling(R0, C), montgomery_addition(R0, R1, C)\n",
    "        else:\n",
    "            R0, R1= montgomery_addition(R0, R1, C), montgomery_doubling(R1, C)\n",
    "    \n",
    "    return R0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1681274011185213290829845720121217409632760205352840302183329577742232042267, 11000771751458358810062840331211779070151159325452818737845305294782515669468)\n"
     ]
    }
   ],
   "source": [
    "p = (1<<255) - 19\n",
    "B = 1\n",
    "A = 486662\n",
    "curve = (A, B, p)\n",
    "x = 4\n",
    "P = lift(x, curve)\n",
    "k = 0xfeed\n",
    "\n",
    "Q = montgomery_binary(P, k, curve)\n",
    "\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correspondence with other models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Short Weierstrass models**  \n",
    "From Montgomery to Weierstrass map\n",
    "    $ M_{A,B}\\rightarrow E_{a,b}$\n",
    "\n",
    "$$(x,y)\\mapsto (t,v)=\\left({\\frac {x}{B}}+{\\frac {A}{3B}},{\\frac {y}{B}}\\right),a={\\frac {3-A^{2}}{3B^{2}}},b={\\frac {2A^{3}-9A}{27B^{3}}}$$\n",
    "\n",
    "The reverse map does not always exist. In fact, in order to transform an Weierstrass form curve into an Montgomery curve we have to satisfy the following conditions:\n",
    "1. $E_{a, b}$ has order divisible by 4\n",
    "1. $z^{3}+az+b=0$ has at least one root $ \\alpha \\in \\mathbb {F}$\n",
    "2. $3\\alpha ^{2}+a$ is a quadratic residue in $ \\mathbb {F}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_on_w(P, C):\n",
    "    \"\"\"\n",
    "    Checks if a point P is on the Weierstrass curve C\n",
    "    \"\"\"\n",
    "    a, b, p = C\n",
    "    x, y = P\n",
    "    return y**2 % p == (x**3 + a * x + b) % p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "def montgomery_to_weierstrass_curve(C):\n",
    "    \"\"\"\n",
    "    Returns the weierstrass form curve from a given Montgomery curve\n",
    "    \"\"\"\n",
    "    A, B, p = C\n",
    "    \n",
    "    a = ((3 - A**2) * inverse(3 * B ** 2, p)) % p\n",
    "    b = ((2 * A**3 - 9*A) * inverse(27 * B**3, p)) % p\n",
    "    \n",
    "    return (a, b, p)\n",
    "\n",
    "def montgomery_to_weierstrass_point(P, CM, CW):\n",
    "    \"\"\"\n",
    "    Maps a point on a montgomery curve CM onto one on the weierstrass curve CW\n",
    "    \"\"\"\n",
    "    \n",
    "    assert(CW == montgomery_to_weierstrass_curve(CM)), f\"The Weierstrass curve {CW} is not mapped from the Montgomery curve {CM}\"\n",
    "    A, B, p = CM\n",
    "    \n",
    "    x, y = P\n",
    "    t = (x * inverse(B, p) + A * inverse(3 * B, p)) % p\n",
    "    v = (y * inverse(B, p)) % p\n",
    "    \n",
    "    Q = (t, v)\n",
    "    assert(is_on_w(Q, CW))\n",
    "    \n",
    "    return Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = (1<<255) - 19\n",
    "B = 1\n",
    "A = 486662\n",
    "curve = (A, B, p)\n",
    "x = 4\n",
    "P = lift(x, curve)\n",
    "k = 0xfeed\n",
    "\n",
    "Q = montgomery_binary(P, k, curve)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(19298681539552699237261830834781317975544997444273427339909597334573241639236, 55751746669818908907645289078257140818241103727901012315294400837956729358436, 57896044618658097711785492504343953926634992332820282019728792003956564819949)\n",
      "(19298681539552699237261830834781317975544997444273427339909597334652188435541, 47499954730490638715091885595963621956489259355261559690379252041373947974816)\n",
      "(20979955550737912528091676554902535385177757649626267642092926912394420477804, 11000771751458358810062840331211779070151159325452818737845305294782515669468)\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "curve_w = montgomery_to_weierstrass_curve(curve)\n",
    "print(curve_w)\n",
    "\n",
    "T = montgomery_to_weierstrass_point(P, curve, curve_w)\n",
    "print(T)\n",
    "\n",
    "S = montgomery_to_weierstrass_point(Q, curve, curve_w)\n",
    "print(S)\n",
    "\n",
    "\n",
    "\n",
    "# Can't be bothered to implement double and add in Weierstrass form\n",
    "from ecdsa.ellipticcurve import CurveFp, Point\n",
    "\n",
    "a, b, p = curve_w\n",
    "CW = CurveFp(p, a, b)\n",
    "T_ = Point(CW, *T)\n",
    "S_ = Point(CW, *S)\n",
    "\n",
    "print(S_ == k * T_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resources\n",
    "- [wikipedia](https://en.wikipedia.org/wiki/Montgomery_curve)\n",
    "- [Bernstein and Lange paper](https://eprint.iacr.org/2017/293.pdf)\n",
    "- [First appeared](https://www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866113-7/S0025-5718-1987-0866113-7.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}