{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a code for point counting on genus 3 curves of type y^2 = x^7 + a*x^4 + bx over prime field Fp\n",
    "#\n",
    "# Tested on SageMath v9.0 with Python 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Auxiliary functions\n",
    "\n",
    "R1.<T> = ZZ['T']\n",
    "R2.<a1,a2,a3> = ZZ[\"a1,a2,a3\"]\n",
    "\n",
    "def EC(f):\n",
    "    l_a2 = f.coefficients(sparse=False)[2]\n",
    "    l_a4 = f.coefficients(sparse=False)[1]\n",
    "    l_a6 = f.coefficients(sparse=False)[0]\n",
    "    try:\n",
    "        res = EllipticCurve([0,l_a2,0,l_a4,l_a6])\n",
    "    except ArithmeticError:\n",
    "        return 'singular'\n",
    "    return res\n",
    "\n",
    "def HEC(f):\n",
    "    try:\n",
    "        res = HyperellipticCurve(f)\n",
    "    except ValueError:\n",
    "        return 'singular'\n",
    "    return res\n",
    "\n",
    "def is_nth_root(number, N):\n",
    "    try:\n",
    "        number.nth_root(N)\n",
    "        return True\n",
    "    except ValueError:\n",
    "        return False\n",
    "\n",
    "# Generation of a random point on a hyperelliptic curve C\n",
    "def HEC_random_point(C):\n",
    "    f = C.hyperelliptic_polynomials()[0]\n",
    "    while True:\n",
    "        x_r = f.base_ring().random_element()\n",
    "        y2 = f(x_r)\n",
    "        if y2.is_square():\n",
    "            return [x_r, y2.sqrt()]\n",
    "\n",
    "# Generation of n random unique points on a hyperelliptic curve C\n",
    "def HEC_random_points_uniq(C, n):\n",
    "    f = C.hyperelliptic_polynomials()[0]\n",
    "    res = []\n",
    "    for i in range(1,n+1):\n",
    "        tries = 100\n",
    "        found = False\n",
    "        while tries > 0 and (not found):\n",
    "            P = HEC_random_point(C)\n",
    "            if not (P in res):\n",
    "                res.append(P)\n",
    "                found = True\n",
    "            tries = tries - 1\n",
    "    return res\n",
    "\n",
    "# Generation of a random element in Jacobian of a curve C\n",
    "def JC_random_element(C):\n",
    "    f = C.hyperelliptic_polynomials()[0]\n",
    "    R.<x> = f.base_ring()['x']\n",
    "    J = C.jacobian()\n",
    "    points = HEC_random_points_uniq(C, C.genus())\n",
    "    u = 1\n",
    "    for point in points:\n",
    "        u = u * (x-point[0])\n",
    "    v = f.parent().lagrange_polynomial(points)\n",
    "    return J([u, v])\n",
    "\n",
    "# Computation of the L-polynomial L_{A,q^2}(T) from L_{A,q}(T) for an abelian surface A\n",
    "def g2_L2_from_L(coeffs, q):\n",
    "    a11 = coeffs[0]\n",
    "    a21 = coeffs[1]\n",
    "    a12 = 2*a21 - a11^2\n",
    "    a22 = a21^2 - 4*q*a21 + 2*q^2 + 2*q*a12\n",
    "    return [a12, a22]\n",
    "\n",
    "# Testing of the Jacobian order candidates\n",
    "def check_order(C, N):\n",
    "    f = C.hyperelliptic_polynomials()[0]\n",
    "    q = f.base_ring().cardinality()\n",
    "    g = C.genus()\n",
    "    # Check Hasse-Weil bound\n",
    "    if not ( (q.sqrt() - 1)^(2*g) <= N and N <= (q.sqrt() + 1)^(2*g) ):\n",
    "        return False\n",
    "    res = True\n",
    "    # Check by multiplying on random elements of the Jacobian\n",
    "    i = 1\n",
    "    while i<10 and res:\n",
    "        P = JC_random_element(C)\n",
    "        if list(P*N)!=list(C.jacobian()(0)):\n",
    "            res = False\n",
    "            break\n",
    "        i = i + 1\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computation of the list of possible characteristic polynomials chi_{A,p}(T) from chi_{A,p}(T) (mod p),\n",
    "# where A is an abelian surface in the decomposition of the Jacobian of curve C: y^2 = x^7 + a*x^4 + b x\n",
    "#\n",
    "# chi_{C,p}(T) (mod p) = (T^2 - t2*T + p) * (T^4 + b1 * T^3 + b2 * T^2 + b1*p + p^2)\n",
    "# chi_{A,p}(T) (mod p) = T^4 + b1t * T^3 + b2t * T^2\n",
    "def g2_chi_from_mod_p(C, t2, b1t, b2t, p):\n",
    "    f = C.hyperelliptic_polynomials()[0]\n",
    "    res = []\n",
    "    b1t = lift(GF(p)(b1t))\n",
    "    b2t = lift(GF(p)(b2t))\n",
    "\n",
    "    l1 = ceil(-4/sqrt(p) - b1t/p)\n",
    "    l2 = floor(4/sqrt(p) - b1t/p)\n",
    "    for k1 in range(l1, l2+1):\n",
    "        b1 = b1t + k1 * p\n",
    "        if p < 64:\n",
    "            m1 = ceil(-6 - b2t/p)\n",
    "            m2 = floor(6 - b2t/p)\n",
    "        else:\n",
    "            m1 = ceil(2*abs(b1) / sqrt(p) - 2 - b2t / p)\n",
    "            m2 = floor(b1^2/(4*p)+ 2 + b2t/p)\n",
    "        for k2 in range(m1, m2+1):\n",
    "            b2 = b2t + k2 * p\n",
    "            N = (1 - t2 + p)*(1 + p^2 + b1*(p+1) + b2)\n",
    "            if check_order(C, N):\n",
    "                res = res + [T^4 + b1*T^3 + b2*T^2 + b1*p*T + p^2]\n",
    "                #print(\"---> found chi_A(T) =\", T^4 + b1*T^3 + b2*T^2 + b1*p*T + p^2)\n",
    "    return list(set(res))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computation of the characteristic polynomial chi_{A,p}(T) from chi_{A,p^2}(T) (mod p),\n",
    "# where A is an abelian surface in the decomposition of the Jacobian of curve C: y^2 = x^7 + a*x^4 + b x\n",
    "#\n",
    "# chi_{C,p}(T) (mod p) = (T^2 - t2*T + p) * (T^4 + b1 * T^3 + b2 * T^2 + b1*p + p^2)\n",
    "# chi_{A,p^2}(T) (mod p) = T^4 + b12t * T^3 + b22t * T^2\n",
    "\n",
    "def g2_chi_from_mod_p_2(C, t2, b12t, b22t, p):\n",
    "    f = C.hyperelliptic_polynomials()[0]\n",
    "    b2t = b22t.sqrt(extend=False)\n",
    "    frobs = []\n",
    "    if (2*b2t - b12t).is_square():\n",
    "        b1t = (2*b2t - b12t).sqrt(extend=False)\n",
    "        frobs = g2_chi_from_mod_p(C, t2, b1t, b2t, p)\n",
    "        if b1t != 0:\n",
    "            frobs = frobs + g2_chi_from_mod_p(C, t2, -b1t, b2t, p)\n",
    "    b2t = -b2t    \n",
    "    if b2t != 0 and (2*b2t - b12t).is_square():\n",
    "        b1t = (2*b2t - b12t).sqrt(extend=False)\n",
    "        frobs = frobs + g2_chi_from_mod_p(C, t2, b1t, b2t, p)\n",
    "        if b1t != 0:\n",
    "            frobs = frobs + g2_chi_from_mod_p(C, t2, -b1t, b2t, p)\n",
    "    return list(set(frobs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computation of the characteristic polynomial chi_{C,p}(T) of a curve C: y^2 = x^7 + a*x^4 + b x\n",
    "\n",
    "def g3_frobenius_polynomial(a, b):\n",
    "    Fq = a.parent()\n",
    "    q = Fq.cardinality()\n",
    "    p = Fq.characteristic()\n",
    "    Fq2 = Fq.extension(2)\n",
    "    R.<x> = Fq['x']\n",
    "    f = x^7 + a*x^4 + b*x\n",
    "    C = HEC(f)\n",
    "    #print(\"-> p (mod 3) =\", Mod(p, 3))\n",
    "    #print(\"-> b is a square root:\", is_nth_root(b, 2))\n",
    "    #print(\"-> b is a 3-root:\", is_nth_root(b, 3))\n",
    "    #print(\"-> -1 is a square:\", is_nth_root(Fq(-1), 2))\n",
    "    \n",
    "    frobs = []\n",
    "    E2 = EC(x^3 + a*x^2 + b*x)\n",
    "    t2 = q + 1 - E2.cardinality()\n",
    "    #print(\"-> t2 =\", t2)\n",
    "    \n",
    "    if is_nth_root(b, 2):\n",
    "        b_sq = b.sqrt()\n",
    "        c = -a / (2*b_sq)\n",
    "        E6 = EC(x^3 - 3*x + 2*c)\n",
    "        t6 = q + 1 - E6.cardinality()\n",
    "        p6 = legendre_symbol(3,p) * t6\n",
    "        if Mod(p,3) == 1:\n",
    "            b6 = b_sq^((p-1)/6)\n",
    "            b1A = -(b6^5+b6)*p6\n",
    "            b2A = b6^6*p6^2\n",
    "            frobs = frobs +  [(T^2-t2*T+p)*fr for fr in g2_chi_from_mod_p(C, t2, b1A, b2A, p)]\n",
    "        if Mod(p,3) == 2:\n",
    "            b1A = 0\n",
    "            b2A = -p6^2\n",
    "            frobs = frobs + [(T^2-t2*T+p)*fr for fr in g2_chi_from_mod_p(C, t2, b1A, b2A, p)]\n",
    "    else:\n",
    "        b_sq = Fq2(b).sqrt()\n",
    "        c = -a / (2*b_sq)\n",
    "        x2 = x.change_ring(Fq2)\n",
    "        E6 = EC(x2^3 - 3*x2 + 2*c)\n",
    "        t6 = p^2 + 1 - E6.cardinality()\n",
    "        #print(\"-> E6:\", E6)\n",
    "        #print(\"--> t6 =\", t6)\n",
    "        p6 = legendre_symbol(3,p) * t6\n",
    "        b6 = b_sq^((p^2-1)/6)\n",
    "        b1A = Fq(-(b6^5+b6)*p6)\n",
    "        b2A = Fq(p6^2)\n",
    "        frobs = frobs + [(T^2-t2*T+p)*fr for fr in g2_chi_from_mod_p_2(C, t2, b1A, b2A, p)]\n",
    "        frobs = frobs + [(T^2-t2*T+p)*fr for fr in g2_chi_from_mod_p_2(C, t2, -b1A, b2A, p)]\n",
    "        c = -c\n",
    "        E6 = EC(x2^3 - 3*x2 + 2*c)\n",
    "        t6_2 = p^2 + 1 - E6.cardinality()\n",
    "        if t6_2 != t6:\n",
    "            #print(\"-> E6:\", E6)\n",
    "            #print(\"--> t6 =\", t6)\n",
    "            p6 = legendre_symbol(3,p) * t6\n",
    "            b6 = b_sq^((p^2-1)/6)\n",
    "            b1A = Fq(-(b6^5+b6)*p6)\n",
    "            b2A = Fq(p6^2)\n",
    "            frobs = frobs + [(T^2-t2*T+p)*fr for fr in g2_chi_from_mod_p_2(C, t2, b1A, b2A, p)]\n",
    "            frobs = frobs + [(T^2-t2*T+p)*fr for fr in g2_chi_from_mod_p_2(C, t2, -b1A, b2A, p)]\n",
    "    return list(set(frobs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Tests on small prime fields\n",
    "primes = [5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, \n",
    "          73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, \n",
    "          151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,\n",
    "          229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,\n",
    "          311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,\n",
    "          397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,\n",
    "          479, 487, 491, 499, 503, 509, 521, 523, 541]\n",
    "g = 3\n",
    "\n",
    "for p in primes:\n",
    "    Fq = FiniteField(p)\n",
    "    R.<x> = Fq['x']\n",
    "    ZZx.<T> = ZZ['T']\n",
    "    for i,a in enumerate(Fq):\n",
    "        for j,b in enumerate(Fq):\n",
    "            if a == 0 or b == 0:\n",
    "                continue\n",
    "            f = x^7 + a *x^4 + b*x\n",
    "            C = HEC(f)\n",
    "            if C == 'singular':\n",
    "                continue\n",
    "            print(C)\n",
    "            frob_C_check = C.frobenius_polynomial().subs({x: T})\n",
    "            print(\"-> Frobenius polynomial of C over Fq (by Sage): \", frob_C_check)\n",
    "            print(\"-> Frobenius polynomial of C over Fq (by Sage, fact.): \", frob_C_check.change_ring(QQ).factor())\n",
    "            print(\"-> Frobenius polynomial of C over Fq (by Sage, mod p): \", frob_C_check.change_ring(GF(p)))\n",
    "            frobs = g3_frobenius_polynomial(a,b)\n",
    "            #print(\"-> Frobenius polynomials of C over Fq (by alg.): \", frobs)\n",
    "            frob_C = frobs[0]\n",
    "            print(\"-> Frobenius polynomial of C over Fq (by alg.): \", frob_C)            \n",
    "            print(\"-> Frobenius polynomial of C over Fq (by alg., fact.): \", frob_C.change_ring(QQ).factor())\n",
    "            print(\"-> Size of the list:\", len(frobs))\n",
    "\n",
    "            if frob_C_check in frobs:\n",
    "                print(\"-> [ok]\")\n",
    "            else:\n",
    "                print(\"-> [error]\")\n",
    "            print(\"--------------------\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Example with big fields\n",
    "for i in range(1,100):\n",
    "    p = random_prime(2^130, lbound=2^128)\n",
    "    Fq = FiniteField(p)\n",
    "    R.<x> = Fq['x']\n",
    "    a = Fq.random_element()\n",
    "    b = Fq.random_element()\n",
    "    frobs = g3_frobenius_polynomial(a,b)\n",
    "    print(\"p =\", p)\n",
    "    print(\"p (mod 3) =\", p % 3)\n",
    "    print(\"a =\", a)\n",
    "    print(\"b =\", b)\n",
    "    print(\"b is a cubic residue:\", is_nth_root(b,3))\n",
    "    print(\"Frobenius polynomial:\", frobs[0])\n",
    "    print(\"List size:\", len(frobs))\n",
    "    print(\"Frobenius polynomial (fact.):\", frobs[0].factor())\n",
    "    print(\"-----------------------------------------------\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.0",
   "language": "sage",
   "name": "sagemath"
  },
  "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": 2
}