rustとjuliaを使ってコラッツ予想を100,000,000回プロットしてみた

今回は、コラッツ予想を1から100,000,000（1億）までの数に対して計算、そしてそれを散布図としてプロットしてみました。

コラッツ予想に関しては、過去にJuliaで取り組んだことがあるので、そちらをご覧ください。

今回は、コラッツ予想の計算をrustで実装して、そこで、1億回計算。その結果をCSVファイルに出力。その後、Juliaを使って、散布図として、描画していきます。

rustでコラッツ予想を計算してCSVファイル出力
juliaでコラッツ予想をプロットする
まとめ

rustでコラッツ予想を計算してCSVファイル出力

まずは、rustでコラッツ予想を計算していきます。

コラッツを計算する

コラッツ予想を計算する心臓部の関数をまずは作りましょう。

コラッツ予想は、ある整数ｎが偶数だったら2で割って、奇数だったら3倍して1を足すっていう操作を繰り返すとどんな数も1になる。っていうやつです。（雑に言うとね）

なので、これをそのままプログラムします。

fn calcollatz(n:u64)->u64{
    if n !=1{
        if n%2 ==0{
            return n/2;
        }else{
            return 3*n+1;
        }
    }
    return 1;
}

これで、偶数なら2で割る、奇数なら3倍して1を足す。っていう操作ができました。

次にこれが1になるまで繰り返す。っていう部分を作ります。

fn collatznums(mut n:u64)->Vec<u64>{
    let mut ans = Vec::new();
    while n != 1{
        ans.push(n);
        n = calcollatz(n);
    }
    ans.push(1);
    return ans;
}

これで、1に至るまでの数の動きがわかります。

それが、ベクターとして出力されるわけね。

例えばｎが3だったら

[3,10,5,16,8,4,2,1]

って感じだね。

今回は回数だけ使うから、こんなベクタで出力しなくても回数だけカウントすればいいんだけど、今後ね、中身も使うかもしれないから、一応、全行程を出力できるようにしておきました。

あとは、collatznums(n)を実行して、そのベクターの長さを求めればOK。

例えばこんな感じにすればいい

let ntimes=collatznums(10).len();

これで、回数が出てくるね。

データの生成

1から1億までのコラッツ予想の計算をしていきます。

出力結果は、一度ベクタに入れて保持します。

    let mut ntimes = Vec::new();
    for n in 1..=num{
        let co = collatznums(n);
        ntimes.push(co.len());
    }

回数をntimesとして、ベクタを作ってて、そこに1から1億まで計算して、その回数を求めた結果を追加していきます。

次に結果をCSVファイルに出力しますよ。

csvファイル出力

CSVファイルを扱うので、csvクレートを使います。

cargo add csv

ってやって、cargo.tomlファイルにcsvを追加しておきます。

let mut wtr = csv::Writer::from_path("collatz.csv").expect("msg");
wtr.write_record(&["num","ntimes"]).expect("msg");

まずは、csvライターを生成して、CSVファイルを作ります。

そのあと、CSVファイルにヘッダーを作ります。

for (i,val) in ntimes.into_iter().enumerate(){
    wtr.serialize((i+1,val)).expect("msg");
}

あとは、それぞれの行をCSVファイルに追加していきます。

インデックス番号＋１がちょうど求めた時の数になるので使います。

コマンドライン引数を受け取る。

今回は1億回やってるけど、何回って、好きな数で求められるようにしておきます。

数を変えるたびにコンパイルしなおすのは面倒なので、コマンドライン引数を受け取ってその回数計算するようにします。

let args:Vec<String> = env::args().collect();
let num :u64= args[1].parse().unwrap();

そんな複雑な処理はしないので、単純に受け取ってそれをu64にパースします。

これで、1から1億までコラッツ予想を計算するプログラムができました。

あとは、これをコンパイルして実行すると1から1億まで計算した結果が出力されます。

ちなみに1億回分の回数だけ記録したCSVファイルの容量はなんと！1ギガ以上あります(笑)

やばいね(笑)

さて、次はjuliaを使ってプロットしていきます。

juliaでコラッツ予想をプロットする

Juliaを使って、1から1億までのコラッツ予想の回数をプロットしていきます。

jupyterでやってもいいんだけど、さすがに1億までやるとね、時間かかるから、コマンドラインからやります。

プロットするためのJuliaプログラムを作ります。

using CSV,DataFrames,Plots,StatsPlots

function main()
    df = CSV.read("collatz.csv",DataFrame);
    describe(df) |> println
    @df df scatter(:num, :ntimes,alpha=0.25,markerstrokewidth=0,size=(1200,800));
    png("plot.png")
end

一応、describeを使って、どんなもんか簡単な集計値を計算して出力する機能も入れておきます。

scatterプロットは、pngファイルにして出力。

これをcollatzana.jlとして保存しておきます。

次にjuliaのREPLを起動して

julia> include("collatzana.jl")
julia> main()

今作ったスクリプトをインクルード、そして、関数の実行。

あとは、じっくり待ってください(笑)

僕は、20分ぐらいかかったかな？

こんな感じにプロットできました(笑)

ちなみに1億まで計算した結果の集計はこんな感じ

 Row │ variable  mean     min    median   max        nmissing  eltype   
     │ Symbol    Float64  Int64  Float64  Int64      Int64     DataType
─────┼──────────────────────────────────────────────────────────────────
   1 │ num         5.0e7      1    5.0e7  100000000         0  Int64
   2 │ ntimes    180.235      1  174.0          950         0  Int64

ntimesのところを見ると、最大値は950、中央値が174回、平均が180回ぐらい

ちなみに最頻値は175回。

最大の950回かかった数は、63,728,127でした。その経路はこんな感じ。

[63728127, 191184382, 95592191, 286776574, 143388287, 430164862, 215082431, 645247294, 322623647, 967870942, 483935471, 1451806414, 725903207, 2177709622, 1088854811, 3266564434, 1633282217, 4899846652, 2449923326, 1224961663, 3674884990, 1837442495, 5512327486, 2756163743, 8268491230, 4134245615, 12402736846, 6201368423, 18604105270, 9302052635, 27906157906, 13953078953, 41859236860, 20929618430, 10464809215, 31394427646, 15697213823, 47091641470, 23545820735, 70637462206, 35318731103, 105956193310, 52978096655, 158934289966, 79467144983, 238401434950, 119200717475, 357602152426, 178801076213, 536403228640, 268201614320, 134100807160, 67050403580, 33525201790, 16762600895, 50287802686, 25143901343, 75431704030, 37715852015, 113147556046, 56573778023, 169721334070, 84860667035, 254582001106, 127291000553, 381873001660, 190936500830, 95468250415, 286404751246, 143202375623, 429607126870, 214803563435, 644410690306, 322205345153, 966616035460, 483308017730, 241654008865, 724962026596, 362481013298, 181240506649, 543721519948, 271860759974, 135930379987, 407791139962, 203895569981, 611686709944, 305843354972, 152921677486, 76460838743, 229382516230, 114691258115, 344073774346, 172036887173, 516110661520, 258055330760, 129027665380, 64513832690, 32256916345, 96770749036, 48385374518, 24192687259, 72578061778, 36289030889, 108867092668, 54433546334, 27216773167, 81650319502, 40825159751, 122475479254, 61237739627, 183713218882, 91856609441, 275569828324, 137784914162, 68892457081, 206677371244, 103338685622, 51669342811, 155008028434, 77504014217, 232512042652, 116256021326, 58128010663, 174384031990, 87192015995, 261576047986, 130788023993, 392364071980, 196182035990, 98091017995, 294273053986, 147136526993, 441409580980, 220704790490, 110352395245, 331057185736, 165528592868, 82764296434, 41382148217, 124146444652, 62073222326, 31036611163, 93109833490, 46554916745, 139664750236, 69832375118, 34916187559, 104748562678, 52374281339, 157122844018, 78561422009, 235684266028, 117842133014, 58921066507, 176763199522, 88381599761, 265144799284, 132572399642, 66286199821, 198858599464, 99429299732, 49714649866, 24857324933, 74571974800, 37285987400, 18642993700, 9321496850, 4660748425, 13982245276, 6991122638, 3495561319, 10486683958, 5243341979, 15730025938, 7865012969, 23595038908, 11797519454, 5898759727, 17696279182, 8848139591, 26544418774, 13272209387, 39816628162, 19908314081, 59724942244, 29862471122, 14931235561, 44793706684, 22396853342, 11198426671, 33595280014, 16797640007, 50392920022, 25196460011, 75589380034, 37794690017, 113384070052, 56692035026, 28346017513, 85038052540, 42519026270, 21259513135, 63778539406, 31889269703, 95667809110, 47833904555, 143501713666, 71750856833, 215252570500, 107626285250, 53813142625, 161439427876, 80719713938, 40359856969, 121079570908, 60539785454, 30269892727, 90809678182, 45404839091, 136214517274, 68107258637, 204321775912, 102160887956, 51080443978, 25540221989, 76620665968, 38310332984, 19155166492, 9577583246, 4788791623, 14366374870, 7183187435, 21549562306, 10774781153, 32324343460, 16162171730, 8081085865, 24243257596, 12121628798, 6060814399, 18182443198, 9091221599, 27273664798, 13636832399, 40910497198, 20455248599, 61365745798, 30682872899, 92048618698, 46024309349, 138072928048, 69036464024, 34518232012, 17259116006, 8629558003, 25888674010, 12944337005, 38833011016, 19416505508, 9708252754, 4854126377, 14562379132, 7281189566, 3640594783, 10921784350, 5460892175, 16382676526, 8191338263, 24574014790, 12287007395, 36861022186, 18430511093, 55291533280, 27645766640, 13822883320, 6911441660, 3455720830, 1727860415, 5183581246, 2591790623, 7775371870, 3887685935, 11663057806, 5831528903, 17494586710, 8747293355, 26241880066, 13120940033, 39362820100, 19681410050, 9840705025, 29522115076, 14761057538, 7380528769, 22141586308, 11070793154, 5535396577, 16606189732, 8303094866, 4151547433, 12454642300, 6227321150, 3113660575, 9340981726, 4670490863, 14011472590, 7005736295, 21017208886, 10508604443, 31525813330, 15762906665, 47288719996, 23644359998, 11822179999, 35466539998, 17733269999, 53199809998, 26599904999, 79799714998, 39899857499, 119699572498, 59849786249, 179549358748, 89774679374, 44887339687, 134662019062, 67331009531, 201993028594, 100996514297, 302989542892, 151494771446, 75747385723, 227242157170, 113621078585, 340863235756, 170431617878, 85215808939, 255647426818, 127823713409, 383471140228, 191735570114, 95867785057, 287603355172, 143801677586, 71900838793, 215702516380, 107851258190, 53925629095, 161776887286, 80888443643, 242665330930, 121332665465, 363997996396, 181998998198, 90999499099, 272998497298, 136499248649, 409497745948, 204748872974, 102374436487, 307123309462, 153561654731, 460684964194, 230342482097, 691027446292, 345513723146, 172756861573, 518270584720, 259135292360, 129567646180, 64783823090, 32391911545, 97175734636, 48587867318, 24293933659, 72881800978, 36440900489, 109322701468, 54661350734, 27330675367, 81992026102, 40996013051, 122988039154, 61494019577, 184482058732, 92241029366, 46120514683, 138361544050, 69180772025, 207542316076, 103771158038, 51885579019, 155656737058, 77828368529, 233485105588, 116742552794, 58371276397, 175113829192, 87556914596, 43778457298, 21889228649, 65667685948, 32833842974, 16416921487, 49250764462, 24625382231, 73876146694, 36938073347, 110814220042, 55407110021, 166221330064, 83110665032, 41555332516, 20777666258, 10388833129, 31166499388, 15583249694, 7791624847, 23374874542, 11687437271, 35062311814, 17531155907, 52593467722, 26296733861, 78890201584, 39445100792, 19722550396, 9861275198, 4930637599, 14791912798, 7395956399, 22187869198, 11093934599, 33281803798, 16640901899, 49922705698, 24961352849, 74884058548, 37442029274, 18721014637, 56163043912, 28081521956, 14040760978, 7020380489, 21061141468, 10530570734, 5265285367, 15795856102, 7897928051, 23693784154, 11846892077, 35540676232, 17770338116, 8885169058, 4442584529, 13327753588, 6663876794, 3331938397, 9995815192, 4997907596, 2498953798, 1249476899, 3748430698, 1874215349, 5622646048, 2811323024, 1405661512, 702830756, 351415378, 175707689, 527123068, 263561534, 131780767, 395342302, 197671151, 593013454, 296506727, 889520182, 444760091, 1334280274, 667140137, 2001420412, 1000710206, 500355103, 1501065310, 750532655, 2251597966, 1125798983, 3377396950, 1688698475, 5066095426, 2533047713, 7599143140, 3799571570, 1899785785, 5699357356, 2849678678, 1424839339, 4274518018, 2137259009, 6411777028, 3205888514, 1602944257, 4808832772, 2404416386, 1202208193, 3606624580, 1803312290, 901656145, 2704968436, 1352484218, 676242109, 2028726328, 1014363164, 507181582, 253590791, 760772374, 380386187, 1141158562, 570579281, 1711737844, 855868922, 427934461, 1283803384, 641901692, 320950846, 160475423, 481426270, 240713135, 722139406, 361069703, 1083209110, 541604555, 1624813666, 812406833, 2437220500, 1218610250, 609305125, 1827915376, 913957688, 456978844, 228489422, 114244711, 342734134, 171367067, 514101202, 257050601, 771151804, 385575902, 192787951, 578363854, 289181927, 867545782, 433772891, 1301318674, 650659337, 1951978012, 975989006, 487994503, 1463983510, 731991755, 2195975266, 1097987633, 3293962900, 1646981450, 823490725, 2470472176, 1235236088, 617618044, 308809022, 154404511, 463213534, 231606767, 694820302, 347410151, 1042230454, 521115227, 1563345682, 781672841, 2345018524, 1172509262, 586254631, 1758763894, 879381947, 2638145842, 1319072921, 3957218764, 1978609382, 989304691, 2967914074, 1483957037, 4451871112, 2225935556, 1112967778, 556483889, 1669451668, 834725834, 417362917, 1252088752, 626044376, 313022188, 156511094, 78255547, 234766642, 117383321, 352149964, 176074982, 88037491, 264112474, 132056237, 396168712, 198084356, 99042178, 49521089, 148563268, 74281634, 37140817, 111422452, 55711226, 27855613, 83566840, 41783420, 20891710, 10445855, 31337566, 15668783, 47006350, 23503175, 70509526, 35254763, 105764290, 52882145, 158646436, 79323218, 39661609, 118984828, 59492414, 29746207, 89238622, 44619311, 133857934, 66928967, 200786902, 100393451, 301180354, 150590177, 451770532, 225885266, 112942633, 338827900, 169413950, 84706975, 254120926, 127060463, 381181390, 190590695, 571772086, 285886043, 857658130, 428829065, 1286487196, 643243598, 321621799, 964865398, 482432699, 1447298098, 723649049, 2170947148, 1085473574, 542736787, 1628210362, 814105181, 2442315544, 1221157772, 610578886, 305289443, 915868330, 457934165, 1373802496, 686901248, 343450624, 171725312, 85862656, 42931328, 21465664, 10732832, 5366416, 2683208, 1341604, 670802, 335401, 1006204, 503102, 251551, 754654, 377327, 1131982, 565991, 1697974, 848987, 2546962, 1273481, 3820444, 1910222, 955111, 2865334, 1432667, 4298002, 2149001, 6447004, 3223502, 1611751, 4835254, 2417627, 7252882, 3626441, 10879324, 5439662, 2719831, 8159494, 4079747, 12239242, 6119621, 18358864, 9179432, 4589716, 2294858, 1147429, 3442288, 1721144, 860572, 430286, 215143, 645430, 322715, 968146, 484073, 1452220, 726110, 363055, 1089166, 544583, 1633750, 816875, 2450626, 1225313, 3675940, 1837970, 918985, 2756956, 1378478, 689239, 2067718, 1033859, 3101578, 1550789, 4652368, 2326184, 1163092, 581546, 290773, 872320, 436160, 218080, 109040, 54520, 27260, 13630, 6815, 20446, 10223, 30670, 15335, 46006, 23003, 69010, 34505, 103516, 51758, 25879, 77638, 38819, 116458, 58229, 174688, 87344, 43672, 21836, 10918, 5459, 16378, 8189, 24568, 12284, 6142, 3071, 9214, 4607, 13822, 6911, 20734, 10367, 31102, 15551, 46654, 23327, 69982, 34991, 104974, 52487, 157462, 78731, 236194, 118097, 354292, 177146, 88573, 265720, 132860, 66430, 33215, 99646, 49823, 149470, 74735, 224206, 112103, 336310, 168155, 504466, 252233, 756700, 378350, 189175, 567526, 283763, 851290, 425645, 1276936, 638468, 319234, 159617, 478852, 239426, 119713, 359140, 179570, 89785, 269356, 134678, 67339, 202018, 101009, 303028, 151514, 75757, 227272, 113636, 56818, 28409, 85228, 42614, 21307, 63922, 31961, 95884, 47942, 23971, 71914, 35957, 107872, 53936, 26968, 13484, 6742, 3371, 10114, 5057, 15172, 7586, 3793, 11380, 5690, 2845, 8536, 4268, 2134, 1067, 3202, 1601, 4804, 2402, 1201, 3604, 1802, 901, 2704, 1352, 676, 338, 169, 508, 254, 127, 382, 191, 574, 287, 862, 431, 1294, 647, 1942, 971, 2914, 1457, 4372, 2186, 1093, 3280, 1640, 820, 410, 205, 616, 308, 154, 77, 232, 116, 58, 29, 88, 44, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1]

大変だねー。

ずーっと回数回数いってるけど、要素数を数えてるから、回数ではなく要素数じゃないか？っていう疑惑もるけど。。。気にするな！

1も1回なの。一応適応してることにして