euler/ipython/EulerProblem012.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Euler Problem 12\n",
    "\n",
    "The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be $1 + 2 + 3 + 4 + 5 + 6 + 7 = 28$. The first ten terms would be:\n",
    "\n",
    "$1, 3, 6, 10, 15, 21, 28, 36, 45, 55, \\dots$\n",
    "\n",
    "Let us list the factors of the first seven triangle numbers:\n",
    "\n",
    "~~~\n",
    " 1: 1\n",
    " 3: 1,3\n",
    " 6: 1,2,3,6\n",
    "10: 1,2,5,10\n",
    "15: 1,3,5,15\n",
    "21: 1,3,7,21\n",
    "28: 1,2,4,7,14,28\n",
    "~~~\n",
    "\n",
    "We can see that 28 is the first triangle number to have over five divisors.\n",
    "\n",
    "What is the value of the first triangle number to have over five hundred divisors?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try brute forcing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_divisors(n):\n",
    "    return [i for i in range(1, int(n/2) + 1) if n % i == 0] + [n]\n",
    "\n",
    "assert(get_divisors(28) == [1,2,4,7,14,28])\n",
    "\n",
    "def triangle_number_generator_function():\n",
    "    c = 0\n",
    "    for i in range(1, 100000000000):\n",
    "        c += i\n",
    "        yield c\n",
    "\n",
    "#ts = triangle_number_generator_function()\n",
    "#for t in ts:\n",
    "#    if len(get_divisors(t)) > 500:\n",
    "#        print(t)\n",
    "#        break\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That failed miserably. We have to come up with something more efficient. I looked at my old Haskell solution which looks like this:\n",
    "\n",
    "~~~\n",
    "divisor_count' x = product . map (succ . length) . group $ prim_factors x\n",
    "~~~\n",
    "\n",
    "Add first I did not understand what it does at all. But the algorithm is as follows:\n",
    "\n",
    "~~~\n",
    "# get prime factors, for example for 28\n",
    "2 * 2 * 7\n",
    "# group the primes\n",
    "[[2, 2], [7]]\n",
    "# get the length of each group and increment it by one\n",
    "[3, 2]\n",
    "# get the product of all values\n",
    "6\n",
    "~~~\n",
    "\n",
    "Honestly, I have no idea why this gives as the  number of divisors. I will implement the solution and then try to come up with an explanation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def is_prime(n, smaller_primes):\n",
    "    for s in smaller_primes:\n",
    "        if n % s == 0:\n",
    "            return False\n",
    "        if s * s > n:\n",
    "            return True\n",
    "    return True\n",
    "\n",
    "def prime_generator_function():\n",
    "    primes = [2, 3, 5, 7]\n",
    "    for p in primes:\n",
    "        yield p\n",
    "    while True:\n",
    "        p += 2\n",
    "        if is_prime(p, primes):\n",
    "            primes.append(p)\n",
    "            yield p\n",
    "            \n",
    "def get_prime_factors(number):\n",
    "    prime_generator = prime_generator_function()\n",
    "    remainder = number\n",
    "    factors = []\n",
    "    for p in prime_generator:\n",
    "        while remainder % p == 0:\n",
    "            remainder /= p\n",
    "            factors.append(p)\n",
    "        if remainder == 1 or p * p > number:\n",
    "            break\n",
    "    if remainder != 1:\n",
    "        factors.append(remainder)\n",
    "    return factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This are the prime numbers related functions we already now. Now we implement a group function, the product function we already know, and based on that the algorithm to get the number of divisors mentioned above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def group(xs):\n",
    "    from functools import reduce\n",
    "    def f(xss, x):\n",
    "        if xss and x in xss[-1]:\n",
    "            xss[-1].append(x)\n",
    "        else:\n",
    "            xss.append([x])\n",
    "        return xss\n",
    "    return reduce(f, xs, [])\n",
    "\n",
    "def product(xs):\n",
    "    from functools import reduce\n",
    "    from operator import mul\n",
    "    return reduce(mul, xs, 1)\n",
    "\n",
    "assert(group([1, 1, 3, 2, 2]) == [[1, 1], [3], [2, 2]])\n",
    "\n",
    "def get_number_of_divisors(n):\n",
    "    # get prime factors, for example for 28\n",
    "    ps = get_prime_factors(n)\n",
    "    # group the primes\n",
    "    ps = group(ps)\n",
    "    # get the length of each group and increment it by one\n",
    "    ps = map(lambda x: len(x) + 1, ps)\n",
    "    # get the product of all values\n",
    "    return product(ps)\n",
    "            \n",
    "assert(get_number_of_divisors(28) == 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to do another brute force attempt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "76576500\n"
     ]
    }
   ],
   "source": [
    "ts = triangle_number_generator_function()\n",
    "for t in ts:\n",
    "    if get_number_of_divisors(t) > 500:\n",
    "        print(t)\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the only question is why this crazy algorithm works. Okay, I got it with the help of [this](https://www.math.upenn.edu/~deturck/m170/wk2/numdivisors.html) page. The problem is actually an instance of the multiplication principle for counting things. Each prime can be used to calculate (n + 1) other divisors. The incrementation by one is required because we can also choose to not use a certain prime. Of course, to get the potential combinations we have to get the product of the potential divisors for all primes. The web page explains it better. The only question is whether I came up with this algorithm myself back then."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "completion_date": "Sun, 31 Aug 2014, 17:07",
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.4"
  },
  "tags": [
   "triangular",
   "number",
   "divisors",
   "factorization",
   "prime",
   "group"
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 0
}