{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a code for point counting on the genus 4 curves of type C: y^2 = x^9 + a*x^5 + b*x over F_q.\n",
    "#\n",
    "# The method is described here: https://arxiv.org/abs/1902.05992\n",
    "\n",
    "# Tested on SageMath 9.0 with Python 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Auxilary functions\n",
    "R1.<T> = ZZ['T']\n",
    "R2.<a1,a2,a3,a4> = ZZ[\"a1,a2,a3,a4\"]\n",
    "\n",
    "def EC(f):\n",
    "    try:\n",
    "        res = EllipticCurve(f)\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",
    "def join(list1, list2):\n",
    "    res = list1\n",
    "    for e in list2:\n",
    "        if not (e in list1):\n",
    "            res = res + [e]\n",
    "    return res\n",
    "\n",
    "# Dickson polynomial of degree n\n",
    "def D_pol(n, x, alpha):\n",
    "    FF = alpha.base_ring()\n",
    "    res = 0\n",
    "    for k in range(0, floor(n/2)+1):\n",
    "        res += FF(n/(n-k) * binomial(n-k,k)) * (-alpha)^k * x^(n-2*k)\n",
    "    return res\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",
    "# Test for the Jacobian order candidates\n",
    "def check_order(C, N, points_check=True):\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",
    "    if points_check:\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\n",
    "\n",
    "def find_k(b):\n",
    "    if is_nth_root(b, 8):\n",
    "        res = 1\n",
    "    elif is_nth_root(b, 4):\n",
    "        res = 2\n",
    "    elif is_nth_root(b, 2):\n",
    "        res = 4\n",
    "    else:\n",
    "        res = 8\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function computes the characteristic polynomial of the curve\n",
    "# X_1: y^2 = (x-2*b^(1/8))*(D_4(x, b^(1/4)) + a) over Fq[b^(1/8)]\n",
    "# using a twist defined over Fq[b^(1/2)] to speed up the computations.\n",
    "#\n",
    "def X1_frobenius_polynomial(a, b, precomp_X1t=None):\n",
    "    F = b.parent()\n",
    "    q = F.cardinality()\n",
    "    F2 = F.extension(2)\n",
    "    R.<z> = F['z']\n",
    "    \n",
    "    k = find_k(b)\n",
    " \n",
    "    if k == 1:\n",
    "        alpha_sq = b.nth_root(8)\n",
    "        alpha = alpha_sq^2\n",
    "        X1 = HEC((z - 2*alpha_sq)*(D_pol(4,z,alpha) + a))\n",
    "        return X1.frobenius_polynomial()\n",
    "\n",
    "    Fk = F.extension(k)\n",
    "    \n",
    "    if not b.is_square():\n",
    "        b = F2(b)\n",
    "        R.<z> = F2['z']\n",
    "    b_sq = b.sqrt()\n",
    "    fX1t = (z+2)*(D_pol(4, z, 1) + a/b_sq)\n",
    "    X1t = HEC(fX1t) # Twist\n",
    "    print(\"-> X1t':\", X1t)\n",
    "    frobX1t = precomp_X1t\n",
    "    if frobX1t == None:\n",
    "        print(\"--> Computing Frob. of X1t over Fq[b^(1/2)] ...\")\n",
    "        frobX1t = X1t.frobenius_polynomial()\n",
    "        print(\"---> \", frobX1t)\n",
    "    else:\n",
    "        print(\"--> Frob. of X1t over Fq[b^(1/2)] (precomputed):\", frobX1t)\n",
    "    \n",
    "    b1 = frobX1t.coefficients(sparse=False)[3]\n",
    "    b2 = frobX1t.coefficients(sparse=False)[2]\n",
    "\n",
    "    i = 1\n",
    "    if k == 8:\n",
    "        i = 2\n",
    "    while i < k:\n",
    "        b1,b2 = g2_L2_from_L([b1,b2], q^i)\n",
    "        i = i * 2\n",
    "    if not is_nth_root(Fk(b_sq), 8):\n",
    "        b1 = -b1\n",
    "    frobX1 = T^4 + b1*T^3 + b2*T^2 + b1*q^k*T + q^(2*k)\n",
    "    return frobX1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computation of the characteristic polynomial chi_{C,q}(T) from chi_{C,q^2}(T) for a genus 4 curve.\n",
    "def L_from_L2(C, a12, a22, a32, a42, q):\n",
    "    g = 4\n",
    "    candidates = []\n",
    "    # case of a1 != 0\n",
    "    c0 = (128*q^4-128*a12*q^3+32*a12^2*q^2+128*a32*q-64*a12*a22*q+16*a12^3*q-64*a42+16*a22^2-8*a12^2*a22+a12^4)^2\n",
    "    c2 = -131072*q^7+163840*a12*q^6-32768*a22*q^5-65536*a12^2*q^5-81920*a32*q^4+45056*a12*a22*q^4+5120*a12^3*q^4+65536*a42*q^3+49152*a12*a32*q^3-16384*a22^2*q^3-12288*a12^2*a22*q^3+2048*a12^4*q^3-49152*a12*a42*q^2+4096*a12^2*a32*q^2+8192*a12*a22^2*q^2-5120*a12^3*a22*q^2+768*a12^5*q^2+16384*a22*a42*q-16384*a32^2*q+4096*a12*a22*a32*q-1024*a12^3*a32*q-4096*a22^3*q+4096*a12^2*a22^2*q-1280*a12^4*a22*q+128*a12^6*q+8192*a32*a42-6144*a12*a22*a42+1536*a12^3*a42+2048*a22^2*a32-1024*a12^2*a22*a32+128*a12^4*a32-512*a12*a22^3+384*a12^3*a22^2-96*a12^5*a22+8*a12^7\n",
    "    c4 = 253952*q^6-233472*a12*q^5+47104*a22*q^4+65024*a12^2*q^4+57344*a32*q^3-26624*a12*a22*q^3-5632*a12^3*q^3-61440*a42*q^2-12288*a12*a32*q^2+7168*a22^2*q^2-2048*a12^2*a22*q^2+1344*a12^4*q^2+22528*a12*a42*q-2048*a22*a32*q-2560*a12^2*a32*q+3584*a12*a22^2*q-1280*a12^3*a22*q+96*a12^5*q-7168*a22*a42+1280*a12^2*a42+4096*a32^2-2048*a12*a22*a32+512*a12^3*a32-256*a22^3+576*a12^2*a22^2-240*a12^4*a22+28*a12^6\n",
    "    c6 = -204800*q^5+136192*a12*q^4-16384*a22*q^3-29696*a12^2*q^3-12288*a32*q^2+1024*a12*a22*q^2+4096*a12^3*q^2+20480*a42*q-1024*a12*a32*q+1024*a22^2*q-320*a12^4*q-2560*a12*a42-1024*a22*a32+768*a12^2*a32+384*a12*a22^2-320*a12^3*a22+56*a12^5\n",
    "    c8 = 79104*q^4-38144*a12*q^3+512*a22*q^2+7104*a12^2*q^2+256*a32*q+640*a12*a22*q-800*a12^3*q-2176*a42+512*a12*a32+96*a22^2-240*a12^2*a22+70*a12^4\n",
    "    c10 = -15360*q^3+5376*a12*q^2+256*a22*q-768*a12^2*q+128*a32-96*a12*a22+56*a12^3\n",
    "    c12 = 1472*q^2-352*a12*q-16*a22+28*a12^2\n",
    "    c14 = 8*a12-64*q\n",
    "    eq = a1^16 + c14 * a1^14 + c12 * a1^12 + c10 * a1^10 + c8 * a1^8 + c6 * a1^6 + c4 * a1^4 + c2 * a1^2 + c0\n",
    "    l = next_prime(4*g*ceil(q.sqrt()) + 1)\n",
    "    Fl = FiniteField(l)\n",
    "    eq = eq.change_ring(Fl).univariate_polynomial()\n",
    "    B = factor(eq)\n",
    "    for t,e in B:\n",
    "        if t.degree() > 1:\n",
    "            continue\n",
    "        t0 = t.monic().coefficients(sparse=False)[0]\n",
    "        a1c = lift(Fl(-t0))\n",
    "        if a1c > floor(2*g*q.sqrt()):\n",
    "            a1c = a1c - l\n",
    "        if a1c == 0:\n",
    "            continue\n",
    "        if (a1c^2 + a12) % 2 != 0:\n",
    "            continue\n",
    "        a2c = ZZ((a1c^2 + a12) / 2)\n",
    "        d = (2*a2c+a1c^2)*q^2-2*a1c^2*a2c*q-a32+a2c*a22-a2c^3+a1c^2*a2c^2\n",
    "        if d.is_square():\n",
    "            d_s = d.sqrt()\n",
    "            a4c = -(2*a1c*d_s+2*a1c^2*q-a22+a2c^2-2*a1c^2*a2c)/2\n",
    "            a3c=(2*a4c-a22+a2c^2)/(2*a1c)\n",
    "            if not QQ(a2c).is_integer() or not QQ(a3c).is_integer() or not QQ(a4c).is_integer():\n",
    "                continue\n",
    "            candidates.append([ZZ(a1c), ZZ(a2c), ZZ(a3c), ZZ(a4c)])\n",
    "            a4c = (2*a1c*d_s-2*a1c^2*q+a22-a2c^2+2*a1c^2*a2c)/2\n",
    "            a3c=(2*a4c-a22+a2c^2)/(2*a1c)\n",
    "            if not QQ(a2c).is_integer() or not QQ(a3c).is_integer() or not QQ(a4c).is_integer():\n",
    "                continue\n",
    "            candidates.append([ZZ(a1c), ZZ(a2c), ZZ(a3c), ZZ(a4c)])\n",
    "    # case of a1 = 0\n",
    "    a1c = 0\n",
    "    if a12 % 2 == 0:\n",
    "        a2c = ZZ(a12 / 2)\n",
    "        if (a22-a2c^2) % 2 == 0:\n",
    "            a4c = ZZ((a22-a2c^2)/2)\n",
    "            a3c_s = 2*a2c*q^2+2*a2c*a4c-a32\n",
    "            if is_square(a3c_s):\n",
    "                a3c = ZZ(sqrt(a3c_s))\n",
    "                candidates.append([a1c, a2c, a3c, a4c])\n",
    "                if a3c != 0:\n",
    "                    candidates.append([a1c, a2c, -a3c, a4c])\n",
    "    # exclude extra solutions\n",
    "    res = []\n",
    "    for a1c,a2c,a3c,a4c in candidates:\n",
    "        if [a1c,a2c,a3c,a4c] in res:\n",
    "            continue\n",
    "        N = ZZ(q^4 + a1c*q^3 + a2c*q^2 + a3c*q + a4c + a3c + a2c + a1c + 1)\n",
    "        if check_order(C,N):\n",
    "            res.append([a1c, a2c, a3c, a4c])\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computation of the characteristic polynomial chi_{C,q}(T) from chi_{C,q^k}(T) for a genus 4 curve C and k = 2^m\n",
    "def L_from_Lk(C, a1k, a2k, a3k, a4k, q, k):\n",
    "    print(\"-> Computation of [a1,a2,a3,a4] from [a1k, a2k, a3k, a4k] =\", [a1k, a2k, a3k, a4k])\n",
    "    f = C.hyperelliptic_polynomials()[0]\n",
    "    if k == 1:\n",
    "        res = [[a1k, a2k, a3k, a4k]]\n",
    "        return res\n",
    "    F = f.base_ring()\n",
    "    L = [[a1k, a2k, a3k, a4k]]\n",
    "    i = k\n",
    "    while i != 1:\n",
    "        S = []\n",
    "        for t in L:\n",
    "            F2 = F.extension(int(i/2))\n",
    "            C2 = HEC(f.change_ring(F2))\n",
    "            S = join(S, L_from_L2(C2, t[0], t[1], t[2], t[3], q^(i/2)))\n",
    "        L = S\n",
    "        i = int(i/2)\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computation of the characteristic polynomial of the curve C: y^2 = x^9 + a*x^5 + b*x over a finite field F_q\n",
    "#\n",
    "# This is an implementation of the Algorithm 3 from https://arxiv.org/abs/1902.05992\n",
    "def g4_frobenius_polynomial(a, b, precomp_X1t=None):\n",
    "    Fq = a.parent()\n",
    "    q = Fq.cardinality()\n",
    "    p = Fq.characteristic()\n",
    "    k = find_k(b)\n",
    "\n",
    "    R.<x> = Fq['x']\n",
    "    C = HEC(x^9 + a*x^5 + b*x)\n",
    "    print(C)\n",
    "    print(\"-> p (mod 8) =\", Mod(p,8))\n",
    "    print(\"-> b is a square root:\", is_nth_root(b, 2))\n",
    "    print(\"-> b is a 4-root:\", is_nth_root(b, 4))\n",
    "    print(\"-> b is a 8-root:\", is_nth_root(b, 8))\n",
    "    print(\"-> b is a 16-root:\", is_nth_root(b, 16))\n",
    "    print(\"-> -1 is a square:\", is_nth_root(Fq(-1), 2))\n",
    "    print(\"-> k =\", k)\n",
    "\n",
    "    frob_X1 = X1_frobenius_polynomial(a, b, precomp_X1t)\n",
    "    print(\"-> Frobenius polynomial of X1 over Fq^k (by fast method): \", frob_X1)\n",
    "    print(\"-> Frobenius polynomial of X1 (fact.): \", frob_X1.change_ring(QQ).factor())\n",
    "    \n",
    "    Fqk = Fq.extension(k)\n",
    "    b1k = frob_X1.coefficients(sparse=False)[3]\n",
    "    b2k = frob_X1.coefficients(sparse=False)[2]\n",
    "    if Fqk(-1).is_square():\n",
    "        a1k = 2*b1k\n",
    "        a2k = 2*b2k+b1k^2\n",
    "        a3k = 2*b1k*q^k+2*b1k*b2k\n",
    "        a4k = 2*q^(2*k)+2*b1k^2*q^k+b2k^2\n",
    "    else:\n",
    "        a1k = 0\n",
    "        a2k = 2*b2k-b1k^2\n",
    "        a3k = 0\n",
    "        a4k = 2*q^(2*k)-2*b1k^2*q^k+b2k^2\n",
    "    res = L_from_Lk(C, a1k, a2k, a3k, a4k, q, k)\n",
    "\n",
    "    if len(res) > 1:\n",
    "        print(\"Warning: multiple variants passed all tests!\")\n",
    "    if len(res) == 0:\n",
    "        print(\"BUG: Coefficients are not found!\")\n",
    "    a11 = ZZ(res[0][0])\n",
    "    a21 = ZZ(res[0][1])\n",
    "    a31 = ZZ(res[0][2])\n",
    "    a41 = ZZ(res[0][3])\n",
    "    frob = T^8 + a11*T^7 + a21*T^6 + a31*T^5 + a41*T^4 + a31*q*T^3 + a21*q^2*T^2 + a11*q^3*T + q^4\n",
    "    return frob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Tests for the algorithm on small primes.\n",
    "primes = [7,11,13,17,19]\n",
    "extensions = [1,2]\n",
    "g = 4\n",
    "\n",
    "for p in primes:\n",
    "    for n in extensions:\n",
    "        q = p^n\n",
    "        Fq = FiniteField(q)\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^9 + a*x^5 + b*x\n",
    "                C = HEC(f)\n",
    "                if C == 'singular':\n",
    "                    continue\n",
    "\n",
    "                frob_C = g4_frobenius_polynomial(a,b)\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",
    "\n",
    "                frob_C_check = C.frobenius_polynomial().subs({x: T})\n",
    "                print(\"-> Frobenius polynomial of C over Fq (by Sage): \", frob_C_check)\n",
    "                if frob_C == frob_C_check:\n",
    "                    print(\"-> [ok]\")\n",
    "                else:\n",
    "                    print(\"-> [error]\")\n",
    "                print(\"-----------\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Big field tests.\n",
    "g = 4\n",
    "\n",
    "for i in range(100):\n",
    "    q = random_prime(2^24, lbound=2^23)\n",
    "    F = FiniteField(q)\n",
    "    K.<x> = F['x']\n",
    "    a = F.random_element()\n",
    "    b = F.random_element()\n",
    "    frob = g4_frobenius_polynomial(a,b)\n",
    "    print(frob)\n",
    "    print(\"Factor:\", frob.change_ring(QQ).factor())\n",
    "    N = frob(1)\n",
    "    print(\"#J(Fq) =\", N, \"(\", floor(log(N)/log(2)), \"bit) =\", factor(N))\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
}