Julia: プログラムを速くするための定石

弱いゴールドバッハ予想の計算
http://commutative.world.coocan.jp/blog5/cat122/

Julia プログラムを速くするために何をすべきかの例を示すための素材として取り上げる。
元のプログラムは，上のリンクをクリック。

冒頭 12 行は，素数リストをファイルから詠み込んでいる。
しかしこれは，Using Primes; list = primes(2, prime(1000)) とすれば簡単である（時間はかからない）。使えるものは何でも使おう。

さて，元のプログラムはトップレベルに書かれているが，Julia のプログラムを速くするための定石は，「どんなに短くて，一回しか使わないプログラムでも関数にする」ということである。
関数にするだけで，元のプログラムが 40 秒ほどかかっているのが 0.4 秒と，100 倍速くなる。なお，結果の出力はコメントアウトし，条件を満たす個数のみを勘定するようにしている。

次いで，結果のチェック時に a <= b && b <= c をしているが，素数リストの添字 a, b, c のループの開始値を変えれば，ループ回数も減らせるし，このチェック自体不要になる。
この最適化によって，実行時間は 0.075 秒ほどになる。

k の計算に while を使っているが，これと同じことを 1 行で行う。
この最適化は，計算時間にはあまり効果を持たない。

最後に，3 重の for ループで，無駄なループをしないようにする（それまでの和がすでに r を超えれば，それ以降のループは不要）。
この最適化により，実行時間は 0.025 秒ほどになる。

結果，元のプログラムより 1500 倍ほど速くなった。39.674786 / 0.025536 ≒ 1553 

なお，入出力は通常の計算より 3 桁ほどの処理時間を要するので，結果を出力し，出力時間も含めた合計処理時間にはほとんど変化はないというような場合もある。

参考までに，以下に，書き換えたプログラムを示しておく。

function f(n=1000)
  list = primes(2, prime(n));
  # max = i-1
  num = 0
  for r = 7:2:n+1
    k = sum(list .<= r) + 1
    for a = 1:k
      list[a] >= r && break
      for b = a:k
        list[a] + list[b] >= r && break
        for c = b:k
          t = list[a] + list[b] + list[c]
          t > r  && break
          if t == r
            num += 1
            # print("$r = ")
            # print(list[a])
            # print("+")
            # print(list[b])
            # print("+")
            # println(list[c])
          end
        end
      end
    end
  end
  println("num = $num")
end

@time f()

実行時間の変化

関数にする

num = 198418
  0.407057 seconds (1.19 k allocations: 29.109 KiB)
  39.674786 / 0.407057 # 97.4673964579899 倍速くなった

b, c の初期値を a, b にする

num = 198418
  0.074374 seconds (1.19 k allocations: 29.109 KiB)

k の計算に while を使わない

num = 198418
  0.065854 seconds (2.68 k allocations: 2.171 MiB)

a, b, c の for loop で無駄なループをしない（break する）

num = 198418
  0.025536 seconds (2.68 k allocations: 2.172 MiB)
  39.674786 / 0.025536 # 1553.6805294486214 倍速くなった

Julia: 長精度計算するときの注意事項

Juliaによる累乗根の計算
http://commutative.world.coocan.jp/blog5/2021/04/julia-9.html

setprecision(100)
x = BigInt(2)
T = BigFloat(0.2)
x = x^T
println("2^0.2 = $x")
x = x^5
println("(2^0.2)^5 = $x")

の結果が

2^0.2 = 1.1486983549970350156384116963047
(2^0.2)^5 = 2.0000000000000000769547959311696だったので，x と T の定義を以下のように

x = BigInt(2)
T = one(BigInt)
T /= 5して

2^(1/5) = 1.1486983549970350067986269467779
(2^(1/5))^5 = 2.0が得られたとして，

> JuliaのBigFloatの計算において、多倍長浮動小数点数の指数は、可能な限り、BigInt/BigIntの分数として与えるべきであり、下手にそれを割り算して多倍長浮動小数点数にしてしまわない方がいいということである。

> なぜそうなのかと問うても仕方ない。それが(現バージョンの)Juliaの仕様ということである。

対応法としては間違いがないが，「BigInt/BigIntの分数として与えるべき」というのは不正確である。

要するに，長精度計算をする場合でも，BigFloat(0.2) は「64bit 浮動小数点数においては不正確な 0.2 を長精度で表しただけということだ。これを避けるには「BigInt/BigIntの分数として与える」のも正しいが，いつも整数の分数として与えることができるとは限らないこともある。そのような場合には， BigFloat の引数を文字列にするのである。上の例だと BigFloat("0.2") そうすれば，長精度の 0.2 が得られるのだ。big"0.2" でもよい（カッコを付けてはいけないので注意）。

BigFloat(0.2)    # 0.20000000000000001110223024625157
BigFloat("0.2")　# 0.20000000000000000000000000000004
big"0.2"         # 0.20000000000000000000000000000004