Диплом.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Алгоритм нахождения j(tau)  [j-инварианта]\n",
    "\n",
    "# Вход: D - дискриминант числового поля\n",
    "#       a, b - коэффициенты квадратичной формы, имеющей дискриминант D.\n",
    "\n",
    "# Вывод: j-инварианта (Комплексное)\n",
    "\n",
    "def find_q(tau):\n",
    "    return exp(2*pi*sqrt(-1)*tau)\n",
    "\n",
    "def find_etha(q_tau, M):\n",
    "    s = int(0)\n",
    "    for m in range(0, M):\n",
    "        \n",
    "        s = s + CC(((-1)^m * (q_tau^(m*(3*m-1)/2))+(q_tau^(m*(3*m+1)/2))))\n",
    "    \n",
    "    s = 1 + s\n",
    "    \n",
    "    s = q_tau^(1/24) * s\n",
    "    \n",
    "    return s\n",
    "\n",
    "def find_delta(etha_tau):\n",
    "    return (2*pi)^12 * etha_tau^24\n",
    "\n",
    "def find_j_tau(D, a, b):\n",
    "    \n",
    "    d = abs(D)\n",
    "    \n",
    "    tau = (-b + sqrt(D)) / (2 * a)\n",
    "    \n",
    "    theta = -pi * b / a\n",
    "    \n",
    "    p_tau = exp(-pi * sqrt(d) / a)\n",
    "    \n",
    "    q_tau = find_q(tau)\n",
    "    \n",
    "    q_2tau = find_q(2*tau)\n",
    "    \n",
    "    precD = ceil((ln(1)+3.14*sqrt(d))/ln(10)*1) + 10\n",
    "    \n",
    "    M = ceil(sqrt(a*2/3*(precD*ln(10)+ln(6))/(pi*sqrt(d))))\n",
    "    \n",
    "    etha_tau = find_etha(q_tau, M)\n",
    "    \n",
    "    etha_2tau = find_etha(q_2tau, M)\n",
    "    \n",
    "    delta_tau = find_delta(etha_tau)\n",
    "    \n",
    "    delta_2tau = find_delta(etha_2tau)\n",
    "    \n",
    "    f = delta_2tau / delta_tau\n",
    "    \n",
    "    j_tau = (256 * f + 1)^3 / (f)\n",
    "    \n",
    "    return j_tau.n(precD)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Алгоритм построения многочлена Гильберта\n",
    "# ! Алгоритм возвращает не такой класс, как стандартный\n",
    "# В документвции алгоритм абсолютно такой же используется\n",
    "# Из-за разных классов возникают трудности с редукцией многочлена,\n",
    "# так что буду использовать стандартный алгоритм, предоставленный SAGE\n",
    "\n",
    "# Вход: D - дискриминант числового поля\n",
    "#       p - модуль\n",
    "#       a, b - коэффициенты ??\n",
    "# Вывод: H_D - многочлен Гильберта\n",
    "\n",
    "def generate_hilbert_polynomial(D, aa, bb):\n",
    "    \n",
    "    d = abs(D)\n",
    "\n",
    "    precD = ceil((ln(1) + pi * sqrt(d)) / ln(10) * 1) + 10\n",
    "    \n",
    "    R.<X> = PolynomialRing(CC, 'X')\n",
    "    \n",
    "    Priv.<X> = PolynomialRing(RR, 'X')\n",
    "    \n",
    "    H_D = 1\n",
    "    \n",
    "    H_D in R\n",
    "    \n",
    "    b = abs(D) % 2\n",
    "    \n",
    "    B = floor(sqrt(abs(D)) / 3)\n",
    "    \n",
    "    while b <= B:\n",
    "        \n",
    "        t = (b^2 - D) / 4\n",
    "        a = max(b, 1)\n",
    "        \n",
    "        while true:\n",
    "            if (t % a == 0):\n",
    "                j_tau = find_j_tau(D, aa, bb)\n",
    "                \n",
    "                if (a == b) or (b == 0):\n",
    "                    H_D = H_D * (X - j_tau)\n",
    "                else:\n",
    "                    H_D = H_D * (X^2 - 2 * real(j_tau) * X + j_tau^2)\n",
    "                    \n",
    "            a = a + 1\n",
    "                \n",
    "            if (a^2 > t):\n",
    "                break\n",
    "            \n",
    "        b = b + 2\n",
    "        \n",
    "    return H_D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Условия выбора коэффициентов кв. форм\n",
    "\n",
    "# Вход: D - дискриминант числового поля\n",
    "#       p - модуль\n",
    "#       a, b - коэффициенты ??\n",
    "# Вывод: H_D - многочлен Гильберта\n",
    "\n",
    "# def find_a_b_c(D):\n",
    "#     b^2-4ac=-D\n",
    "#     |b|<=a<=sqrt(D/3)\n",
    "#     -a<=c\n",
    "#     gcd(a,b,c)=1\n",
    "#     if |b|==a or a==c then b >= 0\n",
    "\n",
    "# Эти значения берутся из таблиц"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Алгоритм построения кривой и её скручивания по j-инварианту\n",
    "\n",
    "# Вход: j-инвариант\n",
    "# Вывод: E - кривая\n",
    "\n",
    "def build_curve(j, p):\n",
    "        \n",
    "    k = j / (1728 - j)\n",
    "    try:\n",
    "        k % p\n",
    "    except ZeroDivisionError:\n",
    "        k = 1\n",
    "    a = 3*k % p\n",
    "    b = 2*k % p\n",
    "    \n",
    "    curve = EllipticCurve(GF(p), [3*a, b])\n",
    "    \n",
    "    return curve\n",
    "\n",
    "def build_curve_twist(j, p):\n",
    "        \n",
    "    k = j / (1728 - j)\n",
    "    try:\n",
    "        k % p\n",
    "    except ZeroDivisionError:\n",
    "        k = 1\n",
    "    a = 3*k % p\n",
    "    b = 2*k % p\n",
    "    \n",
    "    N = [x for x in range(p) if kronecker(x,p) == -1]\n",
    "    \n",
    "    c = N[0]\n",
    "    \n",
    "    curve = EllipticCurve(GF(p), [3*a*c^2, b*c^3])\n",
    "    \n",
    "    return curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Метод КУ на базе многочлена Гильберта\n",
    "\n",
    "# Вывод: Эллиптическая кривая E\n",
    "\n",
    "def count_D(prev_D):\n",
    "    if prev_D == 0:\n",
    "        return -3\n",
    "    else:\n",
    "        while true:\n",
    "            prev_D = prev_D - 1\n",
    "            if (prev_D % 4 == 0):\n",
    "                return prev_D\n",
    "            if (prev_D % 4 == 1):\n",
    "                return prev_D\n",
    "\n",
    "def complex_multiplication_hilbert(t, p = 59):\n",
    "   \n",
    "    u = t - p - 1\n",
    "    \n",
    "    prev_D = 0\n",
    "    v = 0\n",
    "    D = 0\n",
    "    \n",
    "    for v in range(1,p+1):\n",
    "        D = count_D(prev_D)\n",
    "        if (4*p == u^2 + abs(D) * v^2):\n",
    "            break\n",
    "        else:\n",
    "            prev_D = D\n",
    "            \n",
    "#     print u, v, D\n",
    "    \n",
    "    H_D =  hilbert_class_polynomial(D)\n",
    "    \n",
    "#     print H_D.factor_mod(p)\n",
    "    \n",
    "    roots_of_hilbert_polynomial =[]\n",
    "    \n",
    "    for root in H_D.factor_mod(p):\n",
    "        roots_of_hilbert_polynomial.append(root[0].roots())\n",
    "    \n",
    "#     print roots_of_hilbert_polynomial\n",
    "    \n",
    "    curves = []\n",
    "    \n",
    "    for root in roots_of_hilbert_polynomial:\n",
    "        \n",
    "        try:\n",
    "            E = build_curve(root[0][0], p)\n",
    "            curves.append(E)\n",
    "        except ArithmeticError:\n",
    "            E = None\n",
    "        \n",
    "        try:\n",
    "            E_ = build_curve_twist(root[0][0], p)\n",
    "            curves.append(E_)\n",
    "        except ArithmeticError:\n",
    "            E_ = None\n",
    "    \n",
    "    return curves\n",
    "\n",
    "res = complex_multiplication_hilbert(56)\n",
    "\n",
    "for item in res:\n",
    "    print item"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_z(tau):\n",
    "    return exp(2*pi*sqrt(-1)*tau)\n",
    "\n",
    "def find_delta_for_weber(z, M):\n",
    "    s = int(0)\n",
    "    for n in range(1, M):\n",
    "        \n",
    "        s = s + CC(((-1)^n * (z^(n*(3*n-1)/2))+(z^(n*(3*n+1)/2))))\n",
    "    \n",
    "    s = 1 + s\n",
    "    \n",
    "    s = z * s^24\n",
    "    \n",
    "    return s\n",
    "\n",
    "def find_F_z(delta, z):\n",
    "    return (delta / z)^(1/24)\n",
    "\n",
    "def weber_f_0(theta,z):\n",
    "    return theta^(-1/24) * find_F_z(-theta, z) / find_F_z(theta^2,z)\n",
    "\n",
    "def weber_f_1(theta,z):\n",
    "    return theta^(-1/24) * find_F_z(theta,z) / find_F_z(theta^2,z)\n",
    "\n",
    "def weber_f_2(theta,z):\n",
    "    return sqrt(2) * theta^(1/12) * find_F_z(theta^4,z) / find_F_z(theta^2,z)\n",
    "\n",
    "\n",
    "def generate_weber_polynomial(D):\n",
    "    a = 1  #из таблиц берем\n",
    "    b = 1\n",
    "    c = 2\n",
    "    \n",
    "    d = abs(D)\n",
    "    \n",
    "    tau = (-b + sqrt(D)) / (2 * a)\n",
    "    \n",
    "    z = find_z(tau)\n",
    "    \n",
    "    precD = ceil((ln(1)+3.14*sqrt(d))/ln(10)*1) + 10\n",
    "    \n",
    "    M = ceil(sqrt(a*2/3*(precD*ln(10)+ln(6))/(pi*sqrt(d))))\n",
    "    \n",
    "    delta = find_delta_for_weber(z, M)\n",
    "    \n",
    "    theta = z^(-1)\n",
    "    \n",
    "    gamma = (weber_f_0(theta,z)^24 + 8) * (weber_f_1(theta,z)^8 - weber_f_2(theta,z)^8) / weber_f_0(theta,z)^8\n",
    "    \n",
    "    W_D = 1\n",
    "    \n",
    "    if (D % 3 == 0) and (D % 8 != 3):\n",
    "        if (D % 8 == 7):\n",
    "            W_D = hilbert_class_polynomial2(4*D, weber_f_0(theta,z)/sqrt(2))\n",
    "        if (D/4 % 8 == 2) or (D/4 % 8 == 6):\n",
    "            W_D = hilbert_class_polynomial2(D, weber_f_1(theta,z)/sqrt(2))\n",
    "        if (D/4 % 8 == 5):\n",
    "            f_ = weber_f_0(theta,z)\n",
    "            W_D = hilbert_class_polynomial2(D, f_ ^ 4)\n",
    "        if (D/4 % 8 == 1):\n",
    "            f_ = weber_f_0(theta,z)\n",
    "            W_D = hilbert_class_polynomial2(D, (f_ ^ 2) / sqrt(2))\n",
    "    if (D % 6 == 3):\n",
    "        W_D = hilbert_class_polynomial2(D, sqrt(-D) * gamma)\n",
    "    else:\n",
    "        W_D = hilbert_class_polynomial(D)\n",
    "    return W_D\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Метод КУ на базе многочлена Вебера\n",
    "\n",
    "def complex_multiplication_weber(t, p = 59):\n",
    "    D = -19\n",
    "    if (D % 4 != 0 or D % 4 != 3) and D % 8 == 3:\n",
    "        print \"D не соответствует входным условиям!\"\n",
    "        \n",
    "    W_D = generate_weber_polynomial(D)\n",
    "    \n",
    "    u = t - p - 1\n",
    "    \n",
    "    prev_D = 0\n",
    "    v = 0\n",
    "    D = 0\n",
    "    \n",
    "    for v in range(1,p+1):\n",
    "        D = count_D(prev_D)\n",
    "        if (4*p == u^2 + abs(D) * v^2):\n",
    "            break\n",
    "        else:\n",
    "            prev_D = D\n",
    "    \n",
    "    roots_of_weber_polynomial =[]\n",
    "    \n",
    "    for root in W_D.factor_mod(p):\n",
    "        roots_of_weber_polynomial.append(root[0].roots())\n",
    "        \n",
    "    true_roots_weber = []\n",
    "    \n",
    "    for root in roots_of_weber_polynomial:\n",
    "        true_roots_weber.append(root[0])\n",
    "    \n",
    "    roots_of_hilbert_polynomial = []\n",
    "    \n",
    "    for root_w in true_roots_weber:\n",
    "        roots_of_hilbert_polynomial.append(power((power(root_w[0],24) - 16), 3) / power(root_w[0],24))\n",
    "    \n",
    "    curves = []\n",
    "    \n",
    "    for root in roots_of_hilbert_polynomial:\n",
    "        \n",
    "        try:\n",
    "            E = build_curve(root, p)\n",
    "            curves.append(E)\n",
    "        except ArithmeticError:\n",
    "            E = None\n",
    "        \n",
    "        try:\n",
    "            E_ = build_curve_twist(root, p)\n",
    "            curves.append(E_)\n",
    "        except ArithmeticError:\n",
    "            E_ = None\n",
    "        \n",
    "    return curves\n",
    "\n",
    "res = complex_multiplication_weber(56)\n",
    "\n",
    "for item in res:\n",
    "    print item\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Методы КУ на базе многочлена Гильберта с условиями\n",
    "\n",
    "# Общие входные условия:\n",
    "\n",
    "# d – положительное целое число, свободное от квадратов\n",
    "# d (mod 8):      1,2,5         3              6             7\n",
    "# z_d:           sqrt(-d)  (3+sqrt(-d))/2   3+sqrt(-d)  (-3+sqrt(-d))/2   \n",
    "\n",
    "\n",
    "# p >= 5\n",
    "# d – положительное целое число, свободное от квадратов"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Алгоритм 3.1.  [d ≡ 3,2,1 (mod 4)]   d != 1 и 3\n",
    "\n",
    "def f_ot_tau(tau):\n",
    "    dzeta_48 = power(-1,(1/48))\n",
    "    return dzeta_48 * ((tau+1)/2).eta() / tau.eta()    #Возможно тут надо другую брать ЭТУ\n",
    "    \n",
    "def gamma_2(tau):\n",
    "    return (f_ot_tau(tau) ^ 24 - 16) / f_ot_tau(tau) ^ 8\n",
    "\n",
    "def gamma_3(tau):\n",
    "    power((gamma_2(tau) ^ 3 - 1728),(1/2))\n",
    "    \n",
    "def is_return_curve(d, U, V):\n",
    "    if ((d%8 == 7) and ((U-V)%4) == 1) or ((d%8 == 3 and (p%4 == (U+V) and (2*U % 2 == 0)) or ((p%4 == 1 and (2*U % 4 == 3))) or ((p%4 == 3 and (2*V % 4 == 1))))):\n",
    "        return true\n",
    "    else:\n",
    "        return false\n",
    "        \n",
    "\n",
    "def algorithm_3_1(d, z_d):\n",
    "    \n",
    "    our_j_inv = gamma_3(z_d) * sqrt(-d)\n",
    "    \n",
    "    H_D = hilbert_class_polynomial2(-d, our_j_inv)\n",
    "\n",
    "    roots_of_hilbert_polynomial = []\n",
    "    \n",
    "    for root in H_D.factor_mod(p):\n",
    "        roots_of_hilbert_polynomial.append(root[0].roots())\n",
    "    \n",
    "    betha = roots_of_hilbert_polynomial[0]\n",
    "    \n",
    "    alpha = -betha * V / U\n",
    "    \n",
    "    delta = 1728 + power(alpha,2)\n",
    "    \n",
    "    E = EllipticCurve(GF(p), 27 * power(delta, 3), 54 * alpha * power(delta, 4))\n",
    "    \n",
    "    E_ = build_curve_twist(betha, p)\n",
    "    \n",
    "    if  is_return_curve():\n",
    "        return E\n",
    "    else:\n",
    "        return E_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def algorithm_3_2(d, z_d):\n",
    "    \n",
    "    our_j_inv = gamma_3(z_d) * sqrt(-d)\n",
    "    \n",
    "    H_D = hilbert_class_polynomial2(-d, our_j_inv)\n",
    "\n",
    "    roots_of_hilbert_polynomial = []\n",
    "    \n",
    "    for root in H_D.factor_mod(p):\n",
    "        roots_of_hilbert_polynomial.append(root[0].roots())\n",
    "    \n",
    "    betha = roots_of_hilbert_polynomial[0]\n",
    "    \n",
    "    alpha = betha * V / U\n",
    "    \n",
    "    delta = 1728 + power(alpha,2)\n",
    "    \n",
    "    E = EllipticCurve(GF(p), 27 * power(delta, 3), 54 * alpha * power(delta, 4))\n",
    "    \n",
    "    E_ = build_curve_twist(betha, p)\n",
    "    \n",
    "    if (V % 4 == 1) or (V % 4 == (U-1) % 4):\n",
    "        return E\n",
    "    else:\n",
    "        return E_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def algorithm_3_3(d, z_d):\n",
    "    \n",
    "    our_j_inv = elliptic_j(z_d)\n",
    "    \n",
    "    H_D = hilbert_class_polynomial2(-d, our_j_inv)\n",
    "    \n",
    "    H_D_ = hilbert_class_polynomial(-d)\n",
    "    \n",
    "    delta = sqrt(d) % p\n",
    "    \n",
    "    H_D__ = H_D + delta * H_D_\n",
    "    \n",
    "    roots_of_hilbert_polynomial = []\n",
    "    \n",
    "    for root in H_D__.factor_mod(p):\n",
    "        roots_of_hilbert_polynomial.append(root[0].roots())\n",
    "    \n",
    "    betha = roots_of_hilbert_polynomial[0]\n",
    "    \n",
    "    alpha = betha - 1728\n",
    "    \n",
    "    etha = power_mod(alpha,((p-1)/4),p)\n",
    "    \n",
    "    E = EllipticCurve(GF(p), 27 * power(betha, 3) * alpha, 54 * alpha^2 * power(betha, 4))\n",
    "    \n",
    "    E_ = build_curve_twist(betha, p)\n",
    "\n",
    "    if (V % 2 == 0):\n",
    "        epsilon = etha % p\n",
    "        if (U % 4 == epsilon):\n",
    "            return E\n",
    "        else:\n",
    "            return E_\n",
    "    else:\n",
    "        i = delta * V / U % p\n",
    "        if ((etha == i) and (V % 4 == 3)) or ((etha != i) and (V % 4 == 1)):\n",
    "            return E\n",
    "        else:\n",
    "            return E_\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def algorithm_3_4(p):\n",
    "    \n",
    "    F = GF(p)\n",
    "    \n",
    "    A = 0\n",
    "    \n",
    "    if (U % 2 != 0) and (U-1) % 4 == V:\n",
    "        for a in F:\n",
    "            if power_mod(a, (p-1)/4, p) == -1:\n",
    "                A = a\n",
    "                break\n",
    "    \n",
    "        E_a = EllipticCurve(GF(p), A, 0)\n",
    "        return E_a\n",
    "    else:\n",
    "        if (V-1) % 4 == U:\n",
    "            for a in F:\n",
    "                if (power_mod(a, (p-1)/4, p) == (U/V)):\n",
    "                    A = a\n",
    "                    break\n",
    "    \n",
    "            E_a = EllipticCurve(GF(p), A, 0)\n",
    "            return E_a\n",
    "            \n",
    "        else:\n",
    "            V = -V\n",
    "            for a in F:\n",
    "                if power_mod(a, (p-1)/4, p) == (U/V):\n",
    "                    A = a\n",
    "                    break\n",
    "    \n",
    "            E_a = EllipticCurve(GF(p), A, 0)\n",
    "            return E_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def algorithm_3_5(p):\n",
    "    \n",
    "    F = GF(p)\n",
    "    \n",
    "    B = 0\n",
    "    \n",
    "    if (2*V % 3 == 0) and (2* U % 3 == 2):\n",
    "        E_b = EllipticCurve(GF(p), 0, 16)\n",
    "        return E_b\n",
    "    elif((2*V % 3 == 0) and (2* U % 3 == 2)):\n",
    "        if (V-1) % 4 == U:\n",
    "            for b in F:\n",
    "                if power_mod(a, (p-1)/6, p) == (p-1):\n",
    "                    B = b\n",
    "                    break\n",
    "    \n",
    "            E_b = EllipticCurve(GF(p), 0, B*16)\n",
    "            return E_b\n",
    "    elif((2*V % 3 != 0) and (2* U % 3 == 2)):\n",
    "        if (2*V) % 3 == 1:\n",
    "            for b in F:\n",
    "                if (power_mod(b, (p-1)/6, p) == (2*U/(3*V-U))):\n",
    "                    B = b\n",
    "                    break\n",
    "    \n",
    "            E_b = EllipticCurve(GF(p), 0, B*16)\n",
    "            return E_b\n",
    "            \n",
    "        else:\n",
    "            V = -V\n",
    "            for b in F:\n",
    "                if power_mod(b, ((p-1)/6), p) == (2*U/(3*V+U)):\n",
    "                    B = b\n",
    "                    break\n",
    "    \n",
    "            E_b = EllipticCurve(GF(p), 0, B*16)\n",
    "            return E_b\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 19\n",
    "\n",
    "for i in d.factor():\n",
    "    if i[1] > 1:\n",
    "        print \"ОШИБКА: d должно быть свободно от квадратов!\"\n",
    "        break\n",
    "\n",
    "d_ = d % 8\n",
    "z_d = 0\n",
    "\n",
    "if d_ == 1 or d_ == 2 or d_ == 5:\n",
    "    z_d = sqrt(-d)\n",
    "elif (d_ == 3):\n",
    "    z_d = (3 + sqrt(-d))/2\n",
    "elif (d_ == 6):\n",
    "    z_d = (3 + sqrt(-d))\n",
    "elif (d_ == 7):\n",
    "    z_d = (-3 + sqrt(-d))/2\n",
    "\n",
    "print z_d.n()\n",
    "print real(z_d)\n",
    "\n",
    "if d % 4 == 3 and d != 3:\n",
    "    E = algorithm_3_1(d, z_d)\n",
    "\n",
    "    \n",
    "if d % 4 == 2:\n",
    "    E = algorithm_3_2(d, z_d)\n",
    "    \n",
    "    \n",
    "if d % 4 == 1 and d != 1:\n",
    "    algorithm_3_3(d, z_d)\n",
    "\n",
    "    \n",
    "if d == 1:\n",
    "    algorithm_3_4(p)\n",
    "    \n",
    "\n",
    "if d == 3:\n",
    "    algorithm_3_5(p)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hilbert_class_polynomial2(D, our_j_invariant=None):\n",
    "    \n",
    "    D = Integer(D)\n",
    "    if D >= 0:\n",
    "        raise ValueError(\"D (=%s) must be negative\" % D)\n",
    "    if not (D % 4 in [0, 1]):\n",
    "        raise ValueError(\"D (=%s) must be a discriminant\" % D)\n",
    "\n",
    "    # get all primitive reduced quadratic forms, (necessary to exclude\n",
    "    # imprimitive forms when D is not a fundamental discriminant):\n",
    "\n",
    "    rqf = BinaryQF_reduced_representatives(D, primitive_only=True)\n",
    "\n",
    "    h = len(rqf) # class number\n",
    "    c1 = 3.05682737291380 # log(2*10.63)\n",
    "    c2 = sum([1/RR(qf[0]) for qf in rqf], RR(0))\n",
    "    prec =  c2*RR(3.142)*RR(D).abs().sqrt() + h*c1  # bound on log\n",
    "    prec = prec * 1.45   # bound on log_2 (1/log(2) = 1.44..)\n",
    "    prec = 10 + prec.ceil()  # allow for rounding error\n",
    "\n",
    "    # set appropriate precision for further computing\n",
    "\n",
    "    Dsqrt = D.sqrt(prec=prec)\n",
    "    R = ComplexField(prec)['t']\n",
    "    t = R.gen()\n",
    "    pol = R(1)\n",
    "    \n",
    "    if our_j_invariant != None:\n",
    "        if (our_j_invariant not in rqf):\n",
    "            pol *= (t - our_j_invariant)\n",
    "    \n",
    "#     for qf in rqf:\n",
    "#         a, b, c = list(qf)\n",
    "#         tau = (b + Dsqrt) / (a << 1)\n",
    "#         pol *= (t - elliptic_j(tau))\n",
    "\n",
    "    coeffs = [cof.real().round() for cof in pol.coefficients(sparse=False)]\n",
    "    return IntegerRing()['x'](coeffs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "our_j_invariant = 29\n",
    "print hilbert_class_polynomial(-19)\n",
    "print hilbert_class_polynomial2(-19, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = complex_multiplication_hilbert(37, 59)\n",
    "\n",
    "for curve in E:\n",
    "    print curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EE = complex_multiplication_weber(37, 59)\n",
    "\n",
    "for curve in EE:\n",
    "    print curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.6",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}