{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Iterated Dynamical Systems\n\n\nDigraphs from Integer-valued Iterated Functions\n\nSums of cubes on 3N\n-------------------\n\nThe number 153 has a curious property.\n\nLet 3N={3,6,9,12,...} be the set of positive multiples of 3.  Define an\niterative process f:3N->3N as follows: for a given n, take each digit\nof n (in base 10), cube it and then sum the cubes to obtain f(n).\n\nWhen this process is repeated, the resulting series n, f(n), f(f(n)),...\nterminate in 153 after a finite number of iterations (the process ends\nbecause 153 = 1**3 + 5**3 + 3**3).\n\nIn the language of discrete dynamical systems, 153 is the global\nattractor for the iterated map f restricted to the set 3N.\n\nFor example: take the number 108\n\nf(108) = 1**3 + 0**3 + 8**3 = 513\n\nand\n\nf(513) = 5**3 + 1**3 + 3**3 = 153\n\nSo, starting at 108 we reach 153 in two iterations,\nrepresented as:\n\n108->513->153\n\nComputing all orbits of 3N up to 10**5 reveals that the attractor\n153 is reached in a maximum of 14 iterations. In this code we\nshow that 13 cycles is the maximum required for all integers (in 3N)\nless than 10,000.\n\nThe smallest number that requires 13 iterations to reach 153, is 177, i.e.,\n\n177->687->1071->345->216->225->141->66->432->99->1458->702->351->153\n\nThe resulting large digraphs are useful for testing network software.\n\nThe general problem\n-------------------\n\nGiven numbers n, a power p and base b, define F(n; p, b) as the sum of\nthe digits of n (in base b) raised to the power p. The above example\ncorresponds to f(n)=F(n; 3,10), and below F(n; p, b) is implemented as\nthe function powersum(n,p,b). The iterative dynamical system defined by\nthe mapping n:->f(n) above (over 3N) converges to a single fixed point;\n153. Applying the map to all positive integers N, leads to a discrete\ndynamical process with 5 fixed points: 1, 153, 370, 371, 407. Modulo 3\nthose numbers are 1, 0, 1, 2, 2. The function f above has the added\nproperty that it maps a multiple of 3 to another multiple of 3; i.e. it\nis invariant on the subset 3N.\n\n\nThe squaring of digits (in base 10) result in cycles and the\nsingle fixed point 1. I.e., from a certain point on, the process\nstarts repeating itself.\n\nkeywords: \"Recurring Digital Invariant\", \"Narcissistic Number\",\n\"Happy Number\"\n\nThe 3n+1 problem\n----------------\n\nThere is a rich history of mathematical recreations\nassociated with discrete dynamical systems.  The most famous\nis the Collatz 3n+1 problem. See the function\ncollatz_problem_digraph below. The Collatz conjecture\n--- that every orbit returrns to the fixed point 1 in finite time\n--- is still unproven. Even the great Paul Erdos said \"Mathematics\nis not yet ready for such problems\", and offered $500\nfor its solution.\n\nkeywords: \"3n+1\", \"3x+1\", \"Collatz problem\", \"Thwaite's conjecture\"\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "import networkx as nx\n\nnmax = 10000\np = 3\n\n\ndef digitsrep(n, b=10):\n    \"\"\"Return list of digits comprising n represented in base b.\n    n must be a nonnegative integer\"\"\"\n\n    if n <= 0:\n        return [0]\n\n    dlist = []\n    while (n > 0):\n        # Prepend next least-significant digit\n        dlist = [n % b] + dlist\n        # Floor-division\n        n = n // b\n    return dlist\n\n\ndef powersum(n, p, b=10):\n    \"\"\"Return sum of digits of n (in base b) raised to the power p.\"\"\"\n    dlist = digitsrep(n, b)\n    sum = 0\n    for k in dlist:\n        sum += k**p\n    return sum\n\n\ndef attractor153_graph(n, p, multiple=3, b=10):\n    \"\"\"Return digraph of iterations of powersum(n,3,10).\"\"\"\n    G = nx.DiGraph()\n    for k in range(1, n + 1):\n        if k % multiple == 0 and k not in G:\n            k1 = k\n            knext = powersum(k1, p, b)\n            while k1 != knext:\n                G.add_edge(k1, knext)\n                k1 = knext\n                knext = powersum(k1, p, b)\n    return G\n\n\ndef squaring_cycle_graph_old(n, b=10):\n    \"\"\"Return digraph of iterations of powersum(n,2,10).\"\"\"\n    G = nx.DiGraph()\n    for k in range(1, n + 1):\n        k1 = k\n        G.add_node(k1)  # case k1==knext, at least add node\n        knext = powersum(k1, 2, b)\n        G.add_edge(k1, knext)\n        while k1 != knext:  # stop if fixed point\n            k1 = knext\n            knext = powersum(k1, 2, b)\n            G.add_edge(k1, knext)\n            if G.out_degree(knext) >= 1:\n                # knext has already been iterated in and out\n                break\n    return G\n\n\ndef sum_of_digits_graph(nmax, b=10):\n    def f(n): return powersum(n, 1, b)\n    return discrete_dynamics_digraph(nmax, f)\n\n\ndef squaring_cycle_digraph(nmax, b=10):\n    def f(n): return powersum(n, 2, b)\n    return discrete_dynamics_digraph(nmax, f)\n\n\ndef cubing_153_digraph(nmax):\n    def f(n): return powersum(n, 3, 10)\n    return discrete_dynamics_digraph(nmax, f)\n\n\ndef discrete_dynamics_digraph(nmax, f, itermax=50000):\n    G = nx.DiGraph()\n    for k in range(1, nmax + 1):\n        kold = k\n        G.add_node(kold)\n        knew = f(kold)\n        G.add_edge(kold, knew)\n        while kold != knew and kold << itermax:\n            # iterate until fixed point reached or itermax is exceeded\n            kold = knew\n            knew = f(kold)\n            G.add_edge(kold, knew)\n            if G.out_degree(knew) >= 1:\n                # knew has already been iterated in and out\n                break\n    return G\n\n\ndef collatz_problem_digraph(nmax):\n    def f(n):\n        if n % 2 == 0:\n            return n // 2\n        else:\n            return 3 * n + 1\n    return discrete_dynamics_digraph(nmax, f)\n\n\ndef fixed_points(G):\n    \"\"\"Return a list of fixed points for the discrete dynamical\n    system represented by the digraph G.\n    \"\"\"\n    return [n for n in G if G.out_degree(n) == 0]\n\n\nif __name__ == \"__main__\":\n    nmax = 10000\n    print(\"Building cubing_153_digraph(%d)\" % nmax)\n    G = cubing_153_digraph(nmax)\n    print(\"Resulting digraph has\", len(G), \"nodes and\",\n          G.size(), \" edges\")\n    print(\"Shortest path from 177 to 153 is:\")\n    print(nx.shortest_path(G, 177, 153))\n    print(\"fixed points are %s\" % fixed_points(G))"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }
  ], 
  "metadata": {
    "kernelspec": {
      "display_name": "Python 2", 
      "name": "python2", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "2.7.12", 
      "pygments_lexer": "ipython2", 
      "codemirror_mode": {
        "version": 2, 
        "name": "ipython"
      }
    }
  }
}