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

Julia: プログラムを速くするための定石
https://blog.goo.ne.jp/r-de-r/e/2e21bccd94b2f24cdbc2b54030b13a1a

Julia: プログラムを速くするための定石--その2
https://blog.goo.ne.jp/r-de-r/e/4a815e35e7364368ed963e75a966b595

Julia: プログラムを速くするための定石--その3
https://blog.goo.ne.jp/r-de-r/e/bb84137b4936fb4ab5e76a3cbd7ce9e5

Julia: プログラムを速くするための定石--その4 の続きである
https://blog.goo.ne.jp/r-de-r/e/6b538afd6188153bf5e7db614b67f5e5

で，

ボレル予想
http://commutative.world.coocan.jp/blog5/2021/12/post-497.html

√2 の2進表示中の 0 と 1 の個数の勘定

元のプログラムの改善点は，関数にする，global を使わない，BigFloat の使い方，ビット 0/1 は両方数えなくてもよい，Int(floor(X)) より floor(Int, X) のほうが速い，何より，必要な桁数（precision）を正しく見積もる。

元のプログラムの実行結果は，

number of 0s: 500119, number of 1s: 499882
407.799036 seconds (9.55 M allocations: 4.550 TiB, 23.94% gc time)

書き換えたプログラムは以下の通り。

function f(n=1000000; precision=10000000)
  setprecision(precision)
  X = sqrt(big"2") - 1 # 直接求める
  NO = 1 # considering "1...."
  for i=1:n
    X *= 2
    floor(Int, X) == 1 && (NO += 1; X -= 1) # 整数部が 1, Int(floor(X)) より速い
  end
  NZ = n + 1 - NO # NZ は数えなくてよい
  println("number of 0s: $NZ, number of 1s: $NO")
end

@time f(1000000, precision=1000000)

書き換えの経過

このプログラムは，関数化してもたいして速くならない
number of 0s: 500119, number of 1s: 499882
402.086354 seconds (8.55 M allocations: 4.550 TiB, 23.34% gc time)

1 だけを数えれば，計算時間はほぼ半分になると期待される（実際は半分より長い）
number of 0s: 9500119, number of 1s: 499882
271.447004 seconds (7.05 M allocations: 3.413 TiB, 25.05% gc time)

X -= floor(X) するのは無駄
floor(Int, X) は Int(floor(X)) より速い
number of 0s: 500119, number of 1s: 499882
184.315439 seconds (5.05 M allocations: 2.276 TiB, 29.69% gc time)

X は直接 sqrt(big"2") - 1 でよい
for i=1:1000000 で回しているのは，2進数での桁数。よって，precision でも 1000000 で十分(n より多めにとっておく)
number of 0s: 500119, number of 1s: 499882
 18.305939 seconds (5.01 M allocations: 233.246 GiB, 24.00% gc time)

これにより，
402.086354 / 18.305939 = 21.96480355364453 倍速

円周率 π の2進表示中の 0 と 1 の個数の勘定

function g(n=1000000; precision=10000000)
  setprecision(precision)
  X = big(π) - 3 # 直接求める
  NO = 2 # considering "11...."
  for i=1:n
    X *= 2
    floor(Int, X) == 1 && (NO += 1; X -= 1) # 整数部が 1, Int(floor(X)) より速い
  end
  NZ = n + 2 - NO # NZ は数えなくてよい
  println("number of 0s: $NZ, number of 1s: $NO")
end

@time g(1000000, precision=1100000)

元のプログラムは π をファイルで読んでいるので遅いだろう（ファイルを用意できないので再試できない）。
precision は √2 のときより少し多くしておく。

実行結果

number of 0s: 500279, number of 1s: 499723
 21.048134 seconds (5.00 M allocations: 256.281 GiB, 23.65% gc time)

 

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

Julia: プログラムを速くするための定石
https://blog.goo.ne.jp/r-de-r/e/2e21bccd94b2f24cdbc2b54030b13a1a

Julia: プログラムを速くするための定石--その2
https://blog.goo.ne.jp/r-de-r/e/4a815e35e7364368ed963e75a966b595

Julia: プログラムを速くするための定石--その3 の続きである
https://blog.goo.ne.jp/r-de-r/e/bb84137b4936fb4ab5e76a3cbd7ce9e5

この記事では，適切な手法を選択する必要性を述べる。

リストがソートされている場合，逐次探索ではなく，二分探索すべし。

さて，

素数の末尾の桁の多連星
http://commutative.world.coocan.jp/blog5/2022/01/post-510.html

まず，四連星

元プログラムを以下のように書き換えた。

using Primes

function g(n)
  cnt = zeros(Int, 9)
  a, b, c = 2, 3, 5
  for p in primes(7, n)
     d = p%10
     a == b == c == d && (cnt[a] += 1)
     a, b, c = b, c, d
  end
  for i in [1, 3, 7, 9]
    println("$i$i$i$i: $(cnt[i])")
  end
end

@time g(10^10)

元のプログラム

1111: 592917
3333: 584632
7777: 584595
9999: 593299
177.256503 seconds (1.37 G allocations: 27.177 GiB, 0.76% gc time, 0.00% compilation time)

改良後

1111: 592917
3333: 584632
7777: 584595
9999: 593299
 27.439484 seconds (18.69 k allocations: 5.886 GiB, 0.30% gc time, 0.07% compilation time)

6.5 倍ほど速くなった。

 

次に重連星を探すが，

元記事に「プログラムは、四連星の場合の自明な変更でいいので、ここには示さない。」と書いてあるが，かなり面倒なことになるだろう。
上のプログラムならば，変更は僅かである。

using Primes

function g11(n)
  cnt = zeros(Int, 9)
  a, b, c, d, e, f, g, h, i, j = 2, 3, 5, 7, 1, 3, 7, 9, 3, 9
  prime_list = primes(31, n)
  for p in prime_list
      k = p%10
      if a == b == c == d == e == f == g == h == i == j == k
        cnt[a] += 1
        pos = searchsortedfirst(prime_list, p)
        println(prime_list[pos-10:pos])
      end
      a, b, c, d, e, f, g, h, i, j = b, c, d, e, f, g, h, i, j, k
  end
  for i in [1, 3, 7, 9]
    println("$i$i$i$i$i$i$i$i$i$i$i: $(cnt[i])")
  end
end

@time g11(10^10)

見つかった11連星を表示しないバージョン

11111111111: 0
33333333333: 1
77777777777: 2
99999999999: 0
 29.098075 seconds (201 allocations: 5.885 GiB, 0.51% gc time)

見つかった11連星を表示するバージョン

ほとんど追加時間なく処理を終える。
なお，searchsortedfirst() は二分探索するので速いのだが，もし indexin() などを使うと，大変なことになるので注意。

 [452942827, 452942837, 452942857, 452942927, 452942947, 452942977, 452943047, 452943097, 452943107, 452943137, 452943157]
 [2820087797, 2820087817, 2820087847, 2820087917, 2820087947, 2820087967, 2820088067, 2820088097, 2820088127, 2820088147, 2820088157]
 [6879806623, 6879806663, 6879806723, 6879806753, 6879806783, 6879806803, 6879806833, 6879806873, 6879806893, 6879806903, 6879806933]
 11111111111: 0
 33333333333: 1
 77777777777: 2
 99999999999: 0
  29.096040 seconds (1.34 k allocations: 5.885 GiB, 0.55% gc time)

 

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

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

Julia: プログラムを速くするための定石--その2 の続きである

 

さて，

素数の末尾の桁の三連星

http://commutative.world.coocan.jp/blog5/2021/12/post-504.html

元プログラムは，素数リストをファイルから入力する，プログラムを関数化していないなど，実行速度を上げるためには幾つかの改善が必要である。
実行速度は

157.451095 seconds (1.38 G allocations: 27.355 GiB, 1.26% gc time, 0.00% compilation time)
改良プログラムは

関数化，global を除くことにより，
111: 3615401
333: 3553485
777: 3551409
999: 3614711
 45.741210 seconds (455.21 M allocations: 13.581 GiB, 1.24% gc time, 0.03% compilation time)

ファイルからの入力をやめることにより，
111: 3615401
333: 3553485
777: 3551409
999: 3614711
 29.176872 seconds (125 allocations: 5.885 GiB, 0.02% gc time)

&& の代替処理, a, b, c の効率的な使用, カウント用配列
速度には関係ないが表示を for ループで行うなどで，

111: 3615401
333: 3553485
777: 3551409
999: 3614711
 27.276849 seconds (140 allocations: 5.885 GiB, 0.52% gc time)
最終的には 5.8 倍ほど速くなった。

using Primes
function g(n)
  cnt = zeros(Int, 9)
  a, b = 2, 3
  for p in primes(5, n)
     c = p%10
     a == b == c && (cnt[a] += 1)
     a, b = b, c
  end
  for i in [1,3,7,9]
    println("$i$i$i: $(cnt[i])")
  end
end

@time(g(10^10))

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

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

素数の末尾の桁
http://commutative.world.coocan.jp/blog5/2021/12/post-500.html

10^10 までの素数を検査するのに 133.399550 秒かかるが，そのうちの約半分の 60 秒が素数リストをファイルから入力することで費やされている。

global が使われているが，よほど古いバージョンの Julia が動いていない限りは，不要というより，害悪である。

バージョン 1.4 までは，for ループの内部で外部の変数を参照するときに global の指定が必要だった。

Version 1.4.2 (2020-05-23)
julia> count = 0
0

julia> for i in 1:10
         count += 1
       end
ERROR: UndefVarError: count not defined
Stacktrace:
 [1] top-level scope at ./REPL[2]:2

julia>

julia> for i in 1:10
         global count += 1
       end

julia> count
10

しかし，バージョン 1.5 以降では不要である。

Version 1.5.4 (2021-03-11)

julia> count = 0
0

julia> for i in 1:10
       count += 1
     end

julia> count
10

一般的に global が必要な局面は，ほとんど存在しない。関数の引数にするなど，global を使わないで済ますべきだ。

また，別の記事でも書いたが，プログラムをトップレベルに書くのは良くない。どんなプログラムでも，関数にすべきである。

元のプログラムでは，

upto 10^9
last digit = 1: 12711386
last digit = 3: 12712499
last digit = 7: 12712314
last digit = 9: 12711333
 15.177027 seconds (254.25 M allocations: 4.548 GiB, 2.37% gc time)

upto 10^10
last digit = 1: 113761519
last digit = 3: 113765625
last digit = 7: 113764039
last digit = 9: 113761326
133.399550 seconds (2.28 G allocations: 40.703 GiB, 1.33% gc time)

である。

後に示すようにプログラムを書き換えると，次のように，ほぼ 5 倍速くなる。。

upto 10^9
last digit = 1: 12711386
last digit = 3: 12712499
last digit = 7: 12712314
last digit = 9: 12711333
  2.227844 seconds (115 allocations: 643.507 MiB, 0.86% gc time)

upto 10^10
last digit = 1: 113761519
last digit = 3: 113765625
last digit = 7: 113764039
last digit = 9: 113761326
 28.000456 seconds (119 allocations: 5.885 GiB, 0.04% gc time)

 

なお，10^8 までの結果で，元のプログラムでは
last digit = 7: 1440496
になっているが，
last digit = 7: 1440495
ではないだろうか?

using Primes
function g2(n=10^8)
  cnt = zeros(Int, 9)
  prime_list = primes(2, n)
  for p in prime_list
    cnt[p%10] += 1
  end
  println("last digit = 1: $(cnt[1])")
  println("last digit = 3: $(cnt[3])")
  println("last digit = 7: $(cnt[7])")
  println("last digit = 9: $(cnt[9])")
end