今回は簡単な数理パズルを出題します。プログラムを作って解いてもかまいませんが、なかには筆算 (電卓) で解くことができる問題もあるので、興味のある方は挑戦してみてください。
今のパソコンは高性能なので、次のようにプログラムしても瞬時に答えを求めることができます。
julia> function sum_of_multiples(n, m)
a = zero(n)
for x = m : m : n
a += x
end
a
end
sum_of_multiples (generic function with 1 method)
julia> sum_of_multiples(big(1000000), 3)
166666833333
ところが、数列の和を求める公式を使うと、もっと簡単に答えを求めることができます。
上記公式より n 個の 3 の倍数の和は \(3 + 6 + 9 + \cdots + 3n = \frac{3n(n + 1)}{2}\) となります。したがって、1000000 以下の 3 の倍数の和は 1 から div(1000000, 3) = 333333 までの和を 3 倍することで求めることができます。
3 * 333333 * (333333 + 1) / 2 = 166666833333
これをプログラムすると次のようになります。
julia> function sum_of_multiples1(n, m)
x = div(n, m)
div(m * x * (x + 1), 2)
end
sum_of_multiples1 (generic function with 1 method)
julia> sum_of_multiples1(big(1000000), 3)
166666833333
次のように、一定の差で並んだ数列を「等差数列」といいます。
a を「初項」、d を「公差」といいます。等差数列の一般項は次の式で表すことができます。
初項から \(a_n\) までの和 \(S_n\) は次の式で求めることができます。
初項を 1, 公差 を 1 とすると、1 から n までの和は \(\frac{n(n + 1)}{2}\) となります。
この公式は次のように導出することができます。
このように、右辺を逆順に並べ替えて足し算すると、\(2a + (n - 1)d\) が \(n\) 個並ぶことになります。あとは、これを 2 で割り算すればいいわけです。
10000! であれば、次のようなプログラムでも瞬時に答えを求めることができます。
julia> a = factorial(big(10000)); c = 0; while a % 10 == 0; global c += 1; global a = div(a, 10); end; c 2499
単純に n! を求めて、10 で割れる回数を求めているだけです。ところが、この方法では n が大きくなると極端に遅くなります。多倍長整数の場合、除算や余りを求める処理は乗算よりもはるかに時間がかかります。たとえば、1 桁増やした 100000! の場合、階乗の値は短時間で求めることができても、div(a, 10) の回数が増えることにより実行時間が極端に遅くなるのです。
そこで、他の方法を考えてみましょう。階乗を計算するとき、末尾に 0 が付くのは値を 10 倍したときです。これは数字 10 や 100 を乗算するときだけではありません。次の例を見てください。
1 = 1 1 * 2 = 2 1 * 2 * 3 = 6 1 * 2 * 3 * 4 = 24 1 * 2 * 3 * 4 * 5 = 120 1 * 2 * 3 * 4 * 5 * 6 = 720 1 * 2 * 3 * 4 * 5 * 6 * 7 = 5040 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 = 40320 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 = 362880 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 = 3628800
10 は 2 * 5 に素因数分解することができます。つまり、2 と 5 の組があれば、末尾に 0 がひとつ追加されるわけです。また、0 が複数追加されることもあります。次の例を見てください。
24! = 620448401733239439360000 25! = 15511210043330985984000000
25 は 5 * 5 と素因数分解することができます。このとき、2 * 5 の組が 2 つできるので、末尾に 0 が 2 つ付くわけです。階乗を素因数分解したとき、因数 2 の個数は因数 5 の個数よりも多くなるので、2 と 5 は必ず組にすることができます。つまり、因数 5 の個数が末尾に付く 0 の 個数になるわけです。
階乗の場合、因数の個数を求めるのは簡単です。10000! の場合、10000 / 5 で 5 の倍数の個数 2000 を求めることができます。次に、25 (= 5 * 5) の倍数の個数を 10000 / 25 で求めます。さらに、125 (= 5 * 5 * 5) の倍数の個数を 10000 / 125 で求めます、これを \(10000 \gt 5^m\) が成立する m まで繰り返し、その総和が 5 の因子の個数になります。
10000 / 5 = 2000
10000 / 25 = 400
10000 / 125 = 80
10000 / 625 = 16
10000 / 3125 = 3.2 (小数点切捨て)
----------------------
合計 = 2499
プログラムと実行結果を示します。
julia> function solver(n)
a = 0
m = big(5)
while m <= n
a += div(n, m)
m *= 5
end
a
end
solver (generic function with 1 method)
julia> for x = 1 : 10
println(solver(big(10) ^ x))
end
2
24
249
2499
24999
249998
2499999
24999999
249999998
2499999997
Julia の場合、次のように簡単に求めることができます。
julia> (big(7) ^ 654321) % 10 7
この問題は筆算で簡単に求めることができます。\(7^n \ (n \gt 0)\) の末尾の数字は次のように 7, 9, 3, 1, ... と巡回します。
7^1 = 7 7^2 = 49 7^3 = 343 7^4 = 2401 7^5 = 16807 7^6 = 117649 7^7 = 823543 7^8 = 5764801
ここで、\(7^4\) の末尾は 1 になることに注目してください。末尾が 1 の数字を何回乗算しても、その結果の末尾は 1 になります。\(7^{654321}\) は \((7^4)^{163580} \times 7\) なので、末尾の数字は 7 と求めることができます。
米粒の合計値は \(1 + 2 + 2^2 + 2^3 + \cdots + 2^{80}\) になります。これを素直にプログラムすると次のようになります。
julia> sum([big(2) ^ x for x = 0 : 80]) 2417851639229258349412351
とても大きな数になるので、普通の電卓では計算できません。Windows の電卓を使用するときは関数電卓に切り替えてください。もちろん、数学の公式を使うともっと簡単に求めることができます。
次のように、一定の比で並んだ数列を「等比数列」といいます。
a を「初項」、d を「公比」といいます。等比数列の一般項は次の式で表すことができます。
初項から \(a_n\) までの和 \(S_n\) は次の式で求めることができます。
問題は初項 1 で公比 2 なので、米粒の合計は次のようになります。
julia> big(2) ^ 81 - 1 2417851639229258349412351
この公式は次のように導出することができます。
右辺を引き算すると \(ar\) から \(ar^{n-1}\) の項がなくなって、\(a - ar^n\) だけになります。あとは、\(1 - r\) で割り算すればいいわけです。
7 以上の素数で割り切れない正の整数は、素因子が 2, 3, 5 しかない自然数のことです。これを「ハミング数 (Hamming Numbers)」といいます。ハミング数は素因数分解したとき、\(2^i \times 3^j \times 5^k \ (i, j, k \geq 0)\) の形式になります。たとえば、100 以下のハミング数は次のようになります。
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100
それではプログラムを作りましょう。一番簡単な方法は、1 から n までの整数列を生成して、そこからハミング数を取り出していくことです。これを julia でプログラムすると次のようになります。
julia> function check(n)
while n % 2 == 0; n = div(n, 2); end
while n % 3 == 0; n = div(n, 3); end
while n % 5 == 0; n = div(n, 5); end
n == 1
end
check (generic function with 1 method)
julia> hamming(n) = filter(check, 1 : n)
hamming (generic function with 1 method)
julia> for x = 2 : 8
n = 10 ^ x
m = length(hamming(n))
println("$n, $m")
end
100, 34
1000, 86
10000, 175
100000, 313
1000000, 507
10000000, 768
100000000, 1105
関数 check(n) は n がハミング数かチェックします。これは 2, 3, 5 だけで割り切れるか試しているだけです。プログラムはとても簡単ですが、引数 n の値が大きくなると時間がかかるようになります。
n に比べてハミング数の個数は少ないようなので、式 \(2^i \times 3^j \times 5^k \ (i, j, k \geq 0)\) を使ってハミング数を生成したほうがよさそうです。引数 n に対して i, j, k の上限値は \(\log_2 n, \log_3 n, \log_5 n\) で求めることができます。たとえば、100000000 の場合は次のようになります。
i : 0 - 26 j : 0 - 16 k : 0 - 11
全体で 27 * 17 * 12 = 5508 個しかありません。この中から 100000000 以下の数を選べばいいわけです。プログラムは次のようになります。
リスト : ハミングの問題 (2)
function hamming(n)
zs = []
for x = [big(2) ^ i for i = 0 : floor(Int, log(2, n))]
for y = [big(3) ^ j for j = 0 : floor(Int, log(3, n))]
for z = [big(5) ^ k for k = 0 : floor(Int, log(5, n))]
if x * y * z <= n
push!(zs, x * y * z)
end
end
end
end
sort(zs)
end
2, 3, 5 のべき乗の集合を内包表記で生成し、その要素を掛け合わせて、条件を満たす数値を選択していくだけです。実行結果は次のようになりました。
julia> for x = 8 : 12
n = big(10) ^ x
m = length(hamming(n))
println("$n, $m")
end
100000000, 1105
1000000000, 1530
10000000000, 2053
100000000000, 2683
1000000000000, 3429
この方法だと短時間で答えを求めることができます。
クロージャをサポートしているプログラミング言語では、効率を考慮しないでよければ、クロージャを使って「連結リスト」を実装しすることができます。連結リストの操作関数 cons, car, cdr には次の関係が成り立ちます。
x == car(cons(x, y)) => true y == cdr(cons(x, y)) => true
ここで cons(x, y) で生成したオブジェクトがセルではない場合を考えてみましょう。もし、そのオブジェクトに car を適用すれば cons の第 1 引数 x を返し、cdr を適用すれば第 2 引数を返すことができれば、セルと同じことが実現できます。
そこで、cons はセルではなくクロージャを返すことにしましょう。クロージャは引数 x, y の値を保持することができます。そして、このクロージャは引数に関数を受け取ることにします。あとは、この関数に引数 x, y を渡して評価すれば car と cdr を実現することができます。Julia で cons, car, cdr をプログラムすると次のようになります。
julia> cons(x, y) = z -> z(x, y) cons (generic function with 1 method) julia> car(z) = z((x, y) -> x) car (generic function with 1 method) julia> cdr(z) = z((x, y) -> y) cdr (generic function with 1 method)
関数 cons() はクロージャを返します。このクロージャは引数 z に関数を受け取り、その関数に x, y を渡して評価します。car() は引数 x にクロージャを渡して評価し、第 1 引数 a を返します。これで Lisp / Scheme の car と同じ動作になります。同様に、cdr() は引数 x にクロージャを渡して評価し、第 2 引数 b を返します。これで Lisp / Scheme の cdr と同じ動作になります。
それでは実際に試してみましょう。
julia> a = cons(1, 0) #11 (generic function with 1 method) julia> car(a) 1 julia> cdr(a) 0 julia> b = cons(2, a) #11 (generic function with 1 method) julia> car(b) 2 julia> cdr(b) #11 (generic function with 1 method) julia> car(cdr(b)) 1
このように、クロージャを使って連結リストを作成することができます。ご参考までに簡単なリスト操作関数を下記リストに示します。
#
# list2.jl : クロージャによる連結リストの実装
#
# Copyright (C) 2016-2021 Makoto Hiroi
#
# 空リスト
const nil = :nil
# 基本関数
cons(x, y) = z -> z(x, y)
car(z) = z((x, y) -> x)
cdr(z) = z((x, y) -> y)
# 述語
null(x) = x == nil
consp(x) = typeof(x) <: Function
listp(x) = null(x) || consp(x)
atom(x) = !consp(x)
# リストの表示
function printlist(xs)
print("(")
while consp(xs)
if listp(car(xs))
printlist(car(xs))
else
print(car(xs))
end
if consp(cdr(xs)); print(" "); end
xs = cdr(xs)
end
if !null(xs); print(" . $xs"); end
print(")")
end
# リストの生成
function list(args...)
xs = nil
for i = length(args) : -1 : 1
xs = cons(args[i], xs)
end
xs
end
# リストの連結
function append(xs, ys)
if null(xs)
ys
else
cons(car(xs), append(cdr(xs), ys))
end
end
# マッピング
function mapcar(f, xs)
if null(xs)
nil
else
cons(f(car(xs)), mapcar(f, cdr(xs)))
end
end
# フィルター
function remove(f, xs)
if null(xs)
nil
elseif f(car(xs))
remove(f, cdr(xs))
else
cons(car(xs), remove(f, cdr(xs)))
end
end
# 畳み込み
function fold_left(f, a, xs)
while consp(xs)
a = f(a, car(xs))
xs = cdr(xs)
end
a
end
function fold_right(f, a, xs)
if null(xs)
a
else
f(car(xs), fold_right(f, a, cdr(xs)))
end
end
julia> include("list2.jl")
fold_right (generic function with 1 method)
julia> a = list(1, 2, 3, 4, 5)
#3 (generic function with 1 method)
julia> printlist(a)
(1 2 3 4 5)
julia> b = list(6, 7, 8, 9, 10)
#3 (generic function with 1 method)
julia> printlist(append(a, b))
(1 2 3 4 5 6 7 8 9 10)
julia> printlist(mapcar(x -> x * x, a))
(1 4 9 16 25)
julia> printlist(remove(iseven, a))
(1 3 5)
julia> fold_left(+, 0, a)
15
julia> fold_right(+, 0, a)
15
julia> printlist(fold_right(cons, nil, a))
(1 2 3 4 5)
julia> printlist(fold_left((x, y) -> cons(y, x), nil, a))
(5 4 3 2 1)
Julia の場合、関数をデータ型の一つとして扱うことができるので、関数の返り値として関数を返すことができます。この「関数を返す関数」を使うと、関数の引数が一つでも複数の引数を処理することができます。このような関数を「カリー化関数 (curried function)」といいます。
たとえば、(x, y) -> x + y をカリー化関数にする場合、引数 x を受け取ると「引数 y を受け取って x + y を計算する関数」を返し、その関数に引数 y を渡せば x + y を計算することができます。Julia では、無名関数を使って次のように定義することができます。
julia> foo = x -> y -> x + y #3 (generic function with 1 method)
無名関数を定義する -> は右結合なので、x -> y -> x + y は x -> (y -> x + y) と同じになります。これで引数をひとつ受け取ったら、関数 y -> x + y を返すことができます。もちろん、引数を 2 つ与えれば、それらを加算した結果を返します。つまり、最初の引数を受け取って関数を生成し、その関数を 2 番目の引数に適用する、という動作になります。
それでは実際に試してみましょう。
julia> foo1 = foo(1) #4 (generic function with 1 method) julia> foo1(2) 3 julia> foo(10)(20) 30
引数を一つだけ渡すと「引数 y を受け取って x + y を計算する関数」を返します。返り値を変数 foo1 にセットして、foo1(2) を呼び出せば 1 + 2 = 3 を計算することができます。もちろん、引数を 2 つ渡すと、それらを足し算した値を返します。foo(10) の返り値は関数なので、2 番目の引数を渡すときは foo(10)(20) のようにカッコを続けることに注意してください。
また、局所的な関数を定義して、その関数を返すことでもカリー化関数を実現することができます。次の例を見てください。
julia> function mapcar(f)
function _map(xs)
map(f, xs)
end
_map
end
mapcar (generic function with 1 method)
関数 mapcar() はマッピングを 1 引数の関数に直したものです。mapcar() は関数 f を受け取り、その f を呼び出してリストを操作する関数を返します。これでもマッピングの動作ができるのです。簡単な例を示しましょう。
julia> foo2 = mapcar(x -> x * x)
(::var"#_map#1"{var"#2#3"}) (generic function with 1 method)
julia> foo2([1,2,3,4,5])
5-element Vector{Int64}:
1
4
9
16
25
julia> mapcar(x -> x * x)([1,2,3,4,5])
5-element Vector{Int64}:
1
4
9
16
25
最初の例は mapcar() で生成した関数を変数 foo2 にセットし、それから foo2 を関数呼び出しします。次の例は、mapcar() の返り値を直接関数呼び出ししています。カッコが多くなりますが、2 引数の map() と同じように呼び出すことができます。これでもリストの要素を 2 乗することができます。
Julia では function 文でカリー化関数を定義することはできませんが、関数型言語 (たとえば Haskell や ML (SML/NJ, Ocaml) など) では、簡単にカリー化関数を定義できるようになっています。これらのプログラミング言語では、高階関数はカリー化関数として定義されていて、関数を合成して新しい関数を作ることも簡単にできます。興味のある方は関数型言語にも挑戦してみてください。
「ラムダ計算 (lambda calculus)」は、文字 λ を使って関数を表す「λ 記法」という表記法を用いた抽象的な計算モデルで、1930 年代に A. Church 氏によって考案されました。ラムダ計算は Lisp, Scheme, ML, Haskell など多くの関数型言語の基礎理論として、大きな役割を果たしています。ラムダ計算とか計算モデルというと難しい話のように思われるかもしれません。Lisp / Scheme では無名関数のことを「ラムダ式」といいますが、実をいうとこの無名関数の考え方がラムダ計算の基本なのです。
純粋なラムダ計算の定義はとても単純です。ラムダ計算で扱う式は、次に示す 3 通りしかありません。
関数抽象は関数定義、関数適用は関数呼び出し、変数は関数の仮引数のことと考えてください。つまり、純粋なラムダ計算には関数しかないのです。したがって、ラムダ計算では数を表すのにも関数を使います。これを「チャーチ数 (Church numerals)」と呼びます。ラムダ計算と無名関数はまったく同じではありせんが、今回は難しいことを考えずに Julia の無名関数を使って「チャーチ数」を試してみましょう。
チャーチ数は関数 f と x を受け取り、x に f を適用した回数で数 (自然数) を表します。たとえば、自然数 n は f(f(f(...(f(x)) ...))) のように f を n 回呼び出すことで表します。簡単な例を示しましょう。
julia> zero = f -> x -> x #3 (generic function with 1 method) julia> one = f -> x -> f(x) #7 (generic function with 1 method) julia> two = f -> x -> f(f(x)) #11 (generic function with 1 method) julia> three = f -> x -> f(f(f(x))) #15 (generic function with 1 method)
変数 zero に格納された関数は、引数 f を受け取ったら関数 x -> x を返す関数に、引数 f と x を 2 つ受け取ったら x をそのまま返す関数になります。このとき、x に f を適用していないことに注意してください。これで 0 を表すことができます。
同様に、one は f を 1 回適用しているので 1 を、two は 2 回適用しているので 2 を、three は 3 回適用しているので 3 を表すことができます。そうはいっても、このままではよくわからないので、実際に引数として関数 n -> n + 1 と 0 を渡して実行してみましょう。実行結果は次のようになります。
julia> zero(n -> n + 1)(0) 0 julia> one(n -> n + 1)(0) 1 julia> two(n -> n + 1)(0) 2 julia> three(n -> n + 1)(0) 3
n -> n + 1 は引数 n に 1 を加算するので、もう一つの引数に 0 を渡せば n -> n + 1 を適用した回数、つまりチャーチ数を Julia の数に変換することができます。
次は数 n に 1 を加える succ(f)(n)(x) を定義してみましょう。このとき、n はチャーチ数であることに注意してください。プログラムは次のようになります。
julia> succ = f -> n -> x -> f(n(f)(x)) #27 (generic function with 1 method)
n(f)(x) は数 n を表しているので、それに関数 f を再度適用すれば、n に 1 を加えることができます。簡単な実行例を示します。
julia> succ(n -> n + 1)(zero)(0) 1 julia> succ(n -> n + 1)(one)(0) 2 julia> succ(n -> n + 1)(two)(0) 3 julia> succ(n -> n + 1)(three)(0) 4
このように、succ でチャーチ数 zero, one, two, three に 1 を加算することができます。
次は 2 つのチャーチ数 m, n を足し算する plus(f)(m)(n)(x) を作りましょう。この場合、n(f)(x) で n を表すチャーチ数になるので、この結果に m(f) を適用すれば m + n を実現することができます。プログラムは次のようになります。
julia> plus = f -> m -> n -> x -> m(f)(n(f)(x)) #44 (generic function with 1 method) julia> plus(n -> n + 1)(one)(two)(0) 3 julia> plus(n -> n + 1)(two)(two)(0) 4 julia> plus(n -> n + 1)(two)(three)(0) 5
正常に動作していますね。
次はチャーチ数 m と n を掛け算する mult(f)(m)(n)(x) を定義してみましょう。m * n は n を m 回足し算すればいいので、関数 n(f) を m に渡して m(n(f))(x) とするだけです。プログラムと実行結果を示します。
julia> mult = f -> m -> n -> x -> m(n(f))(x) #58 (generic function with 1 method) julia> mult(n -> n + 1)(two)(three)(0) 6 julia> mult(n -> n + 1)(three)(three)(0) 9 julia> mult(n -> n + 1)(three)(zero)(0) 0 julia> mult(n -> n + 1)(zero)(three)(0) 0
このように、チャーチ数の足し算と掛け算は簡単なのですが、引き算はとても難しくなります。本稿の範囲を超える (M.Hiroi が理解できない) ので、チャーチ数はここまでにしておきましょう。興味のある方は調べてみてください。
変数の個数が式の本数よりも多い方程式を「不定方程式」といいます。一次不定方程式 ax + by = c (c != 0) は、a, b, c が整数で、c が a と b の最大公約数の倍数であるとき、整数解を持つことが証明されています。これを「ベズーの等式」といいます。ベズーの等式は「拡張ユークリッドの互除法」を使って簡単に解くことができます。
拡張といっても原理はユークリッドの互除法とまったく同じです。gcd(a, b) を求めるときに一組の解 (x, y) もいっしょに見つける、というものです。ポイントは数を x = xa * a + xb * b の形で表すところです。次の図を見てください。
x = xa * a + xb * b, y = ya * a + yb * b とおき、
x, y に対してユークリッドの互除法を適用する
(初期値は xa = 1, xb = 0, ya = 0, yb = 1)
x と y の商を q, 余りを z とすると
z = x - q * y = xa * a + xb * b - q * (ya * a + yb * b)
= (xa - q * ya) * a + (xb - q * yb) * b
= za * a + zb * b
となる
次は y と z で同様の処理を繰り返す
具体的には変数の値を
x, xa, xb = y, ya, yb
y, ya, yb = z, za, za
のように更新して処理を繰り返す
y が 0 になったとき、x が a と b の最大公約数で (xa, xb) が一つの解となる
具体的な例を示しましょう。4321 * x + 1234 * y = gcd(4321, 1234) を解いてみます。
(x, xa, xb) (y, ya, yb) q (z, za, zb) ------------------------------------------------------ (4321, 1, 0) (1234, 0, 1) 3 (619, 1, -3) (1234, 0, 1) (619, 1, -3) 1 (615, -1, 4) (619, 1, -3) (615, -1, 4) 1 (4, 2, -7) (615, -1, 4) (4, 2, -7) 153 (3, -307, 1075) (4, 2, -7) (3, -307, 1075) 1 (1, 309, -1082) (3, -307, 1075) (1, 309, -1082) 3 (0, -1234, 4321) (1, 309, -1082) (0, -1234, 4321) 4321 * 309 + 1234 * (-1082) = 1 になる
一般解は次のように求めることができます。
4321 * x + 1234 * y = 1 - 4321 * 309 + 1234 * (-1082) = 1 ----------------------------------- 4321 * (x - 309) + 1234 * (y - (-1082)) = 0 となるから 4321 * (x - 309) = -1234 * (y + 1082) が成り立つ -1234 は x - 309 の約数で 4321 は y + 1082 の約数だから、t を整数とすると、 x - 309 = -1234 * t => x = 309 - 1234 * t y + 1082 = 4321 * t => y = -1082 + 4321 * t となる
それでは、実際に値を代入して (x, y) の組を求めてみましょう。
t x y z = 4321 * x + 1234 * y ----------------------------------------- -4 5245 -18366 1 -3 4011 -14045 1 -2 2777 -9724 1 -1 1543 -5403 1 0 309 -1082 1 1 -925 3239 1 2 -2159 7560 1 3 -3393 11881 1 4 -4627 16202 1
プログラムも簡単です。Julia でプログラムを作ると次のようになります。
リスト : 拡張ユークリッドの互除法
function ext_euclid(a::Int, b::Int)
xs = a, 1, 0
ys = b, 0, 1
while ys[1] != 0
q, z = divrem(xs[1], ys[1])
xs, ys = ys, (z, xs[2] - q * ys[2], xs[3] - q * ys[3])
end
xs
end
x, xa, xb をタプル xs に、y, ya, yb をタプル ys に、z, za, zb をタプルにまとめています。関数 divrem(x, y) は x と y の商と余りを求める関数です。あとはアルゴリズムをそのままプログラムしただけなので、難しいところはないと思います。それでは実際に試してみましょう。
julia> ext_euclid(8, 5) (1, 2, -3) julia> ext_euclid(9, 5) (1, -1, 2) julia> ext_euclid(11, 3) (1, -1, 4) julia> ext_euclid(4321, 1234) (1, 309, -1082) julia> ext_euclid(12357, 100102) (1, 30127, -3719)
拙作のページ Puzzle DE Programming に「水差し問題」がありますが、それを解くのに拡張ユークリッドの互除法を使うことができます。
たとえば、8 リットルの容器 A と 5 リットルの容器 B を使う場合、8 * x + 5 * y = gcd(8, 5) を解くと、8 * 2 - 5 * 3 = 1 という解を見つけることができます。これは A で 2 回汲み B で 3 回捨てると、1 リットルの水が残ることを表します。4 リットルの水を汲みだす場合、両辺を 4 倍して 8 * 8 - 5 * 12 = 4 になります。A で 8 回汲み B で 12 回捨ててもいいのですが、一般解を求めるともっと少ない回数の解を見つけることができます。
8 * x + 5 * y = 4 - 8 * 8 - 5 * 12 = 4 --------------------------------- 8 * (x - 8) + 5 * (y + 12) = 0 x - 8 = -5 * t => x = 8 - 5 * t y + 12 = 8 * t => y = -12 + 8 * t t = 0, x = 8, y = -12 t = 1, x = 3, y = -4 t = 2, x = -2, y = 4
上記のように、B で 4 回汲み A で 2 回捨てる解を見つけることができます。
拡張ユークリッドの互除法は、合同式の徐算で逆数を求めるときにも役に立ちます。\(a \times b \equiv 1 \pmod n\) において a の逆数 b を求める場合、\(a \times x + n \times y = 1\) を求めます。ここで、gcd(a, n) = 1 が条件になることに注意してください。両辺に mod n を施すと \(n \times y\) は 0 になるので、\(a \times x \equiv 1 \pmod n\) が成り立ちます。つまり、x が a の逆数になるわけです。簡単な例を示しましょう。
julia> for a = 1 : 10
println(ext_euclid(a, 11))
end
(1, 1, 0)
(1, -5, 1)
(1, 4, -1)
(1, 3, -1)
(1, -2, 1)
(1, 2, -1)
(1, -3, 2)
(1, -4, 3)
(1, 5, -4)
(1, -1, 1)
n が素数の場合、n とその剰余は互いに素になります。合同式で x の負数 -x は (n - x) になるので、mod 11 における 1 から 10 までの逆数は次のようになります。
x 1/x ------------------------- 0 --- 1 1 1 mod 11 = 1 2 6 12 mod 11 = 1 3 4 12 mod 11 = 1 4 3 12 mod 11 = 1 5 9 45 mod 11 = 1 6 2 12 mod 11 = 1 7 8 56 mod 11 = 1 8 7 56 mod 11 = 1 9 5 45 mod 11 = 1 10 10 100 mod 11 = 1
n が素数でない場合、n とその剰余は互いに素になるとは限りません。次の例を見てください。
julia> for a = 1 : 7
println(ext_euclid(a, 8))
end
(1, 1, 0)
(2, 1, 0)
(1, 3, -1)
(4, 1, 0)
(1, -3, 2)
(2, -1, 1)
(1, -1, 1)
2, 4, 6 と 8 の最大公約数は 1 にならないので、その逆数を求めることはできません。つまり、mod 8 の世界では 2, 4, 6 で割り算することはできないわけです。
また、フェルマーの小定理から逆数を求めることもできます。
\(a \times a^{n-2} \equiv 1 \pmod n\) が成り立つので、逆数は \(a^{n-2}\) になります。それでは実際に試してみましょう。
julia> for x = 1 : 10
println("$x $(x^9) $(x^9 % 11)")
end
1 1 1
2 512 6
3 19683 4
4 262144 3
5 1953125 9
6 10077696 2
7 40353607 8
8 134217728 7
9 387420489 5
10 1000000000 10
拡張ユークリッドの互除法と同じ結果になりましたね。なお、簡単な例ということで累乗を計算してから mod をとりましたが、実用的なプログラムではやらないほうがよいでしょう。mod n の世界で累乗を計算する関数を作ったほうがよいと思います。興味のある方は挑戦してみてください。