Julia, pythonで素因数分解

Julia勉強中です。

https://stackoverflow.com/questions/16996217/prime-factorization-list
に載っている素因数分解をするコードをjuliaで書き直してみました。
途中経過がわかるようにprintlnを追加してあります。

```julia:prime_factor.jl
function factors(n)
    gaps = [1,2,2,4,2,4,2,4,6,2,6]
    length, cycle = 12, 4
    f, fs, nxt = 2, typeof(n)[], 1
    while f * f <= n
        while n % f == 0
            push!(fs,f)
            n0 = n
            n ÷= f # integer division
            println("$n0/$f-->$n")
        end
        f += gaps[nxt]
        nxt += 1
        if nxt == length
            nxt = cycle
        end
    end
    if n > 1 
        push!(fs,n)
    end
    return fs
end
```
これを実行すると、

```julia
@time factors(2^51-1)

2251799813685247/7-->321685687669321
321685687669321/103-->3123162016207
3123162016207/2143-->1457378449
1457378449/11119-->131071
  0.022727 seconds (47.87 k allocations: 2.474 MiB)

5-element Array{Int64,1}:
      7
    103
   2143
  11119
 131071
```
となりました。

pythonの方は、

```python:prime_factor.py
def factors_nojit(n):
    gaps = [1,2,2,4,2,4,2,4,6,2,6]
    length, cycle = 11, 3
    f, fs, nxt = 2, [], 0
    while f * f <= n:
        while n % f == 0:
            fs.append(f)
            n0 = n
            n //= f
            print("{0}/{1}-->{2}".format(n0,f,n))
        f += gaps[nxt]
        nxt += 1
        if nxt == length:
            nxt = cycle
    if n > 1: fs.append(n)
    return fs
```

```python
%time factors_nojit(2**51-1)

2251799813685247/7-->321685687669321
321685687669321/103-->3123162016207
3123162016207/2143-->1457378449
1457378449/11119-->131071
CPU times: user 2.72 ms, sys: 27 µs, total: 2.75 ms
Wall time: 2.37 ms

[7, 103, 2143, 11119, 131071]
```

あれ？pythonの方が速い。
pythonではforやwhileのループは遅くなるので使っちゃいけないと言われてますが、それでもjuliaより速いのはなぜだろう。
ちなみにnumbaで@jitデコレータを使うとかえって遅くなってました。

引数がBigIntになると、、、

```julia:julia

y=2^BigInt(67)-1

147573952589676412927

typeof(y)

BigInt

@time factors(y)

147573952589676412927/193707721-->761838257287
 11.160957 seconds (206.72 M allocations: 3.853 GiB, 32.25% gc time)

2-element Array{BigInt,1}:
    193707721
 761838257287

```

```python:python
y=2**67-1

y

147573952589676412927

%time factors_nojit(y)

147573952589676412927/193707721-->761838257287
CPU times: user 8.7 s, sys: 0 ns, total: 8.7 s
Wall time: 8.7 s

[193707721, 761838257287]
```

やっぱりpythonが速い。


素因数分解は当然専用モジュールがあるので、それらを使ってみます。

```julia:julia
using Primes

@time factor(y)

  0.040054 seconds (392.53 k allocations: 6.728 MiB, 44.22% gc time)

193707721 ⋅ 761838257287
```

```python:python
import sympy

%time sympy.primefactors(y)

CPU times: user 57.2 ms, sys: 11 µs, total: 57.2 ms
Wall time: 56.6 ms

[193707721, 761838257287]
```

今度はjulia < python。

どちらも数十ミリ秒なのでそもそも圧倒的に速くなってますが、専門家が書くとjuliaが上回るのか。。。

juliaは高速化させるのにコツがいる、という印象を受けてます。
それはどの言語も同じなんだろうけど、深く考えなくても速いよっていう触れ込みに過剰に期待しすぎなのかも。

## (追記) print文を消す

@septcolorさんの指摘に基づき、単純にprint文をコメントアウトすると、最適化が効くようでInt64の場合に劇的に速くなりました。

```julia:julia
@btime factors(2^51-1)

  12.398 μs (4 allocations: 384 bytes)

```
だけどBigIntになると逆に遅くなる。print文があったときの7倍ほど。
Int64に最適化したので型変換がループ内で起こっているってことなんでしょう。多分。
メモリの消費量も200Mぐらいから1Gまで増えてしまってます。

```julia:julia
y=2^BigInt(67)-1

 147573952589676412927

@time factors(y)

 79.144071 seconds (1.08 G allocations: 18.470 GiB, 29.19% gc time)
```
引数がInt64のときとBigIntのときとで多重定義してみたんですが、25秒ぐらいまでしか改善しない。
なんでprintあるとBigIntに限ってましなんだろう。。。わからない。


python+numbaだとInt64は同様の速度ですが、BigIntは扱えなくなりました。

```python:python
%timeit factors(2**51-1)

 18.8 µs ± 249 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

y=2**67-1
%time factors(y)

---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
<timed eval> in <module>()

OverflowError: int too big to convert

```

個人的には、急に遅くなるよりもエラーはいてくれたほうがありがたいなー。