~ubuntu-branches/ubuntu/trusty/python-networkx/trusty-proposed

« back to all changes in this revision

Viewing changes to doc/examples/iterated_dynamical_systems.py

  • Committer: Bazaar Package Importer
  • Author(s): Cyril Brulebois
  • Date: 2009-02-28 13:36:24 UTC
  • mfrom: (1.2.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090228133624-9l5rwi1ftlzc7b0l
* Upload to unstable now that lenny is released (yay).
* Fix FTBFS with python-support 0.90.3: no longer rely on its internal
  behaviour, and xnow set tests/test.py executable right after “setup.py
  install” (Closes: #517065).
* Drop executable bits from bz2 files.
* Update Vcs-* fields: move from DPMT's svn to collab-maint's git.
* Remote DPMT from Uploaders, following Piotr Ożarowski's request.

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
"""
2
 
Digraphs from Integer-valued Iterated Functions
3
 
===============================================
4
 
 
5
 
 
6
 
Sums of cubes on 3N
7
 
-------------------
8
 
 
9
 
The number 153 has a curious property.
10
 
 
11
 
Let 3N={3,6,9,12,...} be the set of positive multiples of 3.  Define an
12
 
iterative process f:3N->3N as follows: for a given n, take each digit
13
 
of n (in base 10), cube it and then sum the cubes to obtain f(n).
14
 
 
15
 
When this process is repeated, the resulting series n, f(n), f(f(n)),...
16
 
terminate in 153 after a finite number of iterations (the process ends
17
 
because 153 = 1**3 + 5**3 + 3**3).
18
 
 
19
 
In the language of discrete dynamical systems, 153 is the global
20
 
attractor for the iterated map f restricted to the set 3N.
21
 
 
22
 
For example: take the number 108
23
 
 
24
 
f(108) = 1**3 + 0**3 + 8**3 = 513
25
 
 
26
 
and
27
 
 
28
 
f(513) = 5**3 + 1**3 + 3**3 = 153
29
 
 
30
 
So, starting at 108 we reach 153 in two iterations,
31
 
represented as:
32
 
 
33
 
108->513->153
34
 
 
35
 
Computing all orbits of 3N up to 10**5 reveals that the attractor
36
 
153 is reached in a maximum of 14 iterations. In this code we
37
 
show that 13 cycles is the maximum required for all integers (in 3N)
38
 
less than 10,000.
39
 
 
40
 
The smallest number that requires 13 iterations to reach 153, is 177, i.e.,
41
 
 
42
 
177->687->1071->345->216->225->141->66->432->99->1458->702->351->153
43
 
 
44
 
The resulting large digraphs are useful for testing network software. 
45
 
 
46
 
The general problem
47
 
-------------------
48
 
 
49
 
Given numbers n, a power p and base b, define F(n; p, b) as the sum of
50
 
the digits of n (in base b) raised to the power p. The above example
51
 
corresponds to f(n)=F(n; 3,10), and below F(n; p, b) is implemented as
52
 
the function powersum(n,p,b). The iterative dynamical system defined by
53
 
the mapping n:->f(n) above (over 3N) converges to a single fixed point;
54
 
153. Applying the map to all positive integers N, leads to a discrete
55
 
dynamical process with 5 fixed points: 1, 153, 370, 371, 407. Modulo 3
56
 
those numbers are 1, 0, 1, 2, 2. The function f above has the added
57
 
property that it maps a multiple of 3 to another multiple of 3; i.e. it
58
 
is invariant on the subset 3N.
59
 
 
60
 
 
61
 
The squaring of digits (in base 10) result in cycles and the
62
 
single fixed point 1. I.e., from a certain point on, the process
63
 
starts repeating itself.
64
 
 
65
 
keywords: "Recurring Digital Invariant", "Narcissistic Number",
66
 
"Happy Number"
67
 
 
68
 
The 3n+1 problem
69
 
----------------
70
 
 
71
 
There is a rich history of mathematical recreations
72
 
associated with discrete dynamical systems.  The most famous 
73
 
is the Collatz 3n+1 problem. See the function 
74
 
collatz_problem_digraph below. The Collatz conjecture
75
 
--- that every orbit returrns to the fixed point 1 in finite time 
76
 
--- is still unproven. Even the great Paul Erdos said "Mathematics
77
 
is not yet ready for such problems", and offered $500 
78
 
for its solution. 
79
 
 
80
 
keywords: "3n+1", "3x+1", "Collatz problem", "Thwaite's conjecture"
81
 
 
82
 
 
83
 
"""
84
 
from networkx import *
85
 
from math import *
86
 
 
87
 
 
88
 
nmax=10000
89
 
p=3
90
 
mach_eps=0.00000000001
91
 
 
92
 
def digitsrep(n,b=10):
93
 
    """Return list of digits comprising n represented in base b.
94
 
    n must be a nonnegative integer"""
95
 
 
96
 
    # very inefficient if you only work with base 10
97
 
    dlist=[]
98
 
    if n<=0:
99
 
        return [0]
100
 
    maxpow=int(floor( log(n)/log(b) + mach_eps ))
101
 
    pow=maxpow
102
 
    while pow>=0:
103
 
        x=int(floor(n/b**pow))
104
 
        dlist.append(x)
105
 
        n=n-x*b**pow
106
 
        pow=pow-1
107
 
    return dlist
108
 
 
109
 
def powersum(n,p,b=10):
110
 
    """Return sum of digits of n (in base b) raised to the power p."""
111
 
    dlist=digitsrep(n,b)
112
 
    sum=0
113
 
    for k in dlist:
114
 
        sum+=k**p
115
 
    return sum
116
 
 
117
 
def attractor153_graph(n,p,multiple=3,b=10):
118
 
    """Return digraph of iterations of powersum(n,3,10)."""
119
 
    G=DiGraph()
120
 
    for k in xrange(1,n+1):
121
 
        if k%multiple==0 and k not in G:
122
 
            k1=k
123
 
            knext=powersum(k1,p,b)
124
 
            while k1!=knext:
125
 
                G.add_edge(k1,knext)
126
 
                k1=knext
127
 
                knext=powersum(k1,p,b)
128
 
    return G
129
 
 
130
 
def squaring_cycle_graph_old(n,b=10):
131
 
    """Return digraph of iterations of powersum(n,2,10)."""
132
 
    G=DiGraph()
133
 
    for k in xrange(1,n+1):
134
 
        k1=k
135
 
        G.add_node(k1) # case k1==knext, at least add node
136
 
        knext=powersum(k1,2,b)
137
 
        G.add_edge(k1,knext)
138
 
        while k1!=knext: # stop if fixed point 
139
 
             k1=knext
140
 
             knext=powersum(k1,2,b)
141
 
             G.add_edge(k1,knext)
142
 
             if G.out_degree(knext) >=1: 
143
 
                 # knext has already been iterated in and out
144
 
                 break
145
 
    return G
146
 
 
147
 
def sum_of_digits_graph(nmax,b=10):
148
 
    def f(n): return powersum(n,1,b)
149
 
    return discrete_dynamics_digraph(nmax,f)
150
 
 
151
 
def squaring_cycle_digraph(nmax,b=10):
152
 
    def f(n): return powersum(n,2,b)
153
 
    return discrete_dynamics_digraph(nmax,f)
154
 
 
155
 
def cubing_153_digraph(nmax):
156
 
    def f(n): return powersum(n,3,10)
157
 
    return discrete_dynamics_digraph(nmax,f)
158
 
 
159
 
def discrete_dynamics_digraph(nmax,f,itermax=50000):
160
 
    G=DiGraph()
161
 
    for k in xrange(1,nmax+1):
162
 
        kold=k
163
 
        G.add_node(kold)
164
 
        knew=f(kold)
165
 
        G.add_edge(kold,knew)
166
 
        while kold!=knew and kold<<itermax: 
167
 
        # iterate until fixed point reached or itermax is exceeded
168
 
            kold=knew
169
 
            knew=f(kold)
170
 
            G.add_edge(kold,knew)
171
 
            if G.out_degree(knew) >=1: 
172
 
               # knew has already been iterated in and out
173
 
               break
174
 
    return G
175
 
 
176
 
def collatz_problem_digraph(nmax):
177
 
    def f(n):
178
 
        if n%2==0:
179
 
            return n/2
180
 
        else:
181
 
            return 3*n+1
182
 
    return discrete_dynamics_digraph(nmax,f)
183
 
 
184
 
def fixed_points(G):
185
 
    """Return a list of fixed points for the discrete dynamical 
186
 
    system represented by the digraph G.
187
 
    """
188
 
    return [n for n in G if G.out_degree(n)==0]
189
 
 
190
 
    
191
 
if __name__ == "__main__":
192
 
    nmax=10000
193
 
    print "Building cubing_153_digraph(%d)"% nmax
194
 
    G=cubing_153_digraph(nmax)
195
 
    print "Resulting digraph has", len(G), "nodes and", \
196
 
          G.size()," edges" 
197
 
    print "Shortest path from 177 to 153 is:"
198
 
    print shortest_path(G,177,153)
199
 
    print "fixed points are", fixed_points(G)
 
 
b'\\ No newline at end of file'