Sage demonstrating some of the computations of https://eprint.iacr.org/2022/596

ecip_calculations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This demonstrates some of the computations of https://eprint.iacr.org/2022/596. First, that the logarithmic derivative of the Weil reciprocity equality holds over random principal divisors using both mixed and duplicate challenge points. Then, that by taking linear combinations of divisor proofs, can check an ECIP and combine multiplicities for the same basis points. Example ECIP uses base -3 and no CM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize an elliptic curve\n",
    "p = 115792089237316195423570985008687907853269984665640564039457584007908834671663\n",
    "r = 115792089237316195423570985008687907852837564279074904382605163141518161494337\n",
    "Fp = GF(p)  # Base Field\n",
    "Fr = GF(r)  # Scalar Field\n",
    "A = 0\n",
    "B = 7\n",
    "E = EllipticCurve(GF(p), [A,B])\n",
    "assert(E.cardinality() == r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "K.<x> = Fp[]\n",
    "L.<y> = K[]\n",
    "eqn = y^2 - x^3 - A * x - B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Returns line passing through points, works for all points and returns 1 for O + O = O\n",
    "def line(A, B):\n",
    "    if A == 0 and B == 0:\n",
    "        return 1\n",
    "    else:\n",
    "        [a, b, c] = Matrix([A, B, -(A+B)]).transpose().kernel().basis()[0]\n",
    "        return a*x + b*y + c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Works for A == B but not A == -A, as the line has no slope or intercept\n",
    "def slope_intercept(A, B):\n",
    "    [a, b, c] = Matrix([A, B, -(A+B)]).transpose().kernel().basis()[0]\n",
    "    return (-a/b, -c/b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fails at 0\n",
    "def eval_point(f, P):\n",
    "    (x, y) = P.xy()\n",
    "    return f(x=x, y=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# f(x) + y g(x) -> (f(x), g(x)), should reduce mod eqn first\n",
    "def get_polys(D):\n",
    "    return ( K(D(y=0)), K(D(y=1) - D(y=0)) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Accepts arbitrary list of points, including duplicates and inverses, and constructs function\n",
    "# intersecting exactly those points if they form a principal divisor (i.e. sum to zero).\n",
    "def construct_function(Ps):\n",
    "    # List of intermediate sums/principal divisors, removes 0\n",
    "    xs = [(P, line(P, -P)) for P in Ps if P != 0]\n",
    "\n",
    "    while len(xs) != 1:\n",
    "        assert(sum(P for (P, _) in xs) == 0)\n",
    "        xs2 = []\n",
    "\n",
    "        # Carry extra point forward\n",
    "        if mod(len(xs), 2) == 1:\n",
    "            x0 = xs[0]\n",
    "            xs = xs[1:]\n",
    "        else:\n",
    "            x0 = None\n",
    "\n",
    "        # Combine the functions for all pairs\n",
    "        for n in range(0, floor(len(xs)/2)):\n",
    "            (A, aNum) = xs[2*n]\n",
    "            (B, bNum) = xs[2*n+1]\n",
    "\n",
    "            # Divide out intermediate (P, -P) factors\n",
    "            num = L((aNum * bNum * line(A, B)).mod(eqn))\n",
    "            den = line(A, -A) * line(B, -B)\n",
    "            D = num / K(den)\n",
    "            \n",
    "            # Add new element\n",
    "            xs2.append((A+B, D))\n",
    "        \n",
    "        if x0 != None:\n",
    "            xs2.append(x0)\n",
    "\n",
    "        xs = xs2\n",
    "    \n",
    "    assert(xs[0][0] == 0)\n",
    "    \n",
    "    # Normalize, might fail but negl probability for random points. Must be done for zkps\n",
    "    # although free to use any coefficient\n",
    "    return D / D(x=0, y=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random principal divisor with n points\n",
    "def random_principal(n):\n",
    "    # For general elliptic curves, might want to clear cofactor depending on application\n",
    "    # Works for arbitrary curve groups\n",
    "    Ps = [E.random_element() for _ in range(0, n-1)]\n",
    "    Ps.append(-sum(Ps))\n",
    "    return Ps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random principal divisor with points with given multiplicities\n",
    "def random_principal_mults(ms):\n",
    "    # Need to invert the last multiplicity to find the correct final value\n",
    "    m0 = ms[-1]\n",
    "    m0Inv = ZZ(Fr(m0)^(-1))\n",
    "    \n",
    "    Ps = [E.random_element() for _ in range(0, len(ms)-1)]\n",
    "    Q = -m0Inv * sum(m * P for (m, P) in zip(ms[:-1], Ps))\n",
    "    Ps.append(Q)\n",
    "    \n",
    "    assert(sum(m * P for (m, P) in zip(ms, Ps)) == 0)\n",
    "    return sum(( m * [P] for (m, P) in zip(ms, Ps) ), [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test at a random principal divisor\n",
    "Ps = random_principal(33)\n",
    "D = construct_function(Ps)\n",
    "(f, g) = get_polys(D)\n",
    "\n",
    "# Should be the same up to constant\n",
    "assert((f^2 - (x^3 + A * x + B) * g^2) / product(x - P.xy()[0] for P in Ps) in Fp)\n",
    "assert(all(eval_point(D, P) == 0 for P in Ps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test at random principal divisor with multiplicity. For a divisor that does not contain\n",
    "# both P and -P for any P, it is sufficient to check the previous conditions and that \n",
    "# gcd(f, g) = 1\n",
    "Ps = random_principal_mults([1,2,3,4,5,6])\n",
    "D = construct_function(Ps)\n",
    "(f, g) = get_polys(D)\n",
    "\n",
    "assert((f^2 - (x^3 + A * x + B) * g^2) / product(x - P.xy()[0] for P in Ps) in Fp)\n",
    "assert(all(eval_point(D, P) == 0 for P in Ps))\n",
    "assert(gcd(f, g) == 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test to check that a function hits exactly a certain set of points uses Weil reciprocity to check that the product of one function over the points of the divisor of the other is the same quantity, up to leading coefficients. Taking the logarithmic derivative wrt a coordinate of one divisor gives a sum of rational functions. That is what is being checked here. While the proof will evaluate the dlog function of at the points, it is important to note that this is also a rational function in the coefficients of the other function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Return logarithmic derivative wrt x\n",
    "def dlog(D):\n",
    "    # Derivative via partials\n",
    "    Dx = D.differentiate(x)\n",
    "    Dy = D.differentiate(y)\n",
    "    Dz = Dx + Dy * ((3*x^2 + A) / (2*y))\n",
    "    \n",
    "    # This is necessary because Sage fails when diving by D\n",
    "    U = L(2*y * Dz)\n",
    "    V = L(2*y * D)\n",
    "\n",
    "    Den = K((V * V(y=-y)).mod(eqn))\n",
    "    Num = L((U * V(y=-y)).mod(eqn))\n",
    "    \n",
    "    # Must clear the denonimator so mod(eqn) well defined\n",
    "    assert(L(y * (Num * D - Den * Dz)).mod(eqn) == 0)\n",
    "    \n",
    "    return Num/Den # == Dz/D"

   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Given a pair of distinct challenge points/line evaluate the function field element\n",
    "def eval_function_challenge_mixed(A0, A1, D):\n",
    "    assert(A0 != A1)\n",
    "    A2 = -(A0 + A1)\n",
    "    (m, b) = slope_intercept(A0, A1)\n",
    "    DLog = dlog(D)\n",
    "    \n",
    "    # Coefficient per point\n",
    "    coeff = 1/((3 * x^2 + A) / (2 * y) - m)\n",
    "    expr = DLog * coeff\n",
    "    \n",
    "    # From paper, check that expr sum is 0, equals slope derivative wrt intercept\n",
    "    assert(sum(eval_point(coeff, P) for P in [A0, A1, A2]) == 0)\n",
    "    \n",
    "    # Evaluate\n",
    "    return sum(eval_point(expr, P) for P in [A0, A1, A2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Given a duplicated challenge point/line evaluate the function field element\n",
    "def eval_function_challenge_dupl(A0, D):\n",
    "    A2 = -(2*A0)\n",
    "    (m, b) = slope_intercept(A0, A2)\n",
    "    DLog = dlog(D)\n",
    "    \n",
    "    # Coefficient for A2\n",
    "    (xA0, yA0) = A0.xy()\n",
    "    (xA2, yA2) = A2.xy()\n",
    "    coeff2 = (2 * yA2) * (xA0 - xA2) / (3 * xA2^2 + A - 2 * m * yA2)\n",
    "    coeff0 = (coeff2 + 2 * m)\n",
    "    \n",
    "    return eval_point(DLog * coeff0, A0) - eval_point(DLog * coeff2, A2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Given a pair of challenge points, detect if duplicate/mixed and modify numerator\n",
    "def eval_point_challenge(A0, A1, P, mult=1):\n",
    "    (m, b) = slope_intercept(A0, A1)\n",
    "    (xP, yP) = P.xy()\n",
    "    \n",
    "    if A0 == A1:\n",
    "        (xA, _) = A0.xy()\n",
    "        num = (xA - xP)\n",
    "    else:\n",
    "        num = -1\n",
    "    \n",
    "    den = yP - m * xP - b\n",
    "    return mult*num/den"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Both should be true (uses same points as higher mult test)\n",
    "[A0, A1] = [E.random_element() for _ in range(0, 2)]\n",
    "assert(eval_function_challenge_mixed(A0, A1, D) == sum(eval_point_challenge(A0, A1, P) for P in Ps))\n",
    "assert(eval_function_challenge_dupl(A0, D) == sum(eval_point_challenge(A0, A0, P) for P in Ps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ECIP takes advantage of the linearity of the right hand sides of the equations to sum multiplicities of the same point in different functions. The following shows how this works in base $-3$ with scalars that are half the length of the field. Note this is important; if the scalars can exceed field length protocol can fail to be sound. Also works for random linear combinations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# return base -3 digits from {-1, 0, 1} from starting with least signficant\n",
    "def base_neg3(n,k):\n",
    "    ds = []\n",
    "    for i in range(0, k):\n",
    "        q = -floor(n/3)\n",
    "        r = ZZ(mod(n, 3))\n",
    "        if r == 2:\n",
    "            q = q - 1\n",
    "            r = -1\n",
    "        ds.append(r)\n",
    "        n = q\n",
    "    \n",
    "    assert(n == 0)\n",
    "    assert(sum(d * (-3)^i for (i, d) in enumerate(ds)))\n",
    "    \n",
    "    return ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# P and -P are counted separately in basis\n",
    "def pos_neg_mults(ds):\n",
    "    a = sum((-3)^i for (i, d) in enumerate(ds) if d == 1)\n",
    "    b = sum((-3)^i for (i, d) in enumerate(ds) if d == -1)\n",
    "    return (a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the principal divisor for each row given sum from previous row\n",
    "def row_function(ds, Ps, Q):\n",
    "    # Construct divisor for row\n",
    "    Q2 = -3*Q + sum(d * P for (d, P) in zip(ds, Ps))\n",
    "    div_ = [-Q, -Q, -Q, -Q2] + [d * P for (d, P) in zip(ds, Ps)]\n",
    "    div = [P for P in div_ if P != 0]\n",
    "    assert(sum(div) == 0)\n",
    "    \n",
    "    # Check that polynomial for row is correct\n",
    "    D = construct_function(div)\n",
    "    LHS = eval_function_challenge_dupl(A0, D)\n",
    "    RHS = sum(eval_point_challenge(A0, A0, P) for P in div)\n",
    "    assert(LHS == RHS)\n",
    "    \n",
    "    return (D, Q2, div)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the function for each row using Shamir's trick and -3\n",
    "def ecip_functions(dss):\n",
    "    rows = list(dss)\n",
    "    rows.reverse()\n",
    "    \n",
    "    Q = 0\n",
    "    Ds = []\n",
    "    for ds in rows:\n",
    "        (p, Q, _) = row_function(ds, Bs, Q)\n",
    "        Ds.append(p)\n",
    "    \n",
    "    # Want lowest order first\n",
    "    Ds.reverse()\n",
    "    return (Q, Ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct digit vectors, note scalars are smaller than characteristic by construction\n",
    "Bs = [E.random_element() for _ in range(0, 20)]                      # Basis Points\n",
    "es = [ZZ.random_element(-2^127, 2^127) for _ in range(0, len(Bs))]   # Linear combination\n",
    "dss_ = [base_neg3(e, 81) for e in es]                                # Base -3 digits\n",
    "epns = list(map(pos_neg_mults, dss_))                                # list of P and -P mults per e\n",
    "dss = Matrix(dss_).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Kinda slow\n",
    "(Q, Ds) = ecip_functions(dss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Q is the final sum\n",
    "assert(Q == sum(e * B for (e, B) in zip(es, Bs)))\n",
    "assert(Q == sum((ep - en) * B for ((ep, en), B) in zip(epns, Bs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Takes two mults and evaluates both P and -P\n",
    "def eval_point_challenge_signed(A0, A1, P, mp, mn):\n",
    "    return eval_point_challenge(A0, A1, P, mult=mp) + eval_point_challenge(A0, A1, -P, mult=mn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sides should equal, remember to account for result point (-Q)\n",
    "LHS = sum((-3)^i * eval_function_challenge_dupl(A0, D) for (i, D) in enumerate(Ds))\n",
    "basisSum = sum(eval_point_challenge_signed(A0, A0, B, ep, en) for ((ep, en), B) in zip(epns, Bs))\n",
    "RHS = basisSum + eval_point_challenge(A0, A0, -Q)\n",
    "assert(LHS == RHS)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.4",
   "language": "sage",
   "name": "sagemath-9.4"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}