ヴィーフェリッヒ素数
http://commutative.world.coocan.jp/blog5/2020/11/post-182.html
関数化と Primes パッケージの isprime() を使うことによりほぼ 15 倍速となった。
元のプログラムと比較して欲しい。
using Primes
function Wiefelich()
z = BigInt(1)
k = BigInt(2)
for p = 3:100000
isprime(big(p)) && k^(p - 1) % (p^2) == 1 && print("$p ")
end
end
@time Wiefelich()
1093 3511 0.069396 seconds (267.17 k allocations: 59.338 MiB, 3.79% gc time)
Julia によるベルヌーイ数の計算
http://commutative.world.coocan.jp/blog5/2020/11/julia-3.html
factorual と binomial を自前の関数を書いているが，既存の関数（binomial だけでよい）を使うのが吉。
肝の Bernoulli 関数を以下のように書き換えた。
符号 s も元のプログラムでは (-1)^j としているがもったいない。
function Bernoulli(k::BigInt)
B = Rational{BigInt}(0)
for m = 0:k
U = Rational{BigInt}(0)
s = 1
for j = 0:m
U += s * binomial(m, j) * j^k
s = -s
end
B += U // (m + 1)
end
B
end
元の Bernoulli(big(601)) は 93 秒ほどかかったが，書き換えたプログラムでは 0.6 秒で終わった。ほぼ 155 倍速になった。
この関数を 2〜1000 までの偶数に対して呼出す。奇数は 0//1 なので省略。
結果を出力すると余分な処理時間がかかる。
function f()
for k = 2:2:1000
#print("B<sub>$k</sub> = ")
b = Bernoulli2(big(k));
#println(b)
end;
end
@time f()
337.636791 seconds (4.57 G allocations: 327.160 GiB, 5.67% gc time)
黄金比の値をニュートン法により小数点以下 n 桁まで求める。
function newton(n=10000)
prec = floor(Int, n / log10(2))
setprecision(prec+10)
x = BigFloat(2)
err = BigFloat(10)
for i in 1:1000
xn = (x^2+1) / (2x-1)
i > 1 && string(xn)[end-5:end-1] == string(x)[end-5:end-1] && return ("converged", i, x)
x = xn
end
("not converged", 1000, x)
end
戻り値は，収束したかどうか，何回目のループか，整数部（1）を含む長精度での黄金比の値（BigBloat）の三組。
@time message, i, x = newton(10000); # 0.006522 seconds (692 allocations: 1.868 MiB)
ループ一回ごとに精度は倍になるので，10000 桁には 14 回で十分。
message, i # ("converged", 15)
出力関数は別に作った。
function printout(x; index=true)
s = string(x)[3:end]
println(" φ = 1.")
for i in 1:50:length(s)
i+49 > length(s) && break
indx = index ? "$((" " * string(i))[end-4:end]):" : ""
println("$indx $(s[i:i+49])")
end
end
printout(x)
φ = 1.
1: 61803398874989484820458683436563811772030917980576
51: 28621354486227052604628189024497072072041893911374
101: 84754088075386891752126633862223536931793180060766
151: 72635443338908659593958290563832266131992829026788
201: 06752087668925017116962070322210432162695486262963
:
9901: 56345156390526577104264594760405569509598408889037
9951: 62079956638801786185591594411172509231327977113803
http://k-ichikawa.blog.enjoy.jp/blog/2012/09/1-4ba1.html で確認したが，合っているようだ。
ユークリッド－マリン数列
http://commutative.world.coocan.jp/blog5/2022/01/post-523.html
元のプログラムでは第 6 項の 50207 までなら 0.035977 秒で計算できる。しかし，第 7 項まで出そうとすると 1912.027524 秒かかる。しかも出てきた答え 547985393 は間違いである。
関数化，isprime()を使用するという小手先では 1887.464951 秒までしか縮まらずしかも出てくる答えは同じく 547985393 で間違いである。
計算のヒントは，元記事に「そもそも（第 7 項は）、340999 とあるけど、これは、2, 3, 7, 43, 139, 50207 をすべて掛け合わせて 1 を足した数である 12603664039 から、それの最大の素因数として求める」と書いてある。
また，数列のことならいつもの「オンライン整数列大辞典」ということで，検索してみると A000946 で，説明も同じ内容である。
Euclid-Mullin sequence: a(1) = 2, a(n+1) is largest prime factor of Product_{k=1..n} a(k) + 1.
元記事には続けて，「12603664039 から 340999 を求めること自体が大変なのだった」と書いているが，どれほど大変なのかプログラムを書いてみた。
using Primes
function Euclid_Mullin_sequence(limit=10)
n = BigInt(2)
println("1: $n")
cum = n
for i in 2:limit
n = factor(Vector, cum + 1)[end]
cum *= n
println("$i: $n")
end
end
関数の引数 limit を 2, 3, 4, ... と増やしていくと，順調に解が求まる。
第 7 項の 340999 までは 0.000345 秒
第 8 項の 2365347734339 までは 0.001223 秒
第 9 項の 4680225641471129 までは 0.322829 秒
第 10 項の 1368845206580129 までは 2.339921 秒
であったが，第 11 項まで求めようとすると相当時間がかかり一向に終わる気配がない。ので，諦めた。
第 11 項は
BigInt(2)*3*7*43*139*50207*340999*2365347734339*4680225641471129*1368845206580129 + 1
= 65127746440715056169703820896799129325103020802826145086639
の素因数分解をしなくてはならない。
高速計算といえば WolframAlpha であるが，流石に計算時間切れであった。
ということで，この問題を解く速いプログラムの定石は，
大変そうかもしれないが別の解法があるなら一応やってみる。
尤も，使える手段（関数）は惜しげなく使う。
かな。
@time Euclid_Mullin_sequence(10)
1: 2
2: 3
3: 7
4: 43
5: 139
6: 50207
7: 340999
8: 2365347734339
9: 4680225641471129
10: 1368845206580129
1.611882 seconds (51.92 M allocations: 1.024 GiB, 24.76% gc time, 0.10% compilation time)
速いときには 2 秒を切ることもある。
@benchmark では最短 1.5 秒程度
