{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "<style>\n",
    "    p {\n",
    "        margin-top: 1em;\n",
    "        margin-bottom: 0.5em;\n",
    "    }\n",
    "</style>\n",
    "For the sake of comparison we first give formulas for Weierstrass curves in projective coordinates; scroll down to see the formulas for Edwards curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "def add_affine(P,Q,A):\n",
    "    \"\"\" addition algorithm for affine points on an elliptic curve in short Weierstrass form y^2=x^3+Ax+B\"\"\"\n",
    "    if P == O: return Q\n",
    "    if Q == O: return P\n",
    "    assert P[2] == 1 and Q[2] == 1  # we are using affine formulas, so make sure points are in affine form\n",
    "    x1=P[0]; y1=P[1]; x2=Q[0]; y2=Q[1];\n",
    "    if x1 != x2:\n",
    "        m = (y2-y1)/(x2-x1)         # usual case: P and Q are distinct and not opposites\n",
    "    else:\n",
    "        if y1 == -y2: return O      # P = -Q (includes case where P=Q is 2-torsion)\n",
    "        m = (3*x1^2+A) / (2*y1)     # P = Q so we are doubling\n",
    "    x3 = m^2-x1-x2\n",
    "    y3 = m*(x1-x3)-y1\n",
    "    return (x3,y3,1)\n",
    "\n",
    "def dbl_affine(P,A):\n",
    "    \"\"\"doubles a projective non-2-torsion point on an elliptic curve in short Weierstrass form\"\"\"\n",
    "    x1=P[0]; y1=P[1];\n",
    "    m = (3*x1^2+A) / (2*y1)\n",
    "    x3 = m^2-2*x1\n",
    "    y3 = m*(x1-x3)-y1\n",
    "    return (x3,y3,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "def add_projective(P,Q,A):\n",
    "    \"\"\"addition algorithm for projective points on an elliptic curve in short Weierstrass form y^2=x^3+Ax+B\"\"\"\n",
    "    if P[2] == 0: return Q\n",
    "    if Q[2] == 0: return P\n",
    "    x1=P[0]; y1=P[1]; z1=P[2]; x2=Q[0]; y2=Q[1]; z2=Q[2]\n",
    "    # We use formulas from http://hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#addition-add-1998-cmo-2\n",
    "    x1z2 = x1*z2; x2z1 = x2*z1; y1z2 = y1*z2; y2z1 = y2*z1\n",
    "    if x1z2 != x2z1:\n",
    "        # usual case: P and Q are distinct and not opposites\n",
    "        z1z2 = z1*z2\n",
    "        u = y2z1-y1z2; uu = u^2\n",
    "        v = x2z1-x1z2; vv = v^2; vvv = v^3\n",
    "        r = vv*x1z2\n",
    "        s = uu*z1z2 - vvv - 2*r\n",
    "        x3 = v*s\n",
    "        y3 = u*(r-s)-vvv*y1z2\n",
    "        z3 = vvv*z1z2\n",
    "    else:\n",
    "        if y1z2 == -y2z1: return (0,1,0)      # P = -Q (includes case where P=Q is 2-torsion)\n",
    "        # P = Q so we are doubling\n",
    "        xx = x1^2; yy = y1^2; zz = z1^2\n",
    "        w = A*zz + 3*xx\n",
    "        s = 2*y1*z1; ss = s^2; sss = s*ss\n",
    "        r = y1*s; rr = r^2\n",
    "        b = (x1+r)^2-xx-rr\n",
    "        h = w^2-2*b\n",
    "        x3 = h*s\n",
    "        y3 = w*(b-h)-2*rr\n",
    "        z3 = sss\n",
    "    return (x3,y3,z3)\n",
    "\n",
    "def negate(P):\n",
    "    \"\"\"negates the point P on an elliptic curve in short Weierstrass form\"\"\"\n",
    "    if P[2] == 0: return (0,1,0)\n",
    "    return (P[0],-P[1],P[2])\n",
    "    \n",
    "def dbl_projective(P,A):\n",
    "    \"\"\"doubles a projective non-2-torsion point on an elliptic curve in short Weierstrass form\"\"\"\n",
    "    x1 = P[0]; y1 = P[1]; z1 = P[2]\n",
    "    xx = x1^2; yy = y1^2; zz = z1^2\n",
    "    w = A*zz + 3*xx\n",
    "    s = 2*y1*z1; ss = s^2; sss = s*ss\n",
    "    r = y1*s; rr = r^2\n",
    "    b = (x1+r)^2-xx-rr\n",
    "    h = w^2-2*b\n",
    "    x3 = h*s\n",
    "    y3 = w*(b-h)-2*rr\n",
    "    z3 = sss\n",
    "    return (x3,y3,z3)\n",
    "\n",
    "def equal_projective(P,Q):\n",
    "    \"\"\"Test whether two projective points are equal or not\"\"\"\n",
    "    return P[0]*Q[2] == Q[0]*P[2] and P[1]*Q[2] == Q[1]*P[2]\n",
    "    \n",
    "def random_point(A,B):\n",
    "    \"\"\"generates a random affine point on the elliptic curve y^2 = x^2 + Ax + B\"\"\"\n",
    "    F=A.parent()\n",
    "    x0=F.random_element(); t=(x0^3+A*x0+B)\n",
    "    while not t.is_square():\n",
    "        x0=F.random_element(); t=(x0^3+A*x0+B)\n",
    "    return (x0,sqrt(t),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Timing add_affine on random points...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.82 µs ± 344 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n",
      "Timing add_add_affine on equal points...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.73 µs ± 321 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n",
      "Testing dbl_affine...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.33 µs ± 280 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n",
      "Timing add_projective on random points...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.75 µs ± 178 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n",
      "Timing add_projective on equal points...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12.6 µs ± 159 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n",
      "Testing dbl_projective...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.13 µs ± 65.7 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n"
     ]
    }
   ],
   "source": [
    "p=random_prime(2^256,lbound=2^255)\n",
    "F=GF(p)\n",
    "A=F(3); B=F.random_element()\n",
    "P0=random_point(A,B)\n",
    "Q0=random_point(A,B)\n",
    "print(\"Timing add_affine on random points...\")\n",
    "%timeit -n 30000 R=add_affine(P0,Q0,A)\n",
    "print(\"Timing add_add_affine on equal points...\")\n",
    "%timeit -n 30000 R=add_affine(P0,P0,A)\n",
    "print(\"Testing dbl_affine...\")\n",
    "%timeit -n 30000 R=dbl_affine(P0,A)\n",
    "print(\"Timing add_projective on random points...\")\n",
    "%timeit -n 30000 R=add_projective(P0,Q0,A)\n",
    "print(\"Timing add_projective on equal points...\")\n",
    "%timeit -n 30000 R=add_projective(P0,P0,A)\n",
    "print(\"Testing dbl_projective...\")\n",
    "%timeit -n 30000 R=dbl_projective(P0,A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Note that distinct triples of proejctive coordinates may define the same projective point.\n",
    "To compare points for equality we need to use the equal_projective function above (which also works for Edwards curves)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R1 = (39631903559207355134813242690485305706327172391023199616984319804673277252824, 33869231864040439687106036648793069871874765452794702250861206505319640701175, 1271603232244282005892045156571757933303549172868729422943298126929250743038)\n",
      "R2 = (39003209298511745779211371247669215406384172665513230424454802014897447644459, 44765880993678661226918577289361451240836579603741727790577915314251084196108, 77363509625474818908132568781582763179407795883667700618495823692641474154245)\n",
      "R1 == R2: False\n",
      "equal_projective(R1,R2): True\n"
     ]
    }
   ],
   "source": [
    "R1=add_projective(P0,Q0,A)\n",
    "R2=add_projective(Q0,P0,A)\n",
    "print(\"R1 = \" + str(R1))\n",
    "print(\"R2 = \" + str(R2))\n",
    "print(\"R1 == R2: \" + str(R1 == R2))\n",
    "print(\"equal_projective(R1,R2): %s\" % equal_projective(R1,R2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "def add_edwards(P,Q,d):\n",
    "    \"\"\"adds two projective points on an elliptic curve in Edwards form (assumes d non-square, so operation is complete)\"\"\"\n",
    "    # uses unoptimized formula derived in class (could save 1M with optimization)\n",
    "    x1=P[0]; y1=P[1]; z1=P[2]; x2=Q[0]; y2=Q[1]; z2=Q[2]\n",
    "    x1y2 = x1*y2; x2y1 = x2*y1\n",
    "    r = z1*z2; rr = r^2; s = x1y2+x2y1; t = d*x1y2*x2y1; u = y1*y2 - x1*x2\n",
    "    v = rr+t; w = rr-t\n",
    "    x3 = r*s*w\n",
    "    y3 = r*u*v\n",
    "    z3 = v*w\n",
    "    return (x3,y3,z3)\n",
    "    \n",
    "def dbl_edwards(P):\n",
    "    \"\"\"doubles a projective point on an elliptic curve in Edwards form\"\"\"\n",
    "    # uses optimized formula from Bernstein-Lange\n",
    "    x1 = P[0]; y1 = P[1]; z1 = P[2]\n",
    "    B=(x1+y1)^2; C=x1^2; D=y1^2; E = C+D;\n",
    "    J = E-2*z1^2; x3 = (B-E)*J; y3 = E*(C-D); z3 = E*J\n",
    "    return (x3,y3,z3)\n",
    "    \n",
    "def neg_edwards(P):\n",
    "    return (-P[0],P[1],P[2])\n",
    "    \n",
    "def random_edwards_point(d):\n",
    "    \"\"\"generates a random projective point on the Edwards curve x^2 + y^2 = 1+d*x^2*y^2 (assumes d non-square)\"\"\"\n",
    "    F=d.parent()\n",
    "    x0=F.random_element(); t=(1-x0^2)/(1-d*x0^2)\n",
    "    while not t.is_square():\n",
    "        x0=F.random_element(); t=(1-x0^2)/(1-d*x0^2)\n",
    "    return (x0,sqrt(t),1)\n",
    "    \n",
    "def on_edwards_curve (P,d):\n",
    "    x0=P[0]; y0=P[1]; z0=P[2]\n",
    "    return z0^2*(x0^2+y0^2) == z0^4 + d*x0^2*y0^2\n",
    "    \n",
    "def is_identity_edwards (P):\n",
    "    if P[0] == 0 and P[1] == P[2]: return true\n",
    "    else: return false"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Affine points on Edwards curve x^2 + y^2 = 1 + 3x^2y^2 over F_17\n",
      "[0, 1]\n",
      "[0, 16]\n",
      "[1, 0]\n",
      "[2, 5]\n",
      "[2, 12]\n",
      "[3, 4]\n",
      "[3, 13]\n",
      "[4, 3]\n",
      "[4, 14]\n",
      "[5, 2]\n",
      "[5, 15]\n",
      "[7, 7]\n",
      "[7, 10]\n",
      "[10, 7]\n",
      "[10, 10]\n",
      "[12, 2]\n",
      "[12, 15]\n",
      "[13, 3]\n",
      "[13, 14]\n",
      "[14, 4]\n",
      "[14, 13]\n",
      "[15, 5]\n",
      "[15, 12]\n",
      "[16, 0]\n"
     ]
    }
   ],
   "source": [
    "p = 17\n",
    "F=GF(p)\n",
    "d = F.multiplicative_generator()\n",
    "print(\"Affine points on Edwards curve x^2 + y^2 = 1 + %dx^2y^2 over F_%d\"%(d,p))\n",
    "for a in F:\n",
    "    for b in F:\n",
    "        if on_edwards_curve([a,b,1],d): print([a,b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Timing add_edwards on random points...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.44 µs ± 366 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n",
      "Timing add_edwards on equal points...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.46 µs ± 296 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n",
      "Timing dbl_edwards...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.96 µs ± 559 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)\n"
     ]
    }
   ],
   "source": [
    "p=random_prime(2^256,lbound=2^255)\n",
    "F=GF(p)\n",
    "d=F(2);\n",
    "while d.is_square(): d+=1;\n",
    "P0=random_edwards_point(d)\n",
    "Q0=random_edwards_point(d)\n",
    "print(\"Timing add_edwards on random points...\")\n",
    "%timeit -n 30000 add_edwards(P0,Q0,d)\n",
    "print(\"Timing add_edwards on equal points...\")\n",
    "%timeit -n 30000 add_edwards(P0,P0,d)\n",
    "print(\"Timing dbl_edwards...\")\n",
    "%timeit -n 30000 dbl_edwards(P0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sagemath",
   "metadata": {
    "cocalc": {
     "description": "Open-source mathematical software system",
     "priority": 1,
     "url": "https://www.sagemath.org/"
    }
   },
   "name": "sage-9.2",
   "resource_dir": "/ext/jupyter/kernels/sage-9.2"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}