{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9e2b992c-a7de-4f66-89a4-8443ba963895",
   "metadata": {},
   "source": [
    "This notebook illustrates the application of the integer relation detection algorithm `PSLQ`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd8252d9-f3b5-4a61-a58a-dc40472df941",
   "metadata": {},
   "source": [
    "# 1. The Integer Relation Detection Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e3d64e7-71f5-4701-b6ee-25c44f25af9f",
   "metadata": {},
   "source": [
    "Let $x$ be a vector of $n$ real numbers,\n",
    "$x = (x_1,x_2,\\ldots,x_n)$.\n",
    "\n",
    "An integer relation algorithm finds $n$ integers \n",
    "$a = (a_1,a_2,\\ldots,a_n)$:\n",
    "\n",
    "$$\n",
    "   a_1 x_1 + a_2 x_2 + \\cdots + a_n x_n = 0.\n",
    "$$\n",
    "\n",
    "The Euclidean algorithm solves this problem for $n=2$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30fc0252-9490-411f-b450-1abb5af0520d",
   "metadata": {},
   "source": [
    "The PSLQ algorithm by David Bailey and Helaman Ferguson\n",
    "is the best-known integer relation algorithm.\n",
    "PSLQ was named one of ten ``algorithms of the century'' \n",
    "by Computing in Science and Engineering.\n",
    "\n",
    "PSLQ stands for ``Partial Sum of Least Squares.''\n",
    "\n",
    "The algorithm requires high-precision arithmetic software:\n",
    "at least $d \\times n$ digits, where $d$ is the size (in digits)\n",
    "of the largest integer $a_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "addafd28-fd37-4477-a65a-8b203a590dc1",
   "metadata": {},
   "source": [
    "One application of PSLQ is to find polynomials with\n",
    "integer coefficients that vanish at approximate roots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3193ac0a-e502-43f7-90ae-ac31d15184ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpmath import pslq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "35ce562e-9cae-46bb-abe0-4fa8b9f41c02",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5874010519681994"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = 4.0**(1/3)\n",
    "r  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "993624d9-08b4-4ee5-a356-21d71b1a90b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4, 0, 0, -1, 0]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pslq([r**n for n in range(5)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bc0888e-89cc-4c15-b75d-7f3bceff603a",
   "metadata": {},
   "source": [
    "The calculations above simulated the discovery of\n",
    "1.5874010519681994 as $\\sqrt[3]{4}$, as a root of $4 - x^3 = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "163601bd-18d4-4f02-a833-e6f548f3576a",
   "metadata": {},
   "source": [
    "# 2. a hexadecimal expansion for $\\pi$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b5d36be-ee54-422d-a9d2-54ca3f81bb65",
   "metadata": {},
   "source": [
    "One of the remarkable discoveries made by PSLQ\n",
    "is a simple formula that allows to calculating any binary digit\n",
    "of $\\pi$ without calculating the digits preceding it:\n",
    "\n",
    "$$\n",
    "  \\pi = \\sum_{i=0}^\\infty \\frac{1}{16^i}\n",
    "    \\left(\n",
    "      \\frac{4}{8i+1} - \\frac{2}{8i+4} - \\frac{1}{8i+5} - \\frac{1}{8i+6}\n",
    "    \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bed929c-e7f3-4853-b0aa-b48b6555fd4b",
   "metadata": {},
   "source": [
    "* David H. Bailey, Peter B. Borwein, and Simon Plouffe.\n",
    "  **On the Rapid Computation of Various Polylogarithmic Constants.**\n",
    "  *Mathematics of Computation* 66 (218): 903--913, 1997.\n",
    "\n",
    "* David H. Bailey.\n",
    "  **Integer Relation Detection.**\n",
    "  *Computing in Science and Engineering* 2(1):24--28, 2000."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dafa91c1",
   "metadata": {},
   "source": [
    "We will need the package ``mpmath``, which comes with the installation of ``sympy``.\n",
    "Let us verify if we can import ``pslq`` of the ``mpmath`` package ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3b451a11",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpmath import pslq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6c5fd95",
   "metadata": {},
   "source": [
    "The code in the cell below computes the coefficients in the spigot formula for pi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "633d866d-2a6c-45f0-9ba8-60211fcab2ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0d427fb6-1858-484d-915c-57f8f89998e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.007184476410883,\n",
       " 0.5064768766636942,\n",
       " 0.33923024524151824,\n",
       " 0.2554128118793683,\n",
       " 0.2050025576328489,\n",
       " 0.17131707066297353,\n",
       " 0.1472019346691604]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = [sum([1.0/(16**k*(8*k+j)) for k in range(8)]) for j in range(1,8)]\n",
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2a9b8265-e2f2-4182-85e8-5ea8d7da97e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-4, 0, 0, 2, 1, 1, 0, 1]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S.append(math.pi)\n",
    "R = pslq(S)\n",
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33e112d0",
   "metadata": {},
   "source": [
    "We recognize the coefficients of the expansion in the above formula."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}