JuliaでProject Euler(9)

Problem 9

https://projecteuler.net/problem=9

ピタゴラス数は一般に、

と書ける（mとnは互いに素でどちらかが偶数で、）て、直角三角形の周囲の長さはとなるので、周囲の長さ/2を2つの約数に分けて、さらにその一方を2つの約数に分ける必要があるが、mとm+nが互いに素なのを利用すると、余分な計算をしなくて済む。

function factorize_naive(n)
    fs = Dict()
    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
            fs[d] = e
        end
    end
    if n != 1
        fs[n] = 1
    end
    return fs
end

function f_div(fs1, fs2)
    fs3 = Dict()
    for fac in fs1
        p = fac.first
        e1 = fac.second
        if haskey(fs2, p)
            e = e1 - fs2[p]
            if e != 0
                fs3[p] = e
            end
        else
            fs3[p] = e1
        end
    end
    return fs3
end

function f_value(fs)
    return reduce((x, fac) -> x*fac.first^fac.second, fs, init=1)
end

function divisors(fs)
    L = length(fs)
    if L == 0
        return [Dict()]
    elseif L == 1
        for fac in fs
            p = fac.first
            e = fac.second
            divs = [Dict{Int64,Int64}()]
            append!(divs, [ Dict(p => e1) for e1 in 1:e ])
            return divs
        end
    else
        mid = div(L, 2)
        fss1 = divisors(Dict(Iterators.take(fs, mid)))
        fss2 = divisors(Dict(Iterators.drop(fs, mid)))
        return [ merge(fs1, fs2)
                    for (fs1, fs2) in Iterators.product(fss1, fss2) ]
    end
end

function divisors_coprime(fs)
    w = collect(pairs(fs))
    L = length(w)
    if L == 0
        return [1, 1]
    elseif L == 1
        f = first(fs)
        x = f.first^f.second
        return [1, x]
    else
        v::Array{Int} = []
        mid = div(L, 2)
        for x in divisors_coprime(Dict(w[1:mid]))
            for y in divisors_coprime(Dict(w[mid+1:L]))
                push!(v, x*y)
            end
        end
        return v
    end
end

function divisors3(fs)
    triplets::Array{Tuple{Int,Int,Int}} = []
    num = f_value(fs)
    for fs_ in divisors(fs)
        num1 = f_value(fs_)
        if num1%2 != 0
            continue
        end
        l = div(num, num1)
        for m in divisors_coprime(fs_)
            num2 = div(num1, m)
            n = num2 - m
            if 0 < n && n < m
                push!(triplets, (l, m, n))
            end
        end
    end
    return triplets
end

function Pythagoreans(perimeter)
    fs = factorize_naive(div(perimeter, 2))
    return [ (l*(m*m-n*n), 2*l*m*n, l*(m*m+n*n))
                for (l, m, n) in divisors3(fs) ]
end

function e009(n)
    v = Pythagoreans(n)
    println(length(v), v)
    return reduce(*, v[1])
end

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

素因数分解は、Dict(p=>e)で表している。
ここで、mergeはDict同士をmergeしている。

julia> dic1 = Dict(1=>2, 2=>3)
Dict{Int64,Int64} with 2 entries:
  2 => 3
  1 => 2

julia> dic2 = Dict(3=>3, 4=>2)
Dict{Int64,Int64} with 2 entries:
  4 => 2
  3 => 3

julia> merge(dic1, dic2)
Dict{Int64,Int64} with 4 entries:
  4 => 2
  2 => 3
  3 => 3
  1 => 2

ここではこれがたまたまハマっているが、そんなに使わないと思われる。
Iterators.takeとIterators.dropは素因数分解を半分に使っている。それぞれ、第2引数の個数まで取る、あるいは落としている。

JuliaでProject Euler(8)

Problem 8

https://projecteuler.net/problem=8

function make_table()
    a = [
        "73167176531330624919225119674426574742355349194934",
        "96983520312774506326239578318016984801869478851843",
        "85861560789112949495459501737958331952853208805511",
        "12540698747158523863050715693290963295227443043557",
        "66896648950445244523161731856403098711121722383113",
        "62229893423380308135336276614282806444486645238749",
        "30358907296290491560440772390713810515859307960866",
        "70172427121883998797908792274921901699720888093776",
        "65727333001053367881220235421809751254540594752243",
        "52584907711670556013604839586446706324415722155397",
        "53697817977846174064955149290862569321978468622482",
        "83972241375657056057490261407972968652414535100474",
        "82166370484403199890008895243450658541227588666881",
        "16427171479924442928230863465674813919123162824586",
        "17866458359124566529476545682848912883142607690042",
        "24219022671055626321111109370544217506941658960408",
        "07198403850962455444362981230987879927244284909188",
        "84580156166097919133875499200524063689912560717606",
        "05886116467109405077541002256983155200055935729725",
        "71636269561882670428252483600823257530420752963450"
    ]
    
    return [ parse(Int, c) for s in a for c in s ]
end

function e008(n)
    v = make_table()
    gen_products = (reduce(*, v[i:i+n-1]) for i in 1:length(v)-n+1)
    return maximum(gen_products)
end

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

Pythonのgenerator式みたいなのが使える。
数列の最大を求めるのは、maxじゃなくて、maximum。

JuliaでProject Euler(7)

Problem 7

https://projecteuler.net/problem=7

このあと単純にエラトステネスの篩を使う問題があるので、この問題はナイーブに素数判定するのが正しいと思われる。

Pythonでこう書く。

from itertools import *
import sys

def e007(N):
    primes = [2]
    
    def is_prime_naive(n):
        for p in primes:
            if p*p > n:
                return True
            elif n%p == 0:
                return False
    
    for n in count(3):
        if is_prime_naive(n):
            primes.append(n)
            if len(primes) == N:
                return n

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

Juliaにそのまま翻訳すると、

function e007(N)
    primes = [2]
    
    function is_prime_naive(n)
        for p in primes
            if p*p > n
                return true
            elseif n%p == 0
                return false
            end
        end
    end
    
    for n in Iterators.countfrom(3)
        if is_prime_naive(n)
            push!(primes, n)
            if length(primes) == N
                return n
            end
        end
    end
end

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

Arrayに要素を追加するのは、Arrayのメソッドではなくpush!関数を使う。ここで、!が付与されているのは破壊的な関数を示す。pushは無いようだが、例えばsortには2つバージョンがある。

julia> a = [3, 2]
2-element Array{Int64,1}:
 3
 2

julia> b = sort(a)
2-element Array{Int64,1}:
 2
 3

julia> a
2-element Array{Int64,1}:
 3
 2

julia> sort!(a)
2-element Array{Int64,1}:
 2
 3

julia> a
2-element Array{Int64,1}:
 2
 3

このように、sort!は自分を変えるが、sortは元のArrayは変えず、新たなArrayを作ってそれを返す。

上のコードを実行すると、

$ time julia e007a.jl 10000001
179424691

real    0m37.176s
user    0m37.078s
sys     0m0.656s
$ time pypy e007a.py 10000001
179424691

real    2m28.765s
user    2m28.141s
sys     0m0.484s

Juliaのほうが4倍速い。

JuliaでProject Euler(6)

Problem 6

https://projecteuler.net/problem=6

function e006(n)
    return sum(1:n)^2 - sum(k*k for k in 1:n)
end

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

内包表記はPythonと同様に使える。

べき乗は、**ではなく^を使う。PythonやCはXORを^に充てていたが、それではJuliaはビット演算子をどうしているのだろう。

https://docs.julialang.org/en/v1/manual/mathematical-operations/#Bitwise-Operators-1

julia> 1 ⊻ 1
0

よくわからない記号を使っているが、確かに動く。これよりは関数を使ったほうがよいだろう。

julia> xor(1, 1)
0

JuliaでProject Euler(5)

Problem 5

https://projecteuler.net/problem=5

Pythonで次のように書く。

from fractions import gcd
import sys

def lcm(n, m):
    return n/gcd(n, m)*m

def e005(n):
    return reduce(lcm, range(1, n+1), 1)

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

これをそのままJuliaに翻訳する。

function lcm(n, m)
    return div(n, gcd(n, m))*m
end

function e005(n)
    return reduce(lcm, 1:n, init=1)
end

N = parse(Int, ARGS[1])

println(e005(N))

Juliaにもgcdがある。しかも、モジュールなどの指定がなくても使える。名前が衝突、などと考えない。
reduceも使える。しかし、初期値の指定はデフォルトでは無く、上のように必ずinitという名前付きで指定しなければならない。初期値無しは、reduceの邪道の使い方だと思って、あってもなくても同じときでも必ず初期値を指定していた人間からすると、かなり注意が必要だ。

JuliaでProject Euler(4)

Problem 4

https://projecteuler.net/problem=4

Pythonで次のように書く。

import sys

def digits(n):
    ds = []
    while n > 0:
        n, r = divmod(n, 10)
        ds.append(r)
    return ds

def is_palindromic(n):
    ds = digits(n)
    L = len(ds)
    return all(ds[k] == ds[L-k-1] for k in xrange((L+1)/2))

def e004(E):
    max_a = 10**E-1
    max_n = 0
    for a in xrange(max_a, 10**(E-1)-1, -1):
        for b in xrange(max_a, max(a, (max_n+a-1)/a)-1, -1):
            n = a*b
            if is_palindromic(n):
                max_n = n
                break
    return max_n

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

これをそのままJuliaに翻訳する。

function digits(n)
    a = []
    while n > 0
        r = n%10
        n = div(n, 10)
        push!(a, r)
    end
    return a
end

function is_palindromic(n)
    a = digits(n)
    L = length(a)
    for i in 1:div(L+1, 2)
        if a[i] != a[L-i+1]
            return false
        end
    end
    return true
end

function e004(E)
    local max_n = 0
    local max_a = 10^E-1
    for a in max_a:-1:10^(E-1)
        for b in max_a:-1:max(a, div(max_n+a-1, a))
            n = a*b
            if is_palindromic(n)
                max_n = n
                break
            end
        end
    end
    return max_n
end

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

配列のインデックスは、0ではなく1から始まります。まるで地動説から天動説に逆戻りしたかのようですね。

$ time pypy e004.py 8
9999000000009999

real    0m5.356s
user    0m5.281s
sys     0m0.047s

$ time julia e004.jl 8
9999000000009999

real    0m10.907s
user    0m10.922s
sys     0m0.594s

Juliaが遅いのは、たぶん最初のdigitsの中の、

a = []

で、これだと型が、

Array{Any,1}

で、1次元の配列で要素の型は何でもありとなる。
こういうときは、型を指定するとよい。

    a::Array{Int,1} = []

time julia e004a.jl 8
9999000000009999

real    0m8.420s
user    0m8.516s
sys     0m0.547s

あまり速くならない。
divmodが無いのが効いているのだろうか。

Edit:
PyPyでdivmodをやめたが、速度は変わらなかった。

JuliaでProject Euler(3)

Problem 3

https://projecteuler.net/problem=3

Juliaで次のように書く。

function e003(n)
    local p = 2
    while p*p <= n
        if n%p == 0
            while true
                n = div(n, p)
                if n%p != 0
                    break
                end
            end
            if n == 1
                return p
            else
                return e003(n)
            end
        end
        p += 1
    end
    return n
end

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

Juliaでは、Int同士の割り算を普通にするとFloat64になってしまう。

julia> 4/2
2.0

商を整数にするには、div関数を使う。

julia> div(4, 2)
2

divmodに相当する関数は無いようだ。