Codes (julia)

˅ 完全楕円積分 Elliptic integral (任意精度)

参考: 第 1 種および第 2 種完全楕円積分の数値計算プログラム 2018.8.3 鈴木実

julia

function th3(q::BigFloat,n::Int64)::BigFloat
    s = zero(BigFloat)
    for i in 1:n
        s += q^(i^2)
    end
    return s*2+1
end

function eps_to_q(eps::BigFloat,n::Int)::BigFloat
    q = zero(BigFloat)
    for i in 1:n
        q = eps_q_for_n(eps,q,i)
    end
    return q
end

function eps_q_for_n(eps::BigFloat,q::BigFloat,n::Int)::BigFloat
    if iseven(n)
        r = eps + 2*q^4*eps
        for i in 4:2:n
            r += (-1)*q^((i-1)^2) + 2*q^(i^2)*eps
        end
    else
        r = eps
        for i in 3:2:n
            r += 2*q^((i-1)^2)*eps - q^(i^2)
        end
    end
    return r
end

#=
q=e + 2q^4e - q^9 + 2q^16e - q^25
0=     - q        - q^9          - q^25
+ e     +2q^4e       + 2q^16e
=#

function sigm(q::BigFloat,ll::Int64)::BigFloat
    s = zero(BigFloat)
    for i in 1:ll
        s += (q^(i<<1))/((1-q^(i<<1))^2)
    end
    return s
end

function yK(k::BigFloat,ll::Int64;is_k_ = false)::BigFloat
    if k^2 >= 0.5 && !is_k_
        return yK_(k,ll; is_k_ = true)
    end

    if !is_k_
        k_ = sqrt(1-k^2)
    else
        k_ = k
    end

    eps = -(1/2)+1/(1+sqrt(k_))

    q = eps_to_q(eps,ll)

    return (pi/2)*(th3(q,ll))^2
end

function yK_(k::BigFloat,ll::Int64;is_k_ = false)::BigFloat
    if k^2 >= 0.5 && !is_k_
        return yK(k,ll; is_k_ = true)
    end

    if !is_k_
        k_ = sqrt(1-k^2)
    else
        k_ = k
    end   

    eps = -(1/2)+1/(1+sqrt(k_))

    q = eps_to_q(eps,ll)

    return (-log(q)/2)*(th3(q,ll))^2
end

function yE(k::BigFloat,ll::Int64;is_k_ = false)::BigFloat
    if k^2 >= 0.5 && !is_k_
        return yE_(sqrt(1-k^2),ll)
    end

    if !is_k_
        k_ = sqrt(1-k^2)
    else
        k_ = k
    end  

    eps = -(1/2)+1/(1+sqrt(k_))

    q = eps_to_q(eps,ll)

    t = th3(q,ll)

    s = sigm(q,1000)


    return pi*( (1/t^2)*(1/BigFloat(6)-4*s) + t^2*((2-k^2)/BigFloat(6)) )
end

function yE_(k::BigFloat,ll::Int64;is_k_ = false)::BigFloat
    if k^2 >= 0.5 && !is_k_
        return yE(k,ll, is_k_ = true)
    end

    if k <= big"4.5e-154" && !is_k_
        return 1
    end

    if !is_k_
        k_ = sqrt(1-k^2)
    else
        k_ = k
    end  

    eps = -(1/2)+1/(1+sqrt(k_))

    q = eps_to_q(eps,ll)

    t = th3(q,ll)

    s = sigm(q,1000)

    lnq = log(q)

    return -t^2*lnq*((1+k^2)/BigFloat(6)) + (1/t^2)*(lnq*(1/BigFloat(6)-4*s)+1)
end

完全楕円積分は MPFR にはない関数の中でも有名じゃね？
実は、SpecialFunctions.jl の実装は piecewise approximation polynomial を使用していると書いています。実際中身を見ると確かにそうなってるので、任意精度にするには別の実装がいります。
あと Elliptic.jl というのもあって、そっちでは引数が ::Float64 にはなっていますが、反復による評価をしているのでコードだけ引き抜いて ::Float64 を ::BigFloat にして任意精度にするのもいいのかもしれません。

参考文献には ln(1/q) と書いているが -ln(q) とした方が精度がちょっと良くなる。
ll は反復を何回するかあと途中に sigm(q,1000) とあるが、この 1000 も反復を何回するかなので、各自で自分の計算にとって適切な値にしてほしい。

yK の最悪精度を調べてみると、各 ll に対して約 ll*4.556+1.09109 桁の精度(10進数で)が得られました

˅ 平方数判定

参考: GNU MP p.39

julia

libgmp = Base.GMP.libgmp

perfect_square(x::BigInt) = Int(ccall((:__gmpz_perfect_square_p, libgmp), Cint, (Ref{BigInt},), x))

あまり知られていないが、GMPの方に平方数判定専用の関数がある。ccallを使って呼び出せば、アセンブラまで最適化された平方数判定関数が手に入る。同じようにn乗根を求めて整数部のみ取り出す関数もGMP側にあるので、立方数以上の判定はそれを使うとかするといい。

˅ ε-算法 (数列の加速法)

参考: お話：数値解析第 6 回収束の加速法 (中編) 長田直樹

julia

epsilon_array1 = zeros(BigFloat,p)
epsilon_array2 = zeros(BigFloat,p)
for i in 1:p
    epsilon_array2[i] = answer_array[i]
end


for i in 1:Int(floor(p / 2.0)-1)
    epsilon_subarray = 1 ./(circshift(epsilon_array2,-1) .- epsilon_array2)

    epsilon_array1 = circshift(epsilon_array1,-1) .+ epsilon_subarray

    pop!(epsilon_array1)
    pop!(epsilon_array2)

    epsilon_subarray = 1 ./(circshift(epsilon_array1,-1) .- epsilon_array1)

    epsilon_array2 = circshift(epsilon_array2,-1) .+ epsilon_subarray

    pop!(epsilon_array2)
    pop!(epsilon_array1)

end

println()
println("epsilon:")
println(stderr,epsilon_array2[end])

println()
println("epsilon:")
println(stderr,epsilon_array1[end])

詳しいことは参考文献を確認してください。
answer_array が収束させたい数列です。
1 から p まで p 項あるとしています。
高速化はしてないです。
epsilon を二つ出していますが、それはどっちに加速された値が出るのか忘れたからです。

˅ 2-強擬素数判定関数

参考: ミラーラビン法による確率的素数判定 (Qiita)

julia

# 2-強擬素数判定関数
function is_strong_pseudoprime(n::BigInt)
    # n-1 = 2^s * d (d は奇数)
    d = n - 1

    s = trailing_zeros(d)
    d = d>>s

    # x = 2^d mod n
    x = powermod(2, d, n)
    if x == 1 || x == n - 1
        return true
    end
    # 反復平方
    for _ in 2:s
        x = (x * x) % n
        if x == n - 1
            return true
        end
    end
    return false
end

trailing_zerosはd*2^sの形にするときによく使うよね

˅ 約数関数

参考: 約数関数 (Wikipedia)

julia

"""
    sigma(n::Integer) -> Integer

与えられた自然数 `n` のすべての正の約数の和 sigma(n) を返す。
"""
function sigma(n::Integer)::BigInt
    n == 1 && return BigInt(1)
    f = factor(n)
    s = BigInt(1)
    for (p, e) in f
        s *= div((BigInt(p)^(e + 1) - 1) , (BigInt(p) - 1))
    end
    return s
end

"""
    sigma(n::Integer, x::Integer) -> Integer

与えられた自然数 `n` のすべての正の約数のx乗和 sigma(n,x) を返す。
"""
function sigma(n::Integer, x::Integer)::BigInt
    n == 1 && return BigInt(1)
    f = factor(n)
    s = BigInt(1)
    if x == 1
        for (p, e) in f
            s *= (e+1)
        end
        return s
    end
    for (p, e) in f
        s *= div((BigInt(p)^((e + 1)*x) - 1) , (BigInt(p)^x - 1))
    end
    return s
end

素因数分解して積の形でやると早いよね

˅ Pochhammer記号

参考: ポッホハマー記号

julia

using SpecialFunctions
libgmp = Base.GMP.libgmp
function divexact(x::BigInt,y::BigInt)::BigInt
    z = BigInt()
    ccall((:__gmpz_divexact, libgmp), Cvoid, (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}), z, x, y)
    z
end

function pochhammer_bigint(x::BigInt, n::BigInt)::BigInt
    if n <= 0
        throw(DomainError(n, "pochhammer(x, n) is undefined for n ≤ 0"))
    end
    result = one(BigInt)
    for i in 0:BigInt(n - 1)
        result *= x + i
    end
    return result
end

function pochhammer_int(x::Int64, n::Int64)::BigInt
    if n <= 0
        throw(DomainError(n, "pochhammer(x, n) is undefined for n ≤ 0"))
    end
    result = one(BigInt)
    for i in 0:(n - 1)
        result *= x + i
    end
    return result
end

function pochhammer_gmp(x::BigInt, n::BigInt)::BigInt
    if n <= 0
        throw(DomainError(n, "pochhammer(x, n) is undefined for n ≤ 0"))
    end
    return divexact(factorial(x+n-1),factorial(x-1))
end

function pochhammer_mpfr(x::BigFloat, n::BigFloat)::BigFloat
    if (isinteger(x) && x<=0) || (isinteger(x+n) && x+n<=0)
        throw(DomainError(n, "pochhammer(x, n) is undefined for x or x+n is a non-positive integer."))
    end
    return div(SpecialFunctions._gamma(x+n),SpecialFunctions._gamma(x))
end

x に比べて n がまぁまぁ小さいときは pochhammer_bigint を使うべき。
result_200.png は n が 200 に固定されているときの pochhammer_bigint と pochhammer_gmp の実行時間比較。
(他のpngも同様にnがpngの名前に書いている値に固定されている)
x が小さいときは赤(gmp)が早いが n に比べて x がでかくなってくると青(bigint)のほうが早い
そりゃそうすぎる。 (まだ書きかけ)

˅ 1~100までの素数を出すコードゴルフ (外部パッケージなし)

julia

r=2:99;r[r.∉[r*r']]

まずjuliaでは ; が文の区切りを表せるので、これは r=2:99 と r[r.∉[r*r']] の2文に分かれる
r=2:99 は2から99までの整数の行ベクトル [1,2,3,...,99] (正確には違うが)
r[r.∉[r*r']] だが、まず ' というのは adjoint の短縮記法で、 adjoint という関数は
随伴行列(や、エルミート転置共役転置とか言われる)をとる関数。
ここでは r に対して r' としているので、 r' が列ベクトルになる (と思ってくれていい)
よって r*r' は 2~99 までのかけ算表 (98次正方行列) ができる
[r*r'] で「要素に98次正方行列を一つだけ持つ要素数が1の配列」になる。
A∉B はAがBに含まれるなら true 、含まれないなら false を返す関数で、
r.∉[r*r'] とすることで各 r の要素に対する挙動へと変えている。
つまり、 r の各要素 i に対して、 i∉[r*r'] を計算し、その結果をまとめて配列で返している。
r.∉r*r' とはできない。これは各 1≤i,j≤98 に対して r[i] != (r*r')[i,j] の真偽を98次正方行列に入れてる気がする。
( r[3] で行ベクトル r の3列目、 (r*r')[4,5] で行列 r*r' の4行5列目を指す。)
つまりは r.∉[r*r'] は r の各要素について2~98のかけ算表に入っていないかどうかを配列で返す。
r[真か偽かの配列] で真の部分だけ残した配列が得られる(便利)ので、これで素数が得られる。