{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diophantine equation (Euler Problem 66)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Consider quadratic Diophantine equations of the form:\n",
    "\n",
    "$x^2 – Dy^2 = 1$\n",
    "\n",
    "For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.\n",
    "\n",
    "It can be assumed that there are no solutions in positive integers when D is square.\n",
    "\n",
    "By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:\n",
    "\n",
    "$3^2 – 2×2^2 = 1$\n",
    "\n",
    "$2^2 – 3×1^2 = 1$\n",
    "\n",
    "$9^2 – 5×4^2 = 1$\n",
    "\n",
    "$5^2 – 6×2^2 = 1$\n",
    "\n",
    "$8^2 – 7×3^2 = 1$\n",
    "\n",
    "Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.\n",
    "\n",
    "Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def get_floor_sqrt(n):\n",
    "    return math.floor(math.sqrt(n))\n",
    "\n",
    "assert(get_floor_sqrt(5) == 2)\n",
    "assert(get_floor_sqrt(23) == 4)\n",
    "\n",
    "def next_expansion_1(current_a, current_nominator, current_denominator, original_number):\n",
    "    cn = current_nominator\n",
    "    cd = current_denominator\n",
    "    cn, cd = cd, (original_number - cd * cd) // cn\n",
    "    next_a = math.floor((math.sqrt(original_number) + cn) // cd)\n",
    "    cn = cn - next_a * cd\n",
    "    cn *= -1\n",
    "    \n",
    "    return next_a, cn, cd\n",
    "\n",
    "def get_continued_fraction_sequence(n):\n",
    "    \n",
    "    # If number is a square number we return it.\n",
    "    floor_sqrt = get_floor_sqrt(n)\n",
    "    if n == floor_sqrt * floor_sqrt:\n",
    "        return ((floor_sqrt), [])\n",
    "    \n",
    "    # Otherwise, we calculate the next expansion till we\n",
    "    # encounter a step a second time.    \n",
    "    a = floor_sqrt\n",
    "    cn = a\n",
    "    cd = 1\n",
    "    sequence = []\n",
    "    previous_steps = []\n",
    "    \n",
    "    while not (a, cn, cd) in previous_steps :\n",
    "        #print(\"a: {} cn: {} cd: {}\".format(a, cn, cd))\n",
    "        previous_steps.append((a, cn, cd))\n",
    "        a, cn, cd = next_expansion_1(a, cd, cn, n)\n",
    "        sequence.append(a)\n",
    "    sequence.pop()\n",
    "    return ((floor_sqrt), sequence)\n",
    "\n",
    "def get_continued_fraction_sequence(n):\n",
    "    \n",
    "    # If number is a square number we return it.\n",
    "    floor_sqrt = get_floor_sqrt(n)\n",
    "    if n == floor_sqrt * floor_sqrt:\n",
    "        return ((floor_sqrt), [])\n",
    "    \n",
    "    # Otherwise, we calculate the next expansion till we\n",
    "    # encounter a step a second time.    \n",
    "    a = floor_sqrt\n",
    "    cn = a\n",
    "    cd = 1\n",
    "    sequence = []\n",
    "    previous_steps = []\n",
    "    \n",
    "    while not (a, cn, cd) in previous_steps :\n",
    "        #print(\"a: {} cn: {} cd: {}\".format(a, cn, cd))\n",
    "        previous_steps.append((a, cn, cd))\n",
    "        a, cn, cd = next_expansion_1(a, cd, cn, n)\n",
    "        sequence.append(a)\n",
    "    sequence.pop()\n",
    "    return ((floor_sqrt), sequence)\n",
    "\n",
    "assert(get_continued_fraction_sequence(1) == ((1), []))\n",
    "assert(get_continued_fraction_sequence(4) == ((2), []))\n",
    "assert(get_continued_fraction_sequence(25) == ((5), []))\n",
    "assert(get_continued_fraction_sequence(2) == ((1), [2]))\n",
    "assert(get_continued_fraction_sequence(3) == ((1), [1,2]))\n",
    "assert(get_continued_fraction_sequence(5) == ((2), [4]))\n",
    "assert(get_continued_fraction_sequence(13) == ((3), [1,1,1,1,6]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gcd(a, b):\n",
    "    if b > a:\n",
    "        a, b = b, a\n",
    "    while a % b != 0:\n",
    "        a, b = b, a % b\n",
    "    return b\n",
    "\n",
    "def add_fractions(n1, d1, n2, d2):\n",
    "    d = d1 * d2\n",
    "    n1 = n1 * (d // d1)\n",
    "    n2 = n2 * (d // d2)\n",
    "    n = n1 + n2\n",
    "    p = gcd(n, d)\n",
    "    return (n // p, d // p)\n",
    "\n",
    "def next_expansion_2(previous_numerator, previous_denumerator, value):\n",
    "    if previous_numerator == 0:\n",
    "        return (value, 1)\n",
    "    return add_fractions(previous_denumerator, previous_numerator, value, 1)\n",
    "    \n",
    "def get_fractions(n, x):\n",
    "    #print(\"get_fractions(n={}, x={})\".format(n, x))\n",
    "    # Get sequence of a_x\n",
    "    first_value, sequence = get_continued_fraction_sequence(n)\n",
    "    sequence = [first_value] + math.ceil(x / len(sequence)) * sequence\n",
    "    sequence = sequence[:x + 1]\n",
    "    sequence = sequence[::-1]\n",
    "    \n",
    "    n, d = 0, 1\n",
    "    for s in sequence:\n",
    "        n, d = next_expansion_2(n, d, s)\n",
    "    return (n, d)\n",
    "\n",
    "assert(get_fractions(7, 0) == (2, 1))\n",
    "assert(get_fractions(7, 1) == (3, 1))\n",
    "assert(get_fractions(7, 2) == (5, 2))\n",
    "assert(get_fractions(7, 3) == (8, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_minimal_solution(d):\n",
    "    for i in range(0, 100):\n",
    "        x, y = get_fractions(d, i)\n",
    "        if x * x - d * y * y == 1:\n",
    "            return((x, y))\n",
    "\n",
    "assert(get_minimal_solution(2) == (3, 2))\n",
    "assert(get_minimal_solution(3) == (2, 1))\n",
    "assert(get_minimal_solution(85) == (285769, 30996))\n",
    "assert(get_minimal_solution(109) == (158070671986249, 15140424455100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt\n",
    "\n",
    "def is_square(n):\n",
    "    return sqrt(n).is_integer()\n",
    "\n",
    "print(is_square(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d: 661 x: 16421658242965910275055840472270471049\n"
     ]
    }
   ],
   "source": [
    "x_max = 0\n",
    "d_max = 0\n",
    "\n",
    "for d in range(2, 1001):\n",
    "    if is_square(d):\n",
    "        continue\n",
    "    \n",
    "    x, y = get_minimal_solution(d)\n",
    "    if x > x_max:\n",
    "        x_max = x\n",
    "        d_max = d\n",
    "\n",
    "print(\"d: {} x: {}\".format(d_max, x_max))\n",
    "s = d_max\n",
    "assert(s == 661)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "completion_date": "Mon, 28 Jan 2019, 21:08",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python3.6",
   "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.6.5"
  },
  "tags": [
   "diophantine",
   "equation",
   "pell's equation"
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 2
}