ml_based_factor.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyMC5qo/MFefi8aZvhWDv0K7",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "129fbda29c1c444ba41b0f5761e03a73": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "HBoxModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "HBoxModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "HBoxView",
            "box_style": "",
            "children": [
              "IPY_MODEL_2258c532c0fd4d5e8ce4925dc0cc2cb0",
              "IPY_MODEL_ec494eae325f4babbbecbf8e496c2ed4",
              "IPY_MODEL_56f5b7deb29d4dd9b2dced2ebd5df735"
            ],
            "layout": "IPY_MODEL_521d9d1b70924b5192a8ea785f443a10"
          }
        },
        "2258c532c0fd4d5e8ce4925dc0cc2cb0": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "HTMLModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "HTMLModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "HTMLView",
            "description": "",
            "description_tooltip": null,
            "layout": "IPY_MODEL_ab56b76c98294ce28fdd541f9e180d22",
            "placeholder": "​",
            "style": "IPY_MODEL_ebd1c2aacb044564a0c8610b2ad83ecb",
            "value": "100%"
          }
        },
        "ec494eae325f4babbbecbf8e496c2ed4": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "FloatProgressModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "FloatProgressModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "ProgressView",
            "bar_style": "success",
            "description": "",
            "description_tooltip": null,
            "layout": "IPY_MODEL_e2d1913a7cb24163953fbe82d77f8156",
            "max": 10000,
            "min": 0,
            "orientation": "horizontal",
            "style": "IPY_MODEL_233a926e29a14123b21337356176ab4e",
            "value": 10000
          }
        },
        "56f5b7deb29d4dd9b2dced2ebd5df735": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "HTMLModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "HTMLModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "HTMLView",
            "description": "",
            "description_tooltip": null,
            "layout": "IPY_MODEL_5e90d03fa8ba48198c51f8e118d057f5",
            "placeholder": "​",
            "style": "IPY_MODEL_c5f61c0ed0434d5c9cf7bd9fb012914b",
            "value": " 10000/10000 [00:35&lt;00:00, 158.29it/s]"
          }
        },
        "521d9d1b70924b5192a8ea785f443a10": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "ab56b76c98294ce28fdd541f9e180d22": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "ebd1c2aacb044564a0c8610b2ad83ecb": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "DescriptionStyleModel",
          "model_module_version": "1.5.0",
          "state": {
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "DescriptionStyleModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "StyleView",
            "description_width": ""
          }
        },
        "e2d1913a7cb24163953fbe82d77f8156": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "233a926e29a14123b21337356176ab4e": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "ProgressStyleModel",
          "model_module_version": "1.5.0",
          "state": {
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "ProgressStyleModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "StyleView",
            "bar_color": null,
            "description_width": ""
          }
        },
        "5e90d03fa8ba48198c51f8e118d057f5": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "c5f61c0ed0434d5c9cf7bd9fb012914b": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "DescriptionStyleModel",
          "model_module_version": "1.5.0",
          "state": {
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "DescriptionStyleModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "StyleView",
            "description_width": ""
          }
        }
      }
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/stblake/deep_factor/blob/main/ml_based_factor.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Experiments with Fermat's Factorisation Algorithm and Binary Classification"
      ],
      "metadata": {
        "id": "zO2GNrbIFlC9"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Sam Blake, started 1 August, 2022."
      ],
      "metadata": {
        "id": "bYZhU9WpF1PX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## A simple GMP-based implementation of Fermat's integer factorisation algorithm"
      ],
      "metadata": {
        "id": "EhZVgt5MGNV0"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Fermat's method tests various $n$ in the hopes that $a^2 - N^2$ is a perfect square, $b^2$. Then we have $N = a^2 - b^2 = (a-b)(a+b)$, if neither factor is 1, then we have a proper factor of $N$:"
      ],
      "metadata": {
        "id": "9SBANQGyGcnp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%pip install gmpy2\n",
        "import gmpy2\n",
        "from gmpy2 import mpz, mpq, mpfr\n",
        "\n",
        "def factor_fermat(n, max_iter = 65536):\n",
        "    \"\"\"Fermat's factorisation method (using GMP for fast bignum arithmetic.)\"\"\"\n",
        "\n",
        "    a = gmpy2.isqrt(n)\n",
        "    b = a**2 - n\n",
        "\n",
        "    n_iter = 0\n",
        "    while not gmpy2.is_square(b):\n",
        "        a += 1\n",
        "        b = a**2 - n\n",
        "        n_iter += 1\n",
        "        if n_iter > max_iter:\n",
        "            print(f'max_iter of {max_iter} exceeded.')\n",
        "            return mpz(1)\n",
        "\n",
        "    # print(f'n_iterations = {n_iter}')\n",
        "    return a - gmpy2.isqrt(b)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "I5VXMeHgFxFR",
        "outputId": "87fe1408-d00d-4106-d98b-9719f93f46d3"
      },
      "execution_count": 35,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Requirement already satisfied: gmpy2 in /usr/local/lib/python3.10/dist-packages (2.1.5)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "n_bits=500\n",
        "n_lsb_bits=250\n",
        "p = gmpy2.next_prime(2**n_bits + 2**n_lsb_bits)\n",
        "q = gmpy2.next_prime(2**n_bits)\n",
        "N = p*q\n",
        "p,q,N"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rEhjpFQ2GSEp",
        "outputId": "9c8fecfc-76dd-4dc1-d20c-82ec62134319"
      },
      "execution_count": 36,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(mpz(3273390607896141870013189696827599152216642046043064789483291368096133796406483806277603157879397453791647432687767414519617890359001917803452170240213),\n",
              " mpz(3273390607896141870013189696827599152216642046043064789483291368096133796404674554883270092325904157150886684127560071009217256545885393053328527589431),\n",
              " mpz(10715086071862673209484250490600018105614048117055336074437503883703510511255283611453516643897120398782593848908414486810818812365967986162097137667305491643325337101150754402778383578281888490209299974003439160965942719766577633095563804179706669423941805490776263279775622639617075379411613109988803))"
            ]
          },
          "metadata": {},
          "execution_count": 36
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "factor_fermat(N)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BIW0mLGdN5J1",
        "outputId": "108f101c-b46c-4cc0-b31e-05e7633349fe"
      },
      "execution_count": 37,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "mpz(3273390607896141870013189696827599152216642046043064789483291368096133796404674554883270092325904157150886684127560071009217256545885393053328527589431)"
            ]
          },
          "metadata": {},
          "execution_count": 37
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "When factoring $N = p\\,q$, Fermat's factorisation method works best when $p$ and $q$ are close to $\\sqrt{N}$. Otherwise, if we have an approximation $u/v$ to $p/q$, then we can use Fermat's method to factor $u \\, v\\, N$, providing $u/v$ is _sufficiently close_ to $p/q$. Lehman (\"Factoring Large Integers\", _Mathematics of Computation_, vol. 28, no. 126, pp. 637-646, 1974) credits this extension of Fermat's algorithm to Lawrence (\"Factorisation of numbers\", _Messenger of Mathematics_, vol. 24, pp. 100-109, 1895)."
      ],
      "metadata": {
        "id": "ZV34T7BQOmcA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def factor_lawrence(n, u_v_ratio, max_iter = 100000):\n",
        "    \"\"\"Lawrence's extension to Fermat's factoring algorithm.\"\"\"\n",
        "    u,v = u_v_ratio.numerator, u_v_ratio.denominator\n",
        "    return gmpy2.gcd(factor_fermat(u*v*n, max_iter),n)"
      ],
      "metadata": {
        "id": "IXFOWMH7OMwM"
      },
      "execution_count": 38,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "factor_lawrence(748543215795445052722625573101291605706283989, mpq(210381,144089))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ixXRCB91OxdI",
        "outputId": "4629f444-caf7-447d-be70-ab2abe5aeee8"
      },
      "execution_count": 39,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "mpz(33059500175075655435169)"
            ]
          },
          "metadata": {},
          "execution_count": 39
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Generate ML training data - semiprimes"
      ],
      "metadata": {
        "id": "TwjOcywJlBKy"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import random\n",
        "\n",
        "def random_prime_pair(n_bits, ratio_interval):\n",
        "    \"\"\"random_prime_pair returns a semiprime and a pair of random primes with ratio \\\n",
        "    randomly within the ratio_interval, which is given by [mpq(n1,d1), mpq(n2,d2)]. \\\n",
        "    The product of the two semiprimes will be n_bits bits.\"\"\"\n",
        "\n",
        "    interval_lower, interval_upper = ratio_interval\n",
        "    assert interval_lower < interval_upper\n",
        "\n",
        "    while True:\n",
        "      rn = random.randint(0,2**64 - 1)\n",
        "      rd = random.randint(1,2**64 - 1)\n",
        "\n",
        "      if rn > rd:\n",
        "          rn,rd = rd,rn\n",
        "      rq = mpq(rn, rd) # 0 <= rq <= 1\n",
        "\n",
        "      r = interval_lower + (interval_upper - interval_lower)*rq\n",
        "\n",
        "      r_nbits = int(gmpy2.floor(gmpy2.log2(r)) + 1)\n",
        "\n",
        "      rmin, rmax = 2**(n_bits//2 - r_nbits//2 - 1), 2**(n_bits//2 - r_nbits//2)\n",
        "      rint_q = mpz(random.randint(rmin,rmax))\n",
        "      rint_p, _ = gmpy2.t_divmod(r.numerator*rint_q, r.denominator)\n",
        "\n",
        "      p = gmpy2.next_prime(rint_p)\n",
        "      q = gmpy2.next_prime(rint_q)\n",
        "      N = gmpy2.mul(p,q)\n",
        "      if len(N.digits(2)) == n_bits:\n",
        "        break\n",
        "\n",
        "    # We want p/q > 1\n",
        "    if p < q:\n",
        "      p,q = q,p\n",
        "\n",
        "    return N,p,q"
      ],
      "metadata": {
        "id": "cq3p9AA7PA7d"
      },
      "execution_count": 40,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "N,p,q = random_prime_pair(426, [2, 3])\n",
        "N, p, q, p/q, len(p.digits(2)), len(q.digits(2)), len(N.digits(2))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "o5U1frzwlQjC",
        "outputId": "8e178103-3424-42c1-c8fd-f942fc33a769"
      },
      "execution_count": 41,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(mpz(120089546046838015577833287908287976100839426355341480832826875075414395903389250892591544925308138505381476557658929939827648861),\n",
              " mpz(18714422895919543347365009002775639751333304579358778505078802377),\n",
              " mpz(6416951605439145464960302617777176462811630045334299245594101493),\n",
              " mpfr('2.9164039323682602'),\n",
              " 214,\n",
              " 212,\n",
              " 426)"
            ]
          },
          "metadata": {},
          "execution_count": 41
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Generate [semi-prime, p/q-ratio decision] pairs for training:"
      ],
      "metadata": {
        "id": "vs7LuRWiofwj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import tqdm.notebook as tq\n",
        "\n",
        "def generate_training_semiprimes(n_semiprime_bits, \\\n",
        "    min_ratio, max_ratio, ratio_diff_scale, n_training_samples):\n",
        "    \"\"\"generate_training_semiprimes returns a list of /n_training_samples/ pairs of the form [/semiprime/, p/q-decision], used for \\\n",
        "    subsequent training for the binary classification problem.\"\"\"\n",
        "    diff = max_ratio - min_ratio\n",
        "    diff *= ratio_diff_scale\n",
        "    min_extended_ratio_mpq = max(1, gmpy2.f2q(min_ratio - diff))\n",
        "    max_extended_ratio_mpq = max(1, gmpy2.f2q(max_ratio + diff))\n",
        "    pq_ratio_interval = [min_extended_ratio_mpq, max_extended_ratio_mpq]\n",
        "\n",
        "    assert min_extended_ratio_mpq < max_extended_ratio_mpq\n",
        "\n",
        "    n_inside = 0\n",
        "    n_outside = 0\n",
        "\n",
        "    semiprimes_inside = dict()\n",
        "    semiprimes_outside = dict()\n",
        "\n",
        "    pbar = tq.tqdm(total = n_training_samples)\n",
        "\n",
        "    while n_inside + n_outside < n_training_samples:\n",
        "\n",
        "        # generate random prime pair within a specified p/q-ratio interval.\n",
        "        N,p,q = random_prime_pair(n_semiprime_bits, pq_ratio_interval)\n",
        "        ratio = gmpy2.div(p,q)\n",
        "        N = str(N)\n",
        "\n",
        "        if n_inside <= n_training_samples//2 and min_ratio < ratio < max_ratio:\n",
        "            pbar.update(1)\n",
        "            n_inside += 1\n",
        "            semiprimes_inside[N] = True\n",
        "\n",
        "        if n_outside <= n_training_samples//2 and not min_ratio < ratio < max_ratio \\\n",
        "        and min_ratio - diff < ratio < max_ratio + diff:\n",
        "            pbar.update(1)\n",
        "            n_outside += 1\n",
        "            semiprimes_outside[N] = True\n",
        "\n",
        "    pbar.close()\n",
        "\n",
        "    training_data = []\n",
        "\n",
        "    for semiprime in semiprimes_inside.keys():\n",
        "        training_data.append([semiprime, 1])\n",
        "\n",
        "    for semiprime in semiprimes_outside.keys():\n",
        "        training_data.append([semiprime, 0])\n",
        "\n",
        "    random.shuffle(training_data)\n",
        "\n",
        "    return training_data"
      ],
      "metadata": {
        "id": "rxskHAitlcvP"
      },
      "execution_count": 42,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 1 million 426-bit semiprimes, with p/q ratio in [2,3]\n",
        "\n",
        "from os.path import exists\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "\n",
        "min_ratio = 2\n",
        "max_ratio = 3\n",
        "ratio_diff_scale = mpq(1,2)\n",
        "n_semiprime_bits = 426\n",
        "n_train = 10**4\n",
        "\n",
        "training_data = generate_training_semiprimes(\\\n",
        "  n_semiprime_bits = n_semiprime_bits, \\\n",
        "  min_ratio = min_ratio, \\\n",
        "  max_ratio = max_ratio, \\\n",
        "  ratio_diff_scale = ratio_diff_scale, \\\n",
        "  n_training_samples = n_train)\n",
        "\n",
        "# Write training data to file.\n",
        "\n",
        "file_name = f'training_data_{n_semiprime_bits}_ratio_{min_ratio}_{max_ratio}.h5'\n",
        "df_training_data = pd.DataFrame(training_data)\n",
        "\n",
        "if exists(file_name):\n",
        "    df_training_data.to_hdf(\n",
        "    file_name,\n",
        "    key='semiprimes',\n",
        "    append = True,\n",
        "    mode = 'r+',\n",
        "    format = 'table')\n",
        "else:\n",
        "  df_training_data.to_hdf(\n",
        "    file_name,\n",
        "    min_itemsize = n_semiprime_bits,\n",
        "    key = 'semiprimes',\n",
        "    format = 'table')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 49,
          "referenced_widgets": [
            "129fbda29c1c444ba41b0f5761e03a73",
            "2258c532c0fd4d5e8ce4925dc0cc2cb0",
            "ec494eae325f4babbbecbf8e496c2ed4",
            "56f5b7deb29d4dd9b2dced2ebd5df735",
            "521d9d1b70924b5192a8ea785f443a10",
            "ab56b76c98294ce28fdd541f9e180d22",
            "ebd1c2aacb044564a0c8610b2ad83ecb",
            "e2d1913a7cb24163953fbe82d77f8156",
            "233a926e29a14123b21337356176ab4e",
            "5e90d03fa8ba48198c51f8e118d057f5",
            "c5f61c0ed0434d5c9cf7bd9fb012914b"
          ]
        },
        "id": "63J3tEIjrELN",
        "outputId": "38c72674-b1ee-4c2d-cc1d-5198cb842c9a"
      },
      "execution_count": 57,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "  0%|          | 0/10000 [00:00<?, ?it/s]"
            ],
            "application/vnd.jupyter.widget-view+json": {
              "version_major": 2,
              "version_minor": 0,
              "model_id": "129fbda29c1c444ba41b0f5761e03a73"
            }
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "file_name = f'training_data_{n_semiprime_bits}_ratio_{min_ratio}_{max_ratio}.h5'\n",
        "\n",
        "training_data = pd.read_hdf(file_name)\n",
        "training_data[0] = training_data[0].apply(mpz)\n",
        "training_data = training_data.values\n",
        "n_train = training_data.shape[0]\n",
        "n_train"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TsH7uAcgy9xw",
        "outputId": "5b77b0fb-dd6e-46a3-b897-045d2c5bde31"
      },
      "execution_count": 59,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "20000"
            ]
          },
          "metadata": {},
          "execution_count": 59
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Rational base representations:"
      ],
      "metadata": {
        "id": "gesrMblztUkb"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def rat_base(n, b):\n",
        "    \"\"\"rat_base computes the (rational) base b representation of the positive\n",
        "    integer n, where b is a gmpy2 mpq object (a rational number).\"\"\"\n",
        "\n",
        "    if type(b) is int or type(b) is mpz:\n",
        "        return [int(k,b) for k in mpz(n).digits(b)]\n",
        "    elif type(b) is float:\n",
        "        b = gmpy2.f2q(b)\n",
        "\n",
        "    if type(b) is not mpq:\n",
        "        print('ERROR: type(b) == mpq.')\n",
        "\n",
        "    if b < 1:\n",
        "        print('ERROR: b > 1.')\n",
        "        return n\n",
        "\n",
        "    if n < 0:\n",
        "        print('ERROR: n > 0')\n",
        "        return n\n",
        "\n",
        "    if n == 0:\n",
        "        return [0]\n",
        "\n",
        "    m = n\n",
        "    base_rep = []\n",
        "    while m > 0:\n",
        "        d = gmpy2.f_mod(m, b.numerator)\n",
        "        m = gmpy2.f_div(m, b.numerator)*b.denominator\n",
        "        base_rep.append(int(d))\n",
        "\n",
        "    return base_rep[::-1]"
      ],
      "metadata": {
        "id": "5se57zP8sQSy"
      },
      "execution_count": 44,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "rat_base(121201, mpq(3,2))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Q6FTjajutrIR",
        "outputId": "a06f2dc4-d484-4c95-c8d3-2d971d52a639"
      },
      "execution_count": 45,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[2, 1, 2, 2, 1, 1, 1, 2, 0, 2, 1, 1, 2, 2, 2, 0, 2, 0, 0, 2, 0, 0, 0, 1, 1, 1]"
            ]
          },
          "metadata": {},
          "execution_count": 45
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Reshaping multi-feature training data:"
      ],
      "metadata": {
        "id": "TDzjeTgXuHrp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def pad_left(lst, n):\n",
        "    return [0]*(n - len(lst)) + lst"
      ],
      "metadata": {
        "id": "6TNrBHJxtvxs"
      },
      "execution_count": 46,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def reshape_training_data(training_data, base):\n",
        "\n",
        "    n_samples = len(training_data)\n",
        "    max_len = len(rat_base(max(training_data[:,0]), base))\n",
        "\n",
        "    X = np.zeros((n_samples, max_len), dtype = np.float32)\n",
        "    for k, (train_d,_) in tq.tqdm(enumerate(training_data), total = n_samples):\n",
        "        X[k,:] = pad_left(rat_base(train_d, base), max_len)\n",
        "\n",
        "    Y = np.array([classification_d for _,classification_d in training_data])\n",
        "\n",
        "    return X, Y"
      ],
      "metadata": {
        "id": "U8vd5cNGt71G"
      },
      "execution_count": 47,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Train model"
      ],
      "metadata": {
        "id": "plVMALMZuj9n"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def create_baby_model():\n",
        "    # Create model\n",
        "    model = Sequential()\n",
        "\n",
        "    model.add(Dense(n_semiprime_bits//4, \\\n",
        "        kernel_regularizer=regularizers.l2(0.0005), \\\n",
        "        input_shape=(n_semiprime_bits,), \\\n",
        "        activation='relu'))\n",
        "\n",
        "    model.add(Dropout(0.2))\n",
        "    model.add(Dense(n_semiprime_bits//4, \\\n",
        "        kernel_regularizer=regularizers.l2(0.0005), \\\n",
        "        activation='relu'))\n",
        "\n",
        "    model.add(Dropout(0.2))\n",
        "    model.add(Dense(1, activation='sigmoid'))\n",
        "\n",
        "    # Compile model\n",
        "    model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])\n",
        "    return model"
      ],
      "metadata": {
        "id": "oe8lcExRwor_"
      },
      "execution_count": 55,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Convert training data to binary.\n",
        "from sklearn.model_selection import train_test_split\n",
        "\n",
        "base = 2\n",
        "\n",
        "X, Y = reshape_training_data(training_data, base)\n",
        "X_train, X_test, y_train, y_test = train_test_split(\\\n",
        "  X, Y, test_size = 0.33, shuffle = False)"
      ],
      "metadata": {
        "id": "FZan4_6at-UG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Fit model to training data.\n",
        "\n",
        "from tensorflow.keras.models import Sequential\n",
        "from tensorflow.keras.layers import Dense, Dropout, Embedding\n",
        "from tensorflow.keras import regularizers\n",
        "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
        "\n",
        "%pip install scikeras\n",
        "from scikeras.wrappers import KerasClassifier\n",
        "\n",
        "es = EarlyStopping(monitor='accuracy', mode = 'max', restore_best_weights = True, \\\n",
        "  min_delta = 0.001, patience = 25, verbose = 1)\n",
        "\n",
        "mc = ModelCheckpoint(\\\n",
        "  f'baby_model_n_semiprime_bits{n_semiprime_bits}_ratio_{min_ratio}_{max_ratio}_base_{base}.h5', \\\n",
        "  monitor = 'accuracy', mode = 'min', save_best_only = True)\n",
        "\n",
        "estimator = KerasClassifier(model=create_baby_model, epochs=1000, \\\n",
        "  batch_size=100_000, verbose=1, callbacks = [es, mc])\n",
        "\n",
        "fitted_model = estimator.fit(X_train, y_train)"
      ],
      "metadata": {
        "id": "zCBniuSbuiKu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Predict on test data.\n",
        "from sklearn.metrics import accuracy_score\n",
        "\n",
        "yhat = fitted_model.predict(X_test)\n",
        "acc = accuracy_score(y_test, yhat)\n",
        "print(f'accuracy = {acc:.4f}')"
      ],
      "metadata": {
        "id": "nkagB6x6vFnW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.metrics import confusion_matrix\n",
        "cm = confusion_matrix(y_test,yhat)\n",
        "cm/np.sum(cm)"
      ],
      "metadata": {
        "id": "xywFrzHHvFqs"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "bases = [2,3,4,5,6,7,8,9,10]\n",
        "classifications_is = []\n",
        "classifications_oos = []\n",
        "accuracies = []\n",
        "\n",
        "for base in bases:\n",
        "    n_semiprime_bits_base_N = len(rat_base(max(training_data[:,0]), base))\n",
        "\n",
        "    # Convert training data to base.\n",
        "    X, Y = reshape_training_data(training_data, base)\n",
        "    X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = 0.33, shuffle = False)\n",
        "\n",
        "    # Fit model to training data.\n",
        "    es = EarlyStopping(monitor='accuracy', mode = 'max', restore_best_weights = True,\n",
        "      min_delta = 0.001, patience = 25, verbose = 1)\n",
        "    mc = ModelCheckpoint(\n",
        "      f'baby_model_n_semiprime_bits{n_semiprime_bits}_ratio_{min_ratio}_{max_ratio}_base_{base}.h5',\n",
        "      monitor = 'accuracy', mode = 'min', save_best_only = True)\n",
        "    estimator = KerasClassifier(model=create_baby_model, epochs=1000, batch_size=100_000,\n",
        "      verbose=1, callbacks = [es, mc])\n",
        "    fitted_model = estimator.fit(X_train, y_train)\n",
        "\n",
        "    # Predict on test data.\n",
        "    yhat = fitted_model.predict(X_test)\n",
        "    acc = accuracy_score(y_test, yhat)\n",
        "    print(f'accuracy = {acc:.4f}')\n",
        "\n",
        "    # Confusion matrix.\n",
        "    cm = confusion_matrix(y_test,yhat)\n",
        "    print(cm/np.sum(cm))\n",
        "\n",
        "    accuracies.append(acc)\n",
        "    classifications_is.append(fitted_model.predict(X_train))\n",
        "    classifications_oos.append(yhat)"
      ],
      "metadata": {
        "id": "qLWIu5gvvjkJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "[[base,acc] for base,acc in zip(bases,accuracies)]"
      ],
      "metadata": {
        "id": "9V5_e6qDvmVp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Pearson correlation matrix:"
      ],
      "metadata": {
        "id": "KkhGP1vCvzdE"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "correlations_df = classifications_df.drop('target',axis=1).corr(method='pearson')\n",
        "correlations_df.style.background_gradient(cmap='Blues').format(precision=2)"
      ],
      "metadata": {
        "id": "f2fAslZVvmYW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Ensemble classification:"
      ],
      "metadata": {
        "id": "HsTNGF0ovv4C"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%pip install xgboost\n",
        "from xgboost import XGBClassifier\n",
        "\n",
        "X, Y = classifications_df.drop('target', axis=1).values, classifications_df['target'].values\n",
        "\n",
        "model_xgb = XGBClassifier()\n",
        "kfold = StratifiedKFold(n_splits=10, shuffle=True)\n",
        "scores = cross_val_score(model_xgb, X, Y, scoring='accuracy', cv=kfold, n_jobs=-1)\n",
        "\n",
        "print(f'model accuracy: {scores.mean():.2f} +/- {scores.std():.2f}')"
      ],
      "metadata": {
        "id": "cvN53aXlvufj"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}