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を用いる。
なお、初期化で型を指定していないが、除算に時間がかかっているためあまり影響がない。

JuliaでProject Euler(11)

Problem 11

https://projecteuler.net/problem=11

generatorを使うと吉。

function make_table()
    a = [
        "08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08",
        "49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00",
        "81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65",
        "52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91",
        "22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80",
        "24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50",
        "32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70",
        "67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21",
        "24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72",
        "21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95",
        "78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92",
        "16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57",
        "86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58",
        "19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40",
        "04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66",
        "88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69",
        "04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36",
        "20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16",
        "20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54",
        "01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48",
    ]
    
    return [ [ parse(Int, s) for s in split(line) ] for line in a ]
end

function e011(n::Int64)
    g() = Channel() do c
        table = make_table()
        H, W = length(table), length(table[1])
        dirs = [(1, 0), (0, 1), (1, 1), (1, -1)]
        for (dx, dy) in dirs
            for (x0, y0) in Iterators.product(1:W, 1:H)
                x3, y3 = x0+dx*(n-1), y0+dy*(n-1)
                if 1 <= x3 <= W && 1 <= y3 <= H
                    push!(c, [ table[y0+dy*k][x0+dx*k] for k in 0:n-1 ])
                end
            end
        end
    end
    
    return maximum(reduce(*, v) for v in g())
end

N = parse(Int64, ARGS[1])
println(e011(N))

productはPythonと同様に使える。

JuliaでProject Euler(10)

Problem 10

https://projecteuler.net/problem=10

エラトステネスの篩を行うだけ。

function make_prime_table(N::Int64)
    a = [ true for n in 1:N ]
    for p in 2:N
        if p*p > N
            break
        elseif !a[p]
            continue
        end
        for n in p*p:p:N
            a[n] = false
        end
    end
    
    return [ n for n in 2:N if a[n] ]
end

function e010(n::Int64)
    return sum(make_prime_table(n-1))
end

N = parse(Int64, ARGS[1])
println(e010(N))

書き忘れたが、a:b:cは、bがステップ幅。