AIとともにコラッツ予想について考えてみた(2)

概要

前回の続き
$$
T_{n+1}(x) =
\begin{cases}
3x + 1 (xが奇数) \\
x/2 (xが偶数)
\end{cases}
$$
これがコラッツ予想の元になったもので標準コラッツと呼ばれる。
今回はこれらをいくつか可視化していき考察していく。
まず、$3n+d$と1を$d$と一般化した形について考え、次に巨大数にした場合の巨視的な挙動について考える。
最後にはそれに使用したjuliaのコードをつける。

3x+d系のサイクル

まず、$3x+d$と、$1$を$d$とした場合について見ていく
$3x+1$の標準コラッツは以下すべての値が1に落ちているのがわかる

次に$3x+3$が以下であり、多くの初期値、特に大きくなると$3 \rightarrow 12 \rightarrow 6 \rightarrow 3$の3stepの巡回サイクルに入っている。

更に$3x+5$が以下であり複数の巡回サイクルが出てきていて値が大きくても1まで落ちる初期値もある。

最後に$3x+101$とすると以下のようになり、全体が上にシフトしていることが分る。

巡回サイクルの数

巡回サイクルの数を数えると以下のようになった。

これは初期値は1000までしか振ってないので全巡回軌道を数え上げたわけではない点には注意が必要である。
ここで自明なサイクル一つしか持たないものは$3^n$つまり$1, 3, 9, 81, 243, 729$しかありません。
これは$2$と$3^n$の間に強い素性(互いに強い素)があるためと考えられます。

次に見つかった巡回サイクルのdの値と最小値(Cycle Value)の関係は以下のようになります。
ここで色は巡回サイクルのステップ数の対数($\log(N_{loop})$)である。
$y=x$の直線に乗っているのは$d \rightarrow 3d + d = 4d \rightarrow 2d \rightarrow d$の3stepサイクルであり、標準コラッツの$1 \rightarrow 4 \rightarrow 2 \rightarrow 1$に対応するものである。

またそれ以外にも線形に乗っているように見えるものもある。
それを確かめるために最小値をdで割ったものをプロットすると以下のようになる。
このようになり
$$
M_{min}/d \sim 3.8, 4.6, 10
$$
あたりに巡回サイクルが現れるような系列が見える。
それ以外に$M_{min}/d$の小さい、特に$d<1$にランダムに巡回軌道が現れるように見える。
特にそういう軌道は長い周期を持つように見える。

dが巨大数の場合

次に巨大数の場合どうなるかについて実験的に探ってみる。
まず、標準コラッツと似た構造を持つ$d = 3^n$の場合を考える。
その時の$n$とステップ数の関係は以下のようになった。

いくつかの初期値で試しており、(初期値)$ = d$となるところで振る舞いを変え、(初期値)$< d$ではnに対してある程度、線形に動いているように見える。

サイクルの最大値は以下のようになった。

(初期値)$ > d$では初期値付近に張り付き、その後はdに張り付くように見える。

これを考察するために$d=3^{150}$でのサイクルの動きをみると次のようになった。

(初期値)$ \ll d$の場合、一度$d$付近まで上がってそこかで巡回サイクルに入って止まっている。
これは、サイクルが落ちる$1/2$の力よりも、上に持ち上げる$d = 3^n$の力が強いからと考えられる。
(初期値)$ \gg d$の場合、指数的減少している(y軸が対数であることに注意)。
これはドリフトなどと呼んだ、サイクルが落ちる力、$\exp(-3/4)$が効いていると考えられる。
そして$d$付近ではそれら二つの力がせめぎあいその付近を動いていると考えられる。

まとめ

今回はいく$3n + d$型のコラッツ系ををいくつか可視化してみてどのようになるか見てみた。
その中で$d$を巨大数、つまり巨視的に見た場合は$\exp(-3/4)$のドリフトの効果と$d$が値を押し上げる力があり$d$付近ではちょうど釣りあっているという構造を確認した。

また小さい$d$の領域では巡回サイクルについてサーチし、列挙した。
$M_{min}/d = 1$には自明な巡回サイクルが見えた、当たり前だが
それ以外に

• $M_{min}/d = 3.8$付近
• $M_{min}/d = 4.6$付近
• $M_{min}/d = 10$付近

に非連続ながら線になるような、ある程度周期が安定した巡回サイクルが見られた。

それ以外に、$d$が小さい領域にランダムに現れるような巡回サイクルが見られ、これらは時に大きな周期を持つように見えた。

次に考えられること

• 小さな$d$の領域で構造を持つであろう巡回サイクル(M_{min}/d = 3.8, 4.6, 10等)を列挙、中身を具体的に解析することでより精密な構造を探る、そのうえで巨大数の領域でどうなるかのサーチ等。
• $M_{min}/d<1$でのランダムに見える巡回サイクルについても列挙、法則性を探る等

が考えられる。

使ったコードなど

サイクル描画

plot.jl

import Colors

using Printf
using Plots
using Colors

const q = 3
const p = 243
const k = 2

# const i::Int = 27
const max_value::Int = 10^3
const interval::Int = 1

const next = x -> x % k == 0 ? x÷k : q*x + p
const reverse = T -> (T-p) % q == 0 ? (T-p)÷q : k*T

function main()
    @printf "%dx%+d\n" q p
    @printf "x/%d\n" k

    cycles = Vector{Vector{Int}}()
    @time begin
        for i in 1:interval:max_value
            push!(cycles, calc_cycle(i))
        end
    end
    println("Calc Finish")

    @time begin
        plt = plot(1: length(cycles[1]), cycles[1],
                   color = HSV(360/max_value, 1.0, 1.0),
                   xlabel = "Step", ylabel = "Value", yscale = :log10, legend = false, marker = :circle
                   )
        for i in 2:length(cycles)
            plot!(plt, 1: length(cycles[i]), cycles[i], color = HSV(i*360/length(cycles), 1.0, 1.0), marker = :circle)
        end
        display(plt)
    end
    savefig(plt, "pic/3n+243.png")
    readline()
end

function calc_cycle(i::Int)::Vector{Int}
    cycle = Int[ i ]
    local x = i
    while true
        x = next(x)
#        println(x)
        if x in cycle
            push!(cycle, x)
            println("!!! Non-Trivial Cycle : ")
            break
        end

        push!(cycle, x)
        if x == 1
#            println("Trivial Cycle")
            break
        elseif x > 10^17
            println("!!! n is overflow !!!")
            break
        end
    end
    return cycle
end

main()

巡回サイクルカウント、列挙

n_cycle.jl

import Colors

using Printf
using Plots
using Colors

const q::Int = 3
const k::Int = 2
const max_d::Int = 10^3+1

const max_init::Int = 100

const generate_next = (q::Int, p::Int, k::Int)-> n -> n % k ==0 ? n÷k : q*n + p
const generate_reverse = (q::Int, p::Int, k::Int)-> T -> (T-p) % k == 0 ? (T-p)÷q : k*T

function my_savefig(plt, name)
    display(plt)
    println("output figure : ", name)
    savefig(name)
    readline()
end

function plot_scatter(cycles_arr)
    xs::Vector{Int} = []
    ys::Vector{Int} = []
    y2s::Vector{Float64} = []
    counts::Vector{Int} = []
    for (d_idx, cycles) in enumerate(cycles_arr)
        mins::Vector{Int} = []
        for cycle in cycles
            loop = get_loop(cycle)
            if min in mins
            else
                push!(mins, minimum(loop))
                push!(xs, 2*d_idx-1)
                push!(ys, minimum(loop))
                push!(y2s, minimum(loop)/(2*d_idx-1))
                push!(counts, length(loop))
            end
        end
    end

    plt = scatter(xs, ys, marker_z = log10.(counts),
                  title = "Found Loop Cycle", xlabel = "d", ylabel = "Cycle Value",
                  markercolor = :viridis, markerstrokewidth = 0, colorbar = true,
                  label = false)
    my_savefig(plt, "pic/cycle_sscatter_1000.png")

    plt = scatter(xs, y2s, marker_z = log10.(counts),
                  title = "Found Loop Cycle", xlabel = "d", ylabel = "(Cycle Value)/d",
                  markercolor = :viridis, markerstrokewidth = 0, colorbar = true,
                  label = false)
    my_savefig(plt, "pic/cycle_sscatter_d_1000.png")
end

function plot_n_cycle(cycles_arr)
    ds = Int[]
    loop_mins = Vector{Vector{Int}}()
    for (d_idx, cycles) in enumerate(cycles_arr)
        push!(ds, 2*d_idx-1)
        mins::Vector{Int} = []
        for cycle in cycles
            push!(mins, get_loop_min(cycle))
        end
        push!(loop_mins, mins)
    end

    println(ds)
    println(map(x -> length(unique(x)), loop_mins))
    plt = plot(ds, map(x -> length(unique(x)), loop_mins), xlabel = "d", ylabel = "n-cycle", label = false, marker = :circle)
    my_savefig(plt, "pic/n_cycle.png")
end

function main()
    cycles_arr = Vector{Vector{Vector{Int}}}()
    for d in 1:2:max_d
        println("Calc d = ", d, " START")
        cycles = Vector{Vector{Int}}()
        for i in 1:max_init
            push!(cycles, calc_cycle(d, i))
        end
        push!(cycles_arr, cycles)
    end

    plot_scatter(cycles_arr)
end

function get_loop(cycle)
    if cycle[end] == 1
        return cycle[end:end]
    else
        idx = findlast(==(cycle[end]), cycle[1:end-1])
        return cycle[idx:end-1]
    end
end

function calc_cycle(i, initial_value)
    local next = generate_next(q, i, k)

    local x = initial_value
    local cycle = Int[ x ]
    while true
        x = next(x)
        if x in cycle
            push!(cycle, x)
            #                println("!!! Non-Trivial Cycle : ", length(cycle))
            break
        end

        push!(cycle, x)
    end
    return cycle
end

main()

巨大数の実験

import Colors

using Printf
using Plots
using Colors

const q::BigInt = 3
const k::BigInt = 2

const N::BigInt = 500
const initial_values::Vector{BigInt} = [ 0, 100, 200, 300 ]

const generate_next = (q::BigInt, p::BigInt, k::BigInt)-> n -> n % k ==0 ? n÷k : q*n + p
const generate_reverse = (q::BigInt, p::BigInt, k::BigInt)-> T -> (T-p) % k == 0 ? (T-p)÷q : k*T

function my_savefig(plt, name)
    display(plt)
    println("output figure : ", name)
    savefig(name)
    readline()
end

function plot_step(cycles_arr)
    plt = plot(1: length(cycles_arr[1]), map(x-> length(x), cycles_arr[1]),
                label = "Initail value:1", xlabel = "n", ylabel = "Step")
    for i in 2:length(cycles_arr)
        plot!(plt, 1: length(cycles_arr[i]), map(x-> length(x), cycles_arr[i]), label = "Initial value:3^$(initial_values[i])")
    end
    my_savefig(plt, "pic/3n_step.png")
end

function plot_max(cycles_arr)
    plt = plot(1: length(cycles_arr[1]), map(x-> maximum(x), cycles_arr[1]),
#               label = "Initail value:1", xlabel = "n", ylabel = "log(Max)")  # Log Scale
               label = "Initail value:1", xlabel = "n", ylabel = "Max", yscale = :log10)
    for i in 2:length(cycles_arr)
        plot!(plt, 1: length(cycles_arr[i]), map(x-> maximum(x), cycles_arr[i]), label = "Initial value:3^$(initial_values[i])")
    end
    my_savefig(plt, "pic/3n_max.png")
end

function plot_cycle(cycles_arr, n)
    plt = plot(1: length(cycles_arr[1][n]), cycles_arr[1][n], label = "Inital values:1", title = "3n + 3^$(n)",
               xlabel = "Step", ylabel = "Value", yscale = :log10)
    for i in 2:length(cycles_arr)
        plot!(plt, 1: length(cycles_arr[i][n]), cycles_arr[i][n],
              xlabel = "Step", ylabel = "Value", label = "Initial value:3^$(initial_values[i])")
    end
    my_savefig(plt, "pic/3n_cyvle.png")
end

function main()
    cycles_arr = Vector{Vector{Vector{BigInt}}}()
    for initial_value in initial_values
        cycles = Vector{Vector{BigInt}}()
        println("Inital value START 3^", initial_value)
        for i in 1:N
            push!(cycles, calc_cycle(3^i, 3^initial_value))
        end
        push!(cycles_arr, cycles)
    end

    plot_max(cycles_arr)
    plot_step(cycles_arr)
    plot_cycle(cycles_arr, 150)
end

function calc_cycle(i, initial_value)
    local next = generate_next(q, i, k)

    local x = initial_value
    local cycle = BigInt[ x ]
    while true
        x = next(x)
        if x in cycle
            push!(cycle, x)
            #                println("!!! Non-Trivial Cycle : ", length(cycle))
            break
        end

        push!(cycle, x)
        if x == 1
            println("1-Cycle")
            break
        end
    end
    return cycle
end

main()