{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spiral primes (Euler Problem 58)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "[https://projecteuler.net/problem=58](https://projecteuler.net/problem=58)\n",
    "\n",
    "Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.\n",
    "\n",
    "    37 36 35 34 33 32 31\n",
    "    38 17 16 15 14 13 30\n",
    "    39 18  5  4  3 12 29\n",
    "    40 19  6  1  2 11 28\n",
    "    41 20  7  8  9 10 27\n",
    "    42 21 22 23 24 25 26\n",
    "    43 44 45 46 47 48 49\n",
    "\n",
    "It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.\n",
    "\n",
    "If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_last_corner_value(side_length):\n",
    "    return side_length * side_length\n",
    "\n",
    "assert(get_last_corner_value(1) == 1)\n",
    "assert(get_last_corner_value(3) == 9)\n",
    "assert(get_last_corner_value(5) == 25)\n",
    "\n",
    "\n",
    "def get_corner_values(side_length):\n",
    "    if side_length == 1:\n",
    "        return [1]\n",
    "    return [get_last_corner_value(side_length) - i * (side_length - 1)\n",
    "            for i in range(0, 4)][::-1]\n",
    "\n",
    "assert(get_corner_values(1) == [1])\n",
    "assert(get_corner_values(3) == [3, 5, 7, 9])\n",
    "assert(get_corner_values(5) == [13, 17, 21, 25])\n",
    "assert(get_corner_values(7) == [31, 37, 43, 49])\n",
    "\n",
    "\n",
    "def get_diagonal_values(side_length):\n",
    "    return [corner_value\n",
    "        for length in range(1, side_length + 1, 2)\n",
    "        for corner_value in get_corner_values(length)\n",
    "    ]\n",
    "\n",
    "assert(get_diagonal_values(1) == [1])\n",
    "assert(get_diagonal_values(3) == [1, 3, 5, 7, 9])\n",
    "assert(get_diagonal_values(5) == [1, 3, 5, 7, 9, 13, 17, 21, 25])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expmod(base, exp, m):\n",
    "    if exp == 0:\n",
    "        return 1\n",
    "    if (exp % 2 == 0):\n",
    "        return (expmod(base, exp // 2, m) ** 2 % m)\n",
    "    return (base * expmod(base, exp - 1, m) % m)\n",
    "    \n",
    "\n",
    "def fermat_test_sicp(n):\n",
    "    if n < 40:\n",
    "        return n in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 39]\n",
    "    a = n - 3\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 5\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 7\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 11\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 13\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 17\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 19\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 23\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 29\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 31\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 37\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    a = n - 39\n",
    "    if not expmod(a, n, n) == a:\n",
    "        return False\n",
    "    return True\n",
    "\n",
    "def trial_division(n):\n",
    "    a = []  \n",
    "    if n % 2 == 0:\n",
    "        a.append(2)\n",
    "    while n % 2 == 0:\n",
    "        n //= 2\n",
    "    f = 3\n",
    "    while f * f <= n:\n",
    "        if n % f == 0:\n",
    "            a.append(f)\n",
    "            while n % f == 0:\n",
    "                n //= f\n",
    "        else:\n",
    "            f += 2   \n",
    "    if n != 1:\n",
    "        a.append(n)\n",
    "    return a\n",
    "\n",
    "def is_prime(n):\n",
    "    if n == 1:\n",
    "        return False\n",
    "    return trial_division(n)[0] == n\n",
    "\n",
    "assert(fermat_test_sicp(3))\n",
    "assert(fermat_test_sicp(107))\n",
    "assert(fermat_test_sicp(108) == False)\n",
    "assert(fermat_test_sicp(109))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n: 26241 count_total: 52481 count_primes: 5248 ratio: 0.09999809454850327\n"
     ]
    }
   ],
   "source": [
    "def get_solution():\n",
    "    n = 1\n",
    "    count_primes = 0\n",
    "    count_total = 0\n",
    "    while True:\n",
    "        for v in get_corner_values(n):\n",
    "            count_total += 1\n",
    "            if is_prime(v):\n",
    "                count_primes += 1\n",
    "        ratio = count_primes / count_total\n",
    "        if ratio != 0 and ratio < 0.10:\n",
    "            print(\"n: {} count_total: {} count_primes: {} ratio: {}\".format(n, count_total, count_primes, ratio))\n",
    "            return n\n",
    "        n += 2\n",
    "\n",
    "s = get_solution()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26241\n"
     ]
    }
   ],
   "source": [
    "print(s)\n",
    "assert(s == 26241)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I actually got different results for the Fermat test and for the prime test which relies on reliable computation. I will actually try to solve the problem with ramdonized Fermat test now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n: 26641 count_total: 53281 count_primes: 5328 ratio: 0.09999812315834913\n",
      "26641\n",
      "Expected n is 26241.\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "def is_prime(n):\n",
    "    if n == 1:\n",
    "        return False\n",
    "    for _ in range(10):\n",
    "        a = n - random.randint(1, n - 1)\n",
    "        if not expmod(a, n, n) == a:\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "print(get_solution())\n",
    "print(\"Expected n is 26241.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Something seems to be off with my Fermat test... \n",
    "\n",
    "Seems like there are systematic errors with the Fermat tests. Certain primes cannot be deteced.\n",
    "\n",
    "Try this algorithm instead [https://en.wikipedia.org/wiki/Miller–Rabin_primality_test](https://en.wikipedia.org/wiki/Miller–Rabin_primality_test)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n: 26241 count_total: 52481 count_primes: 5248 ratio: 0.09999809454850327\n",
      "26241\n",
      "For 561 sicp says True and real test says False.\n",
      "For 1105 sicp says True and real test says False.\n",
      "For 1729 sicp says True and real test says False.\n",
      "For 2465 sicp says True and real test says False.\n",
      "For 2821 sicp says True and real test says False.\n",
      "For 6601 sicp says True and real test says False.\n"
     ]
    }
   ],
   "source": [
    "carmichael_number = [561, 1105, 1729, 2465, 2821, 6601]\n",
    "\n",
    "def little_fermat(a, p): # SICP uses this\n",
    "    return (a ** p) % p == a % p\n",
    "\n",
    "def little_fermat_simple(a, p): # Wikipedia uses this\n",
    "    assert(a < p)\n",
    "    return (a ** (p - 1)) % p == 1\n",
    "\n",
    "def expmod(base, exp, m):\n",
    "    if exp == 0:\n",
    "        return 1\n",
    "    if (exp % 2 == 0):\n",
    "        return (expmod(base, exp // 2, m) ** 2 % m)\n",
    "    return (base * expmod(base, exp - 1, m) % m)\n",
    "\n",
    "def fermat_test(n):\n",
    "    ps = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 39]\n",
    "    if n <= ps[-1]:\n",
    "        return n in ps\n",
    "    for a in ps:\n",
    "        #if not little_fermat_simple(a, n):\n",
    "        if expmod(a, n - 1, n) != 1:\n",
    "            #print(\"Fermat witness for {} is {}.\".format(n + 1, a))\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "def is_prime(n):\n",
    "    return fermat_test(n)\n",
    "\n",
    "print(get_solution())\n",
    "\n",
    "for n in carmichael_number:\n",
    "    r_sicp = fermat_test_sicp(n)\n",
    "    r = fermat_test(n)\n",
    "    print(\"For {} sicp says {} and real test says {}.\".format(n, r_sicp, r))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "completion_date": "Sat, 29 Dec 2018, 09:00",
  "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": [
   "prime",
   "fermat",
   "spiral",
   "diagonal",
   "todo"
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 2
}