ここに記事を書くのは 6年ぶり。
はじめに
Change-Making Problem で硬貨が3種類の場合の効率的なアルゴリズムを思いついたので実装してみました。
（2023年には思いついて Haskell-jpもくもく会で実装していたのですが記事を書くのをずるずる延ばしていたら 1年以上経ってしまいました）
Change-Making Problem（略称 CMP: お釣り生成問題） とは $m$種類の正整数値の硬貨 $u_1 \gt u_2 \gt ... \gt u_m = 1$ を使って与えられた金額 $n$ を支払う際の最も少ない総枚数はいくつかを求める問題です。
すべての金額を支払えるように $u_m$ を 1 に限定していますが 1 に限定しない問題を考えることもあるようです。
この問題は競技プログラミングで動的計画法を使って解く入門問題としてよく取り上げられますが、そこで紹介される解法の計算量は $O(n m)$ です。
今回思いついたアルゴリズムは硬貨が $u_1 \gt u_2 \gt u_3 = 1$ の3種類の場合限定用ですが計算量は $O(\log u_2)$ です。
ちなみに私がネット検索して見つけた中で最も計算量が少ないのは計算量 $O(n \ log(n) \ log(log(n)))$ のアルゴリズム。 On the Change-Making Problem
（これは硬貨の種類が 3つに限らずいくつでも使える一般的な方法です。
$n \ge u_1 u_2$ の問題は定数オーダーで $n \lt u_1 u_2$ の問題に帰着できるので $O(u_1 u_2 \ log(u_1 u_2) \ log(log(u_1 u_2)))$ とも評価できます）
今回思いついたアルゴリズムと Python と Ruby と Haskel での実装を載せますが、「$u_1,u_2,n$ がこれこれの値のとき間違っている」とか「このアルゴリズムはすでに発表されている」などの情報がありましたらぜひ教えてください。
アルゴリズム
掲載アルゴリズムの改修履歴
- 2024年10月29日: ステップ3の $\lfloor x / r_{2k-1} \rfloor \gt \lfloor y / q_{2k-1} \rfloor$ のときの処理内容の記述に間違いがあったのを修正（実装コードの方は間違っていませんでした）
- 2024年10月29日: $x$ が $0$ のときすぐ処理を終了するように改修（バグだったわけではありません）
(ステップ 1)
$u_1$ と $u_2$ から以下の初期値と漸化式で 4つの整数列 ${a_i}, {r_i}, {p_i}, {q_i}$ を求める。（拡張ユークリッドの互除法的処理）
\begin{align}
初期値 \\
&r_0 = u_1, p_0 = 0, q_0 = 1 （a_0 は定義しない）\\
&r_1 = u_2, p_1 = 1, q_1 = 0 \\
\\
漸化式 \\
&a_i = r_{i-1} を r_i で割った商 \\
&r_{i+1} = r_{i-1} を r_i で割った余り \\
&（つまり r_{i-1} = a_i r_i + r_{i+1} という関係） \\
&p_{i+1} = a_i p_i + p_{i-1} \\
&q_{i+1} = a_i q_i + q_{i-1} \\
\\
終了条件 \\
&r_{i-1} が r_i で割り切れたところで終わりとし、その i を l で表す
\end{align}
(ステップ 2)
変数 $y,x,z$ の初期化
$n$ を $u_1$ で割った商を $y$ に、余りを $x$ に設定し、$z$ を $0$ にする。
（つまり $n = u_1 y + x, z = 0$ という関係）
変数 $y,x,z$ の意味
$y$ は $u_1$ の枚数、$x$ は $1$ の枚数、$z$ は $u_2$ の枚数。
この後の処理で、合計金額 $u_1 y + x + u_2 z = n$ という条件を保ちつつ $y + x + z$ がより小さくなるように $y,x,z$ を更新していく。
以下、$k$ を $1$ から処理していく。
(ステップ 3)
$x = 0$ または $2k-1 \gt l$ ならば
$u_1$ を $y$枚、$u_2$ を $z$枚、$1$ を $x$枚使うのが最適として処理を終える。
$\lfloor x / r_{2k-1} \rfloor$ と $\lfloor y / q_{2k-1} \rfloor$ （ただし $q_{2k-1} = 0$ のときは $\infty$ とする）を計算し
小さい方を $i$ とする。
$\lfloor x / r_{2k-1} \rfloor \gt \lfloor y / q_{2k-1} \rfloor$ のときをケース1、そうでないときをケース2 と呼ぶ。
$y,x,z$ を以下のように更新する。
\begin{align}
&y := y - i \ q_{2k-1} \\
&x := x - i \ r_{2k-1} \\
&z := z + i \ p_{2k-1}
\end{align}
ケース1だったならば
$u_1$ を $y$枚、$u_2$ を $z$枚、$1$ を $x$枚使うのが最適として処理を終える。
この時点で $x \lt r_{2k-1}$ になっている。
(ステップ 4)
$x = 0$ または $2k \gt l$ ならば
$u_2$ を $y$枚、$u_2$ を $z$枚、$1$ を $x$枚使うのが最適として処理を終える。
$x \lt r_{2k}$ ならば $k$ を $1$ 増やして (ステップ 3) の処理に飛ぶ。
$x - r_{2k-1}$ はステップ 3 の最後に「この時点で $x \lt r_{2k-1}$ になっている」と書いたとおりこの値は負になる。
この $x - r_{2k-1}$ を $r_{2k}$ で割った商を $j$、余りを $x'$ とする。
（つまり $x - r_{2k-1} = j r_{2k} + x'$ ここで $j$ は負, $0 \leq x' \lt r_{2k}$ という関係）
$j$ は負で後の処理でわかりづらいので符号を変えて正に更新する。
$j := -j$
$r_{2k-1} - p_{2k-1} + q_{2k-1}$ を $r_{2k} + p_{2k} - q_{2k}$ で割った商を $a$、余りを $b$ とする。
（つまり $r_{2k-1} - p_{2k-1} + q_{2k-1} = a (r_{2k} + p_{2k} - q_{2k}) + b$ という関係）
ただし $a = a_{2k}$ かつ $b = 0$ のときは $a := a - 1$ に修正する。
$j \gt a$ または $y \lt j q_{2k} + q_{2k-1}$ のときは
$u_2$ を $y$枚、$u_2$ を $z$枚、$1$ を $x$枚使うのが最適として処理を終える。
そうでなければ $y,x,z$ を以下のように更新する。
\begin{align}
&y := y - j \ q_{2k} - q_{2k-1} \\
&x := x' （x := x + j \ r_{2k} - r_{2k-1} と同じ） \\
&z := z + j \ p_{2k} + p_{2k-1}
\end{align}
もし $a \lt a_{2k}$ のときは
$u_2$ を $y$枚、$u_2$ を $z$枚、$1$ を $x$枚使うのが最適として処理を終える。
そうでなければ $k$ を $1$ 増やし (ステップ 3) の処理に飛ぶ。
アルゴリズムは以上。
補足
連分数を知っている人向け説明になりますが
\frac {u_1}{u_2} = [a_1; a_2, a_3, \cdots, a_l]
つまり
\frac {u_1}{u_2} =
a_1 + \cfrac{1}{
\,a_2 + \cfrac{1}{
\,a_3 + \cfrac{1}{a_4 + \ddots}
}
}
で
\frac {p_i}{q_i} = [a_1; a_2, a_3, \cdots, a_i]
が $i$-次主近似分数、ステップ 4 に出てきた
\frac {j \ p_{2k} + p_{2k-1}}{j \ q_{2k} + q_{2k-1}} = [a_1; a_2, a_3, \cdots, a_{2k-1}, j]
が中間近似分数になります。
アルゴリズムの適用例
$u_1 = 166, u_2 = 69, n = 2431$ のときを求める。
ステップ 1 の処理
$u_1 = 166, u_2 = 69$ のとき $l$ は $5$ で $r_i, a_i, p_i, q_i$ は以下のようになる。
| $i$ | $0$ | $1$ | $2$ | $3$ | $4$ | $5$ |
|---|---|---|---|---|---|---|
| $r_i$ | $166$ | $69$ | $28$ | $13$ | $2$ | $1$ |
| $a_i$ | $2$ | $2$ | $2$ | $6$ | $2$ | |
| $p_i$ | $0$ | $1$ | $2$ | $5$ | $12$ | $77$ |
| $q_i$ | $1$ | $0$ | $1$ | $2$ | $5$ | $32$ |
ステップ 2 の処理
$n: 2431$ を $u_1: 166$ で割った商 $14$ を $y$、余り $107$ を $x$ とする。
$z$ は $0$ とする。
現在 $(y,x,z) = (14, 107, 0)$
ステップ 3 の処理（k = 1）
$\lfloor x / r_1 \rfloor = \lfloor 107 / 69 \rfloor = 1$ と $\lfloor y / q_1 \rfloor = \lfloor 14 / 0 \rfloor = \infty$ なので
$i = \lfloor x / r_1 \rfloor = 1$ とし、 $y,x,z$ を以下のように更新する。
\begin{align}
&y := y - i q_1 = 14 - 1 * 0 = 14 \\
&x := x - i r_1 = 107 - 1 * 69 = 38 \\
&z := z + i p_1 = 0 + 1 * 1 = 1 \\
\end{align}
現在 $(y,x,z) = (14, 38, 1)$
ステップ 4 の処理（k = 1）
$x - r_1 = 38 - 69 = -31$ を $r_2 = 28$ で割った商 $-2$ を $j$、余り $25$ を $x'$ とする。
$j$ の符号を変えて $j := 2$ とする。
$r_1 - p_1 + q_1 = 69 - 1 - 0 = 68$ を $r_2 + p_2 - q_2 = 28 + 2 - 1 = 29$ で割った商 $2$ を $a$、余り $10$を $b$ とする。
$j \le a$ かつ $y \ge j q_2 + q_1$ なので終了せず $y,x,z$ を以下のように更新する。
\begin{align}
&y := y - j q_2 - q_1 = 14 - 2 * 1 - 0 = 12 \\
&x := x' = 25 \\
&z := z + j p_2 + p_1 = 1 + 2 * 2 + 1 = 6 \\
\end{align}
現在 $(y,x,z) = (12, 25, 6)$
ステップ 3 の処理（k = 2）
$\lfloor x / r_3 \rfloor = \lfloor 25 / 13 \rfloor = 1$ と $\lfloor y / q_3 \rfloor = \lfloor 12 / 2 \rfloor = 6$ なので
$i = \lfloor x / r_3 \rfloor = 1$ とし、$y,x,z$ を以下のように更新する。
\begin{align}
&y := y - i q_3 = 12 - 1 * 2 = 10 \\
&x := x - i r_3 = 25 - 1 * 13 = 12 \\
&z := z + i p_3 = 6 + 1 * 5 = 11 \\
\end{align}
現在 $(y,x,z) = (10, 12, 11)$
ステップ 4 の処理（k = 2）
$x - r_3 = 12 - 13 = -1$ を $r_4 = 2$ で割った商 $-1$ を $j$、余り $1$ を $x'$ とする。
$j$ の符号を変えて $j := 1$ とする。
$r_3 - p_3 + q_3 = 13 - 5 + 2 = 10$ を $r_4 + p_4 - q_4 = 2 + 12 - 5 = 9$ で割った商 $1$ を $a$、余り $1$を $b$ とする。
$j \le a$ かつ $y \ge j q_4 + q_3$ なので終了せず $y,x,z$ を以下のように更新する。
\begin{align}
&y := y - j q_4 - q_3 = 10 - 1 * 5 - 2 = 3 \\
&x := x' = 1 \\
&z := z + j p_4 + p_3 = 11 + 1 * 12 + 5 = 28 \\
\end{align}
現在 $(y,x,z) = (3, 1, 28)$
$a = 1 \lt 6 = a_{4}$ のケースなのでここで処理終了。
最終結果 $(y,x,z) = (3, 1, 28)$
つまり $166$ を $3$枚、$69$ を $28$枚、$1$ を $1$枚 使うのが最適で合計$32$枚。
実装
掲載実装の改修履歴
- 2024年10月29日: x が 0 のときすぐ処理を終了するように改修（バグだったわけではありません）
Python
def cf(u,v):
r1, r2 = u, v
p1, p2, q1, q2 = 0, 1, 1, 0
ds = []
while r2 > 0:
a, r3 = divmod(r1, r2)
ds.append([r2,a,p2,q2])
p1, p2 = p2, a * p2 + p1
q1, q2 = q2, a * q2 + q1
r1, r2 = r2, r3
return [ds[i:i+2] for i in range(0, len(ds), 2)]
def cmp(y, x, ds):
z = 0
for [d1,d2] in ds:
r1, a1, p1, q1 = d1
r2, a2, p2, q2 = d2
if x == 0:
break
i = x // r1
i1 = y / q1 if y < i * q1 else i
y -= i1 * q1
x -= i1 * r1
z += i1 * p1
if x == 0 or i1 < i or not r2:
break
a, b = divmod(r1 - p1 + q1, r2 + p2 - q2)
if a == a2 and b == 0:
a -= 1
if x >= r2:
j2, x2 = divmod(x - r1, r2)
j2 = -j2
y2 = y - j2 * q2 - q1
z2 = z + j2 * p2 + p1
if j2 > a or y2 < 0:
break
y, x, z = y2, x2, z2
if a2 > a:
break
return [y, z, x]
def sol(u1, u2, n):
ds = cf(u1, u2)
y, x = divmod(n, u1)
ns = cmp(y, x, ds)
return [sum(ns), ns]
# u1 = 166, u2 = 69, n = 2431 のとき
# print(sol(166, 69, 2431)) # => [32, [3, 28, 1]]
Ruby
def cf(u,v)
r1, r2 = u, v
p1, p2, q1, q2 = 0, 1, 1, 0
ds = []
while r2 > 0
a, r3 = r1.divmod(r2)
ds.push([r2,a,p2,q2])
p1, p2 = p2, a * p2 + p1
q1, q2 = q2, a * q2 + q1
r1, r2 = r2, r3
end
ds.each_slice(2).to_a
end
def cmp(y, x, ds)
z = 0
ds.each do |d1,d2|
r1, a1, p1, q1 = d1
r2, a2, p2, q2 = d2
break if x == 0
i = x / r1
i1 = if y < i * q1 then y / q1 else i end
y -= i1 * q1
x -= i1 * r1
z += i1 * p1
break if x == 0 || i1 < i
break unless r2
a, b = (r1 - p1 + q1).divmod(r2 + p2 - q2)
a -= 1 if a == a2 && b == 0
if x >= r2
j2, x2 = (x - r1).divmod(r2)
j2 = -j2
y2 = y - j2 * q2 - q1
z2 = z + j2 * p2 + p1
break if j2 > a || y2 < 0
y, x, z = y2, x2, z2
end
break if a2 > a
end
[y, z, x]
end
def sol(u1, u2, n)
ds = cf(u1, u2)
y, x = n.divmod(u1)
ns = cmp(y, x, ds)
[ns.sum, ns]
end
# u1 = 166, u2 = 69, n = 2431 のとき
# p sol(166, 69, 2431) # => [32, [3, 28, 1]]
Haskell
Python版や Ruby版とは少し実装を変えています。（本質的アルゴリズムは同じ）
import Data.List (unfoldr)
import Test.QuickCheck
type YXZ = (Int,Int,Int)
cf :: Int -> Int -> [(Int,Int,Int,Int)]
cf u v = unfoldr f (u,v,0,1,1,0) where
f (_,0,_,_,_,_) = Nothing
f (r1,r2,p1,q1,p2,q2) = Just((r2,a,p2,q2), (r2,r3,p2,q2,p,q)) where
(a, r3) = r1 `divMod` r2
p = a * p2 + p1
q = a * q2 + q1
chunk :: [(Int,Int,Int,Int)] -> [(Int,Int,Int, Int,Int,Int, Int)]
chunk [] = []
chunk [(r1,_,p1,q1)] = [(r1,p1,q1, 0,undefined,undefined, undefined)]
chunk ((r1,_,p1,q1):(r2,a2,p2,q2):ds)
| a' < a2 = [pair]
| otherwise = pair: chunk ds
where
(a, l) = (r1 - p1 + q1) `divMod` (r2 + p2 - q2)
a' = if l == 0 then a-1 else a
pair = (r1,p1,q1, r2,p2,q2, a')
cmp :: Int -> Int -> Int -> (Int, YXZ)
cmp u1 u2 n = (y'+x'+z',(y',x',z')) where
(y, x) = n `divMod` u1
ds = chunk $ cf u1 u2
(y',x',z') = cmp' ds (y,x,0)
cmp' :: [(Int,Int,Int, Int,Int,Int, Int)] -> YXZ -> YXZ
cmp' [] yxz = yxz
cmp' ((r1,p1,q1, r2,p2,q2, a): ds) yxz@(y,x,z)
| x == 0 = yxz
| i1 < i = yxz1
| x1 == 0 = yxz1
| r2 == 0 = yxz1
| x1 < r2 = cmp' ds yxz1
| j2 <= a && y2 >= 0 = cmp' ds yxz2
| otherwise = yxz1
where
i = x `div` r1
i1 = if y < i * q1 then y `div` q1 else i
yxz1@(y1,x1,z1) = (y - i1 * q1, x - i1 * r1, z + i1 * p1) -- u2硬貨を i1*p枚使う
-- x1 < r1 になっている
(j,x2) = (x1 - r1) `divMod` r2
j2 = -j
yxz2@(y2,_,_) = (y1 - j2 * q2 - q1, x2, z1 + j2 * p2 + p1) -- u2硬貨を j2*p2+p1枚使うのがよさそうなら使いたい
-- u1 = 166, u2 = 69, n = 2431 のとき
-- main = print $ cmp 166 69 2431 -- => (32,(3,1,28))
おまけのプロパティベーステスト実装
Haskell での QuickCheck を使ったプロパティベーステスト（素朴な動的計画法での結果と一致するかのテスト）のコードです。
u1: 11万程度以下、u2: 1万程度以下、n: 100万以下という条件でテスト回数 1000万回にしてます。（私のパソコン環境では処理時間 6分くらい）
import Data.List (unfoldr)
import Test.QuickCheck
type YXZ = (Int,Int,Int)
cf :: Int -> Int -> [(Int,Int,Int,Int)]
cf u v = unfoldr f (u,v,0,1,1,0) where
f (_,0,_,_,_,_) = Nothing
f (r1,r2,p1,q1,p2,q2) = Just((r2,a,p2,q2), (r2,r3,p2,q2,p,q)) where
(a, r3) = r1 `divMod` r2
p = a * p2 + p1
q = a * q2 + q1
chunk :: [(Int,Int,Int,Int)] -> [(Int,Int,Int, Int,Int,Int, Int)]
chunk [] = []
chunk [(r1,_,p1,q1)] = [(r1,p1,q1, 0,undefined,undefined, undefined)]
chunk ((r1,_,p1,q1):(r2,a2,p2,q2):ds)
| a' < a2 = [pair]
| otherwise = pair: chunk ds
where
(a, l) = (r1 - p1 + q1) `divMod` (r2 + p2 - q2)
a' = if l == 0 then a-1 else a
pair = (r1,p1,q1, r2,p2,q2, a')
cmp :: Int -> Int -> Int -> (Int, YXZ)
cmp u1 u2 n = (y'+x'+z',(y',x',z')) where
(y, x) = n `divMod` u1
ds = chunk $ cf u1 u2
(y',x',z') = cmp' ds (y,x,0)
cmp' :: [(Int,Int,Int, Int,Int,Int, Int)] -> YXZ -> YXZ
cmp' [] yxz = yxz
cmp' ((r1,p1,q1, r2,p2,q2, a): ds) yxz@(y,x,z)
| x == 0 = yxz
| i1 < i = yxz1
| x1 == 0 = yxz1
| r2 == 0 = yxz1
| x1 < r2 = cmp' ds yxz1
| j2 <= a && y2 >= 0 = cmp' ds yxz2
| otherwise = yxz1
where
i = x `div` r1
i1 = if y < i * q1 then y `div` q1 else i
yxz1@(y1,x1,z1) = (y - i1 * q1, x - i1 * r1, z + i1 * p1) -- u2硬貨を i1*p枚使う
-- x1 < r1 になっている
(j,x2) = (x1 - r1) `divMod` r2
j2 = -j
yxz2@(y2,_,_) = (y1 - j2 * q2 - q1, x2, z1 + j2 * p2 + p1) -- u2硬貨を j2*p2+p1枚使うのがよさそうなら使いたい
cmp_naive :: Int -> Int -> Int -> (Int, YXZ)
cmp_naive u1 u2 n = let (y, x) = n `divMod` u1
in minimum [(i+j+k, (i,k,j)) | i <- [(max 0 (y-u2))..y], let (j,k) = (n - i*u1) `divMod` u2]
u1u2n :: Gen (Int,Int,Int)
u1u2n = do
d1 <- chooseInt (1, 10000)
d2 <- chooseInt (1, 100000)
n <- chooseInt (0, 1000000)
let (u1, u2) = (u2 + d2, 1 + d1)
return (u1, u2, n)
quickCheckN n = quickCheckWith stdArgs { maxSuccess = n }
main = do
quickCheckN 10000000 $ forAll u1u2n $ \(u1,u2,n) -> fst (cmp u1 u2 n) == fst (cmp_naive u1 u2 n)
May Monad be with you