硬貨3種類の Change-Making Problem（お釣り生成問題）を解く計算量 log オーダーのアルゴリズム

ここに記事を書くのは 6年ぶり。

はじめに

Change-Making Problem で硬貨が3種類の場合の効率的なアルゴリズムを思いついたので実装してみました。
（2023年には思いついていたのだけれど記事を書くのをずるずる延ばしていたら 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$ がこれこれの値のとき間違っている」とか「このアルゴリズムはすでに発表されている」などの情報がありましたらぜひ教えてください。

アルゴリズム

(ステップ 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)

$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$ とする）を計算する。

$\lfloor x / r_{2k-1} \rfloor \gt \lfloor y / q_{2k-1} \rfloor$ ならば
$u_1$ を $y$枚、$u_2$ を $z$枚、$1$ を $x$枚使うのが最適として処理を終える。

そうでなければ $i = \lfloor x / r_{2k-1} \rfloor$ とし $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}

この時点で $x \lt r_{2k-1}$ になっている。

(ステップ 4)

$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$枚。

実装

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

    i = x // r1
    i1 = y / q1 if y < i * q1 else i
    y -= i1 * q1
    x -= i1 * r1
    z += i1 * p1

    if 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

    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 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)

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) (y,x,z)
    | i1 < i = 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) (y,x,z)
    | i1 < i = 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