JuliaでProject Euler(18)

Problem 20

https://projecteuler.net/problem=20

Problem16とだいたい同じ。

Pythonで書くと、

from itertools import *
import sys

def multiply(v, n):
    v2 = []
    carry = 0
    for d in v:
        carry, r = divmod(d*n+carry, 10)
        v2.append(r)
    while carry > 0:
        carry, r = divmod(carry, 10)
        v2.append(r)
    return v2

def e020(N):
    a = [1]
    for n in range(1, N+1):
        a = multiply(a, n)
    return sum(a)

N = int(sys.argv[1])
print e020(N)

これをJuliaに翻訳して、

function multiply(v, n)
    v2::Array{Int} = []
    carry = 0
    for d in v
        m = d*n+carry
        r = m%10
        carry = div(m, 10)
        push!(v2, r)
    end
    while carry > 0
        r = carry%10
        carry = div(carry, 10)
        push!(v2, r)
    end
    return v2
end

function e020(N)
    a = [1]
    for n in 1:N
        a = multiply(a, n)
    end
    return sum(a)
end

N = parse(Int, ARGS[1])
println(e020(N))

やはり、2行目で型指定していないと遅い。
PyPyで、

$ time pypy e020.py 30000
511470

real    0m15.933s
user    0m13.875s
sys     0m2.016s

Juliaで型指定しないと、

$ time julia e020.jl 30000
511470

real    2m6.301s
user    2m3.984s
sys     0m0.594s

上のように型指定すると、

$ time julia e020a.jl 30000
511470

real    0m10.654s
user    0m10.734s
sys     0m0.500s

PyPyより速くなった。

JuliaでProject Euler(17)

Problem 17

https://projecteuler.net/problem=17

何度も解いているが、一度も最初で正しい答えを出した記憶がない。
できるだけ再帰を使う。

function e017(N)
    function num_length(n)
        if haskey(num_letters, n) && n < 100
            return length(num_letters[n])
        elseif n < 100
            d = n%10
            return num_length(n-d) + num_length(d)
        elseif n < 1000
            d = div(n, 100)
            r = n%100
            if r == 0
                return num_length(d) + length(num_letters[100])
            else
                return num_length(d*100) + 3 + num_length(r)
            end
        else
            d = div(n, 1000)
            r = n%1000
            if r == 0
                return num_length(d) + length(num_letters[1000])
            else
                return num_length(d*1000) + 3 + num_length(r)
            end
        end
    end
    
    num_letters = Dict(
        0 => "", 1 => "one", 2 => "two", 3 => "three", 4 => "four", 5 => "five",
        ...
    )
    
    return sum(num_length(n) for n in 1:N)
end

N = parse(Int, ARGS[1])
println(e017(N))

JuliaでProject Euler(16)

Problem 16

https://projecteuler.net/problem=16

本当は速い方法(メルセンヌ数を10進表示する)があるが、ここでは単に2倍ずつしていく。

Pythonでは、

import sys

def twice(v):
    v2 = []
    carry = 0
    for d in v:
        carry, r = divmod(d*2+carry, 10)
        v2.append(r)
    if carry > 0:
        v2.append(carry)
    return v2

def e016(N):
    a = reduce(lambda x, y: twice(x), xrange(N), [1])
    return sum(a)

N = int(sys.argv[1])
print e016(N)

Juliaにそのまま翻訳。

function twice(v)
    v2 = []
    carry = 0
    for d in v
        n = d*2+carry
        r = n%10
        carry = div(n, 10)
        push!(v2, r)
    end
    if carry > 0
        push!(v2, carry)
    end
    return v2
end

function e016(N)
    a = reduce((x, y) -> twice(x), 1:N, init=[1])
    return sum(a)
end

N = parse(Int, ARGS[1])
println(e016(N))

実行してみると、

% time julia e016.jl 100000
135178

real    1m27.103s
user    1m27.047s
sys     0m0.531s

% time pypy e016.py 100000
135178

real    0m9.283s
user    0m8.547s
sys     0m0.719s

2行目はやはり、型を指定してやらないといけない。

    v2::Array{Int} = []

こうすると、PyPyと同様になる。

time julia e016a.jl 100000
135178

real    0m9.041s
user    0m9.203s
sys     0m0.453s

PyPyでは型が推定できているのだろう。

JuliaでProject Euler(15)

Problem 15

https://projecteuler.net/problem=15

DPで解く。

function e015(N)
    a = zeros(Int, N+1, N+1)
    a[1,1] = 1
    for x in 2:N+1
        a[1,x] = 1
    end
    for y in 2:N+1
        a[y,1] = 1
    end
    for y in 2:N+1
        for x in 2:N+1
            a[y,x] = a[y,x-1] + a[y-1,x]
        end
    end
    
    return a[N+1,N+1]
end

N = parse(Int, ARGS[1])
println(e015(N))

Juliaでは、配列の配列だけでなく多次元配列が使える。
アクセス方法は上の通り。

JuliaでProject Euler(14)

Problem 14

https://projecteuler.net/problem=14

Collatzはこう書ける。

function Collatz_length(n)
    len = 1
    while n != 1
        if n%2 == 0
            n >>= 1
        else
            n = n*3+1
        end
        len += 1
    end
    return len
end

function e014(N::Int)
    max_length, max_n = maximum((Collatz_length(n), n) for n in 1:N-1)
    return max_n
end

N = parse(Int, ARGS[1])
println(e014(N))

PyPyでもそうだが、なぜかメモ化すると遅くなる。なぜだろう。
再帰を使わなければ、メモ化が効く。

function e014(N::Int)
    memo = [ 0 for _ in 1:N-1 ]
    memo[1] = 1
    for n in 2:N-1
        counter = 0
        m = n
        while m >= N || memo[m] == 0
            if m%2 == 0
                m = m >> 1
            else
                m = m*3 + 1
            end
            counter += 1
        end
        memo[n] = counter + memo[m]
    end
    max_length, max_n = maximum((memo[n], n) for n in 1:N-1)
    return max_n
end

N = parse(Int, ARGS[1])
println(e014(N))

10^8で上のコードより5倍くらい速くなった。

JuliaでProject Euler(13)

Problem 13

https://projecteuler.net/problem=13

多倍長整数を自作してもいいが、ここではBigIntを用いる。

function make_big_numbers()
    a = [
        "37107287533902102798797998220837590246510135740250",
        "46376937677490009712648124896970078050417018260538",
        ...
        "53503534226472524250874054075591789781264330331690",
    ]
    return [ parse(BigInt, s) for s in a ]
end

function e013()
    v = make_big_numbers()
    return string(sum(v))[1:10]
end

println(e013())

Juliaでは、勝手に多倍長整数になったりしない。

julia> 10^9 * 10^10
-8446744073709551616

IntのBigIntへの変換は次のように行う。

julia> convert(BigInt, 10^9) * 10^10
10000000000000000000

JuliaでProject Euler(12)

Problem 12

https://projecteuler.net/problem=12

本当は速い方法があるが、ここでは愚直に順番に約数の個数を調べていく。

function factorize_naive(n)
    fs = []
    for d in Iterators.countfrom(2)
        if d*d > n
            break
        elseif n%d == 0
            e = 1
            n = div(n, d)
            while n%d == 0
                e += 1
                n = div(n, d)
            end
            push!(fs, Pair(d, e))
        end
    end
    if n != 1
        push!(fs, Pair(n, 1))
    end
    return fs
end

function num_divisors(fs)
    return reduce((x, f) -> x*(f.second+1), fs, init=1)
end

function T(n)
    return div(n*(n+1), 2)
end

function e012(N)
    for n in Iterators.countfrom(1)
        fs = factorize_naive(T(n))
        num = num_divisors(fs)
        if num > N
            return T(n)
        end
    end
end

N = parse(Int, ARGS[1])
println(e012(N))

素因数分解の結果をArray{Pair{Int64,Int64}}で表している。Pairの要素へのアクセスは、first, secondを用いる。
なお、初期化で型を指定していないが、除算に時間がかかっているためあまり影響がない。