迭代动力系统

整数值迭代函数的有向图

3N上的立方和

数字153有一个奇特的性质。

设3N={3,6,9,12,…}为3的正倍数集合。定义一个迭代过程f:3N->3N如下：对于给定的n，取n的每一位数字（以10为基数），将其立方，然后求立方和得到f(n)。

当这个过程重复进行时，得到的序列n, f(n), f(f(n)),…在有限次迭代后终止于153（因为153 = 1**3 + 5**3 + 3**3，过程结束）。

在离散动力系统的语言中，153是迭代映射f在集合3N上的全局吸引子。

例如：取数字108

f(108) = 1**3 + 0**3 + 8**3 = 513

并且

f(513) = 5**3 + 1**3 + 3**3 = 153

所以，从108开始，我们在两次迭代后到达153，表示为：

108->513->153

计算3N中所有轨道的上限为10**5，发现吸引子153在最多次数为14次迭代后达到。在这个代码中，我们展示了对于所有小于10,000的整数（在3N中），最多需要13次迭代。

需要13次迭代才能到达153的最小数字是177，即，

177->687->1071->345->216->225->141->66->432->99->1458->702->351->153

由此产生的大型有向图对于测试网络软件很有用。

一般问题

给定数字n、幂p和基数b，定义F(n; p, b)为n的各位数字（以b为基数）的p次幂之和。上述例子对应于f(n)=F(n; 3,10)，下面F(n; p, b)实现为函数powersum(n,p,b)。由上述映射n:->f(n)定义的迭代动力系统（在3N上）收敛到一个固定点；153。将映射应用于所有正整数N，导致一个具有5个固定点的离散动力过程：1, 153, 370, 371, 407。模3这些数字是1, 0, 1, 2, 2。上述函数f具有将3的倍数映射到另一个3的倍数的附加性质；即它在子集3N上是不变的。

数字的平方（以10为基数）导致循环和单个固定点1。即，从某个点开始，过程开始重复自身。

关键词：“循环数字不变量”，“自恋数”，“快乐数”

3n+1问题

与离散动力系统相关的数学娱乐有着丰富的历史。最著名的是科拉茨3n+1问题。参见下面的函数collatz_problem_digraph。科拉茨猜想——每个轨道在有限时间内返回到固定点1——尚未被证明。即使是伟大的保罗·埃尔德什也说“数学还没有准备好解决这样的问题”，并为它的解决方案提供了500美元。

关键词：“3n+1”，“3x+1”，“科拉茨问题”，“斯威特猜想”

Building cubing_153_digraph(10000)
Resulting digraph has 10000 nodes and 10000  edges
Shortest path from 177 to 153 is:
[177, 687, 1071, 345, 216, 225, 141, 66, 432, 99, 1458, 702, 351, 153]
fixed points are []

import networkx as nx

nmax = 10000
p = 3


def digitsrep(n, b=10):
    """返回以基数b表示的非负整数n的数字列表。
"""

    if n <= 0:
        return [0]

    dlist = []
    while n > 0:
        # 前置下一个次低位数字
        dlist = [n % b] + dlist
        # Floor-division
        n = n // b
    return dlist


def powersum(n, p, b=10):
    """返回数字 n（以 b 为基数）的 p 次幂的各位数字之和。"""
    dlist = digitsrep(n, b)
    sum = 0
    for k in dlist:
        sum += k**p
    return sum


def attractor153_graph(n, p, multiple=3, b=10):
    """返回 powersum(n,3,10) 迭代的有向图。"""
    G = nx.DiGraph()
    for k in range(1, n + 1):
        if k % multiple == 0 and k not in G:
            k1 = k
            knext = powersum(k1, p, b)
            while k1 != knext:
                G.add_edge(k1, knext)
                k1 = knext
                knext = powersum(k1, p, b)
    return G


def squaring_cycle_graph_old(n, b=10):
    """返回 powersum(n,2,10) 迭代的有向图。"""
    G = nx.DiGraph()
    for k in range(1, n + 1):
        k1 = k
        G.add_node(k1)  # 当k1等于knext时，至少添加节点
        knext = powersum(k1, 2, b)
        G.add_edge(k1, knext)
        while k1 != knext:  # stop if fixed point
            k1 = knext
            knext = powersum(k1, 2, b)
            G.add_edge(k1, knext)
            if G.out_degree(knext) >= 1:
                # knext 已经进行了迭代输入和输出
                break
    return G


def sum_of_digits_graph(nmax, b=10):
    def f(n):
        return powersum(n, 1, b)

    return discrete_dynamics_digraph(nmax, f)


def squaring_cycle_digraph(nmax, b=10):
    def f(n):
        return powersum(n, 2, b)

    return discrete_dynamics_digraph(nmax, f)


def cubing_153_digraph(nmax):
    def f(n):
        return powersum(n, 3, 10)

    return discrete_dynamics_digraph(nmax, f)


def discrete_dynamics_digraph(nmax, f, itermax=50000):
    G = nx.DiGraph()
    for k in range(1, nmax + 1):
        kold = k
        G.add_node(kold)
        knew = f(kold)
        G.add_edge(kold, knew)
        while kold != knew and kold << itermax:
            # 迭代直到达到固定点或超过itermax
            kold = knew
            knew = f(kold)
            G.add_edge(kold, knew)
            if G.out_degree(knew) >= 1:
                # 已经迭代进出过
                break
    return G


def collatz_problem_digraph(nmax):
    def f(n):
        if n % 2 == 0:
            return n // 2
        else:
            return 3 * n + 1

    return discrete_dynamics_digraph(nmax, f)


def fixed_points(G):
    """返回由有向图G表示的离散动力系统的不动点列表。
"""
    return [n for n in G if G.out_degree(n) == 0]


nmax = 10000
print(f"Building cubing_153_digraph({nmax})")
G = cubing_153_digraph(nmax)
print("Resulting digraph has", len(G), "nodes and", G.size(), " edges")
print("Shortest path from 177 to 153 is:")
print(nx.shortest_path(G, 177, 153))
print(f"fixed points are {fixed_points(G)}")