#フィズバズエクストリーム 問題 拙解答の解説 #codeiq

はじめに

CodeIQで出題されていた、通称「フィズ・バズ・エクストリーム」問題（←問題の公開は終了したので左記リンクでは問題閲覧はできません）。
うまい漸化式を見付けて、それを元にコンパクトな再帰で解けたのですが、意外とこの漸化式で解いている人が少なそうなので、（誰か解説書いてくれるだろうと踏んでいたのですけれどこの際だから自分で）それを紹介・解説してみます。
ついでに、そこから Ruby でコードゴルフしたのでそのコードも簡単に解説。

問題（概要）

問題全文は、↓こちらを参照。
CodeIQ 「フィズ・バズ・エクストリーム」 問題 [PDF]

文章を引用する代わりに私の言葉で意訳したのが↓こちら。

与えられた $n$ に対して、$n$ 以下の正の整数のうち、3, 5, 7, 11, 13, 17, 19, 23, 29, 31 の少なくともどれか一つの倍数となるものの総和を求めよ。

• Q1(Normal) $n = 10^5$
• Q2(Hard) $n = 10^9$

包除原理を使って解くところまではすぐに思いついた¹のですが、より発展させて次節のような綺麗な漸化式を導き出しました。

漸化式

$D_m = \{ d_i | d_i > 1 \ {\rm for} \ \forall i: 1 \leq i \leq m, gcd(d_i, d_j) = 1 (i \ne j) \}$ とする。
つまり、$d_1, d_2, \dots , d_m$ を、どの2つも互いに素な整数（>1）とし、その集合を $D_m$ とする。

$F(n, D_m)$ を、「$n$ 以下の正の整数のうち、$D_m$ に属する数（$d_1,d_2,\dots,d_m$）のいずれかの倍数となるものの総和」とする。
（元々の問題は $F(10^5, \{3, 5, 7, 11, 13, 17, 19, 23, 29, 31\})$ および $F(10^9, \{3, 5, 7, 11, 13, 17, 19, 23, 29, 31\})$ を求めるもの、となる）

また $f(n, d_i)$ を、「$n$ 以下の正の整数のうち、$d_i$ の倍数となるものの総和」とする。
（等差級数の公式から、$f(n, d_i) = \frac{d_i n' (n' + 1)}{2}$（ただし $n' := \lfloor n / d_i \rfloor$）である。）

このとき、以下が成り立つ（$m \ge 2$ のとき）：

$$F(n, D_m) = f(n, d_m) + F(n, D_{m-1}) - d_m F(\lfloor n / d_m \rfloor, D_{m-1})$$

証明（概略）

$F(n, D_m)$ は、$f(n, d_i)$ を利用して以下のように書ける（包除原理）：

$$\begin{eqnarray}
F(n, D_m) &=& \ \ f(n, d_1) + f(n, d_2) + \dots + f(n, d_m)\\
& & - (f(n, d_1 d_2) + f(n, d_1 d_3) + \dots + f(n, d_{m-1} d_m))\\
& & + (f(n, d_1 d_2 d_3) + f(n, d_1 d_2 d_4) + \dots + f(n, d_{m-2} d_{m-1} d_m))\\
& & -+ \dots\\
& & + (-1)^{m+1}f(n, d_1 d_2 \dots d_m)
\end{eqnarray}$$

また $d_i, d_j (i \ne j)$ に対して、以下が成立する：

$$f(n, d_i d_j) = d_j f(\lfloor n / d_j \rfloor, d_i)$$

$\because$ $n' := \lfloor n / d_i d_j \rfloor$ とすると、
$$\begin{eqnarray}
n' &=& \lfloor n / d_i d_j \rfloor\\
&=& \lfloor n / d_j / d_i \rfloor\\
&=& \lfloor \lfloor n / d_j \rfloor / d_i \rfloor
\end{eqnarray}$$
（$\lfloor n / d_j / d_i \rfloor \geq \lfloor \lfloor n / d_j \rfloor / d_i \rfloor$ は明らか、$\gt$ と仮定すると $n / d_j \geq \lfloor n / d_j \rfloor + 1$ とならなければならず矛盾）
より
$$\begin{eqnarray}
f(n, d_i d_j) &=& \frac{d_i d_j n' (n' + 1)}{2}\\
&=& d_j \frac{d_i n' (n' + 1)}{2}\\
&=& d_j f(\lfloor n / d_j \rfloor, d_i)
\end{eqnarray}$$

以上より、

$$\begin{eqnarray}
F(n, D_m) &=& \ \{f(n, d_1) + f(n, d_2) + \dots + f(n, d_{m-1})\} + \ f(n, d_m)\\
& & - [\{f(n, d_1 d_2) + f(n, d_1 d_3) + \dots + f(n, d_{m-2} d_{m-1})\} + \{f(n, d_1 d_m) + \dots + f(n, d_{m-1} d_m)\}]\\
& & + [\{f(n, d_1 d_2 d_3) + f(n, d_1 d_2 d_4) + \dots + f(n, d_{m-3} d_{m-2} d_{m-1})\} + \{f(n, d_1 d_2 d_m) + \dots + f(n, d_{m-2} d_{m-1} d_m)\}]\\
& & -+ \dots\\
& & + (-1)^{m+1}f(n, d_1 d_2 \dots d_{m-1} d_m)\\
&=& f(n, d_m)\\
& & + \left(\begin{array}\ \ f(n, d_1) + f(n, d_2) + \dots + f(n, d_{m-1})\\

• (f(n, d_1 d_2) + f(n, d_1 d_3) + \dots + f(n, d_{m-2} d_{m-1}))\\

• (f(n, d_1 d_2 d_3) + f(n, d_1 d_2 d_4) + \dots + f(n, d_{m-3} d_{m-2} d_{m-1}))\\
  -+ \dots\\
• (-1)^mf(n, d_1 d_2 \dots d_{m-1}) \end{array}\right)\\
  & & - \left(\begin{array}\ \ d_m (f(\lfloor n / d_m \rfloor, d_1) + \dots + f(\lfloor n / d_m \rfloor, d_{m-1}))\\

• d_m (f(\lfloor n / d_m \rfloor, d_1 d_2) + \dots + f(\lfloor n / d_m \rfloor, d_{m-2} d_{m-1}))\\
  +- \dots\\

• (-1)^m d_m f(\lfloor n / d_m \rfloor, d_1 d_2 \dots d_{m-1}) \end{array}\right)\\
  &=& f(n, d_m) + F(n, D_{m-1}) - d_m F(\lfloor n / d_m \rfloor, D_{m-1})
  \end{eqnarray}$$

補足（追記 2015/06/29 23:23）

途中の数式はややこしいですが、この漸化式そのものは、
「『$d_1$〜$d_m$ のいずれかの倍数の総和』は、『$d_m$ の倍数の総和』と『$d_1$〜$d_{m-1}$ のいずれかの倍数の総和』から『$d_m$ の倍数でかつ $d_1$〜$d_{m-1}$ のいずれかの倍数 でもある数の総和』を引いたもの」
を意味しています。つまりこれも包除原理を表している、ということ。

提出解（Haskell）

これを利用して、Haskell で解いて提出したのが↓。

answer_q1560.hs

-- answer_q1560.hs
import Data.Time

-- calc
calc :: Integer -> Integer -> Integer
calc n d = 
    let m = n `div` d
    in  d * (div (m * (m + 1)) 2)

-- solve
solve :: Integer -> [Integer] -> Integer
solve 0 _ = 0
solve n (d:[]) = calc n d
solve n (d:ds) = (calc n d) + (solve n ds) - d * (solve (div n d) ds)

main :: IO()
main = do
    st <- getCurrentTime
    print $ solve (10^5) [3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
    -- => 3469796305
    print $ solve (10^9) [3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
    -- => 347147851533059033
    et <- getCurrentTime
    print $ diffUTCTime et st
    -- => 0.00037s

ideone.com での実行結果。

calc 関数が上述の $f$ 関数、solve 関数が $F$ 関数ですね。
まんまなので特に難しいことは無いと思います。

コードゴルフ解（Ruby(79)）

で、コレを元にコードゴルフ²してみました↓。

answer_q1560_golf_79.rb

a=->n,d=1{d<31?a[n,d+=2]-d*52078[d/2]*(a[n/=d,d]+n*~n/2):0}
p a[10**5],a[10**9]

ideone.com での実行結果。

ポイントは2点：

• 漸化式をいかに短く書くか（変数を極力減らす工夫）
• (奇)素数列をいかにごまかしてコード量を減らすか

なお $f$ 関数（Haskell 版の calc 関数）は、n/=d とした上で -d*n*~n/2 としてインラインで表現。d は奇数です（素数とは限らない）。

注目していただきたいのは、52078[d/2] の部分。
Ruby では、Integer#[] は「指定した位置のビット（最下位ビットはインデックス0）が立っていれば1、そうでなければ0」です。
52078は2進表記すると0b1100101101101110。つまり、1, 2, 3, 5, 6, 8, 9, 11, 14, 15 番目の奇数（ただし0番目を1とする、つまり1番目とは3）が素数、というわけ。
で、まずは本来の再帰（a[n,d+=2]の部分）をしてから、d が素数ならば $f(n,d) - d F(\lfloor n/d\rfloor,D_{m-1})$ を計算する（より正確には、その計算結果も考慮する）。
という感じをごちゃごちゃすることで、このコードに到りました。

まとめ

数学系の問題は、数式をこねくり回すと、時として綺麗なカタチに帰着して幸せになれます³。
あと久々にまともに $\rm\TeX$ 数式書いたけれど疲れた。

追記 (2015/06/29 23:23)

同じ漸化式で解いている方の解答をいくつか見付けました。
紹介してみます。

@angel_p_57 さんの解答：Ruby(87)、Perl(117)、C(126)。
特に C でこの短さはすごい。こんな式で3〜31の素数が全て現れるなんてなんかすごい。あと単純に、色々な言語でコードゴルフしてること自体がすごい（私にはその元気がない(^-^;）

id:haruya12 さんの解答・解説：CodeIQ 「フィズ・バズ・エクストリーム」問題。
冒頭に「誰か解説書いてくれるだろう」て書いたけれど投稿時間から言えばこの記事より30分余り早いですね。あとコードゴルフ私のより短い！ 10**5とかを直書きか標準入力かの違いがあるけれど、それを統一しても私のコードの方が2バイト多い！ すごい！ すばらしい！

その他、皆さんの解答が ↓togetter でまとめられています。こちらもご参照あれ。
CodeIQ「フィズ・バズ・エクストリーム」問題 みんなのコード

追記 (2015/07/01 05:55⁴)

haruya12 さんの記事 が追記されて、拙コードゴルフ解の言及と改良案が。

あんちもん2さんの #フィズバズエクストリーム 問題 拙解答の解説 #codeiq のゴルフ解。今は小さい素数表を参照してるわけだけど、これをフェルマーの小定理の逆を使って

a=->n,d=1{d<31?a[n,d+=2]-(2**d%d==2?d*(a[n/=d,d]+n*~n/2):0):0}
p a[105],a[109]

とすると、3バイト長くなるけど素数の上限をもっと伸ばせることに気がついた。
ideone.com で 67 まで行けた。http://ideone.com/9pfafJ

なるほどです。
私のゴルフ解は短くする代わりにパフォーマンスを犠牲（計算しなくても良い項も考慮して無駄な再帰発生）してるのですがそれもクリア（三項演算子?:を使用）してますしね。

逆に、素数の上限をのばすことが目的なら、n の条件を考慮しさらに素数列を降順にしちゃうとさらに幸せだと思うのです。↓ならさらに +2 バイトで済みます。

answer_q1560_golf_84.rb

a=->n,d=33{n*d>0?a[n,d-=2]-(2**d%d==2?d*(a[n/=d,d]+n*~n/2):0):0}
p a[10**5],a[10**9]

d=99（素数の最大値 97）でも ideone で 0.92s と余裕です♪ http://ideone.com/t2wQHT

さらに追記 (2015/07/01 10:35⁵)

↑のを眺めてたら効率無視である程度汎用的でより短いコード↓を思いついてしまった。

answer_q1560_golf_78.rb

a=->n,d=33{d<5?0:a[n,d-=2]-d*4[2**d%d]*(a[n/=d,d]+n*~n/2)}
p a[10**5],a[10**9]

取り敢えず自己記録1バイト縮めた♪
ideone だと d=45（最大素数 43）が限界だけれど。→ http://ideone.com/glewxt

1. | ω・)｡oO （「包除原理」という呼び方そのものはフィードバックもらうまですっかり忘れてた（現役時代『組合せ数学』を専攻していたにもかかわらず）なんて言えない…） ↩

2. | ω・)｡oO （「身体壊すからコードゴルフは控えよう」と心に誓ったはずなのに半年と保たなかったなんて言えない…） ↩

3. | ω・)｡oO （本当はいろんなコード書いて試して「ひょっとしてこゆ漸化式あるんじゃね？」と目星付けてから証明は後付け（しかも提出後）なんて言えない…） ↩

4. | ω・)｡oO （おはようございます…） ↩

5. | ω・)｡oO （はい仕事に戻ります…） ↩