2018-03-21

第５問「10/16数」問題解答

「10/16数」問題の解法です。
問題文はこちら。
他の問題の解法は、第１問、第２問、第３問、第４問。

自然数 n を 16 進数のように読んで 10 進数で表したものを F(n) とするときの、
F(n) が n の整数倍となる n を探す問題です。

F(n) を単純にコーディングするだけなら、さほど難しくありません。
n を一つ一つ調べていくコードは次のようになります。

def f(n):
    r, t = 0, 1
    while n > 0:
        r += t * (n % 10)
        n, t = n / 10, t * 16
    return r

def g(m):
    c, n = 0, 1
    while c < m:
        if f(n) % n == 0:
            c += 1
            if c == m: return n
        n += 1

print g(100)

ただこれだと、100 番目の 10/16 数を求めるには長い時間がかかってしまいます。
10/16 数というのがなかなかにレアな存在なためです。

うまく候補を絞り込んで探索する必要があります。
そこでまず、桁数を d とし、n の各桁を a₁a₂・・a_d と表します。
さらに F(n)÷n の値を r とします。

例えば d = 4，r = 3 のとき、a₁, a₂, a₃, a₄ が満たすべき等式は何でしょうか。
それは次の式です。

$3(1000a_1 + 100a_2 + 10a_3+ a_4) = 4096a_1 + 256a_2 + 16a_3 + a_4$

整理するとこうなります。

$1096a_1 - 64a_2 - 14a_3 - 2a_4 = 0$

a₁ の係数が大きいことがポイントです。
例えば a₁ = 1 としてみます。第 1 項の値は 1096 です。
a₂, a₃, a₄ は 0～9 の自然数ですから、左辺の第 2 項以降は、どれだけ小さくとも－64・9－14・9－2・9 = －720 にしかなりません。第 1 項の値を打ち消すには足りません。
したがって、具体的な a₁, a₂, a₃, a₄ の値を実際に試してみるまでもなく、この等式は決して解を持たないと分かります。

このように、上の桁（a₁）から数字を仮決めして、以下の項でどうやっても等式を成立（辺の値をゼロにする）させられないと分かれば、以降の探索を中断する、という枝刈りのアプローチです。

コード例は以下です。ちょっとややこしいです。
始めに桁数 d と、F(n)÷n の値 r を固定します。
次に深さ優先探索で、上の桁から順に数字を仮決めしながら探索していきます。
探索の途中で、都度、途中までの値と、残りの項で作れる最大の数（配列 u）とを比較し、後者の方が小さければ、残りの項でどうやってもゼロにすることは不可能とみなします。

def search(x, n, i, s, v, u):
    if i == len(v):
        if x == 0: s.add(n)
        return
    if abs(x) > u[i]: return # 探索を打ち切り
    for a in range(10):
        if i == 0 and a == 0: continue
        search(x + v[i] * a, 10 * n + a, i + 1, s, v, u)

def g(m):
    s, d = set(), 1
    while True:
        r = 1
        while True:
            # v[i]: 等式の係数, u[i]: 第i項以下で作れる最大の数
            v, u = [0]*d, [0]*d
            for i in range(d):
                v[i] = r * 10**(d - i - 1) - 16**(d - i - 1)
            if v[0] > 0: break
            for i in range(d):
                u[d - i - 1] = 10 * v[d - i - 1]
                if i != 0: u[d - i - 1] += u[d - i]
            for i in range(d):
                u[i] = abs(u[i])
            search(0, 0, 0, s, v, u)
            r += 1
        if len(s) >= m:
            return sorted(list(s))[m-1]
        d += 1

print g(100)

短時間勝負のコンテストでは、この手の全探索＋枝刈りの問題はほとんど出題されないかもしれません。
ゆえにか結構難しめだったと思います。

2018-03-21

第４問「ラディカル和」問題解答

「ヒストリー和」問題の解法です。
問題文はこちら。
他の問題の解法は、第１問、第２問、第３問、第５問。

多数の m について、m²＋1 の素因数の積を求める問題です。
素朴にm²＋1 の素因数をひとつひとつ求めていたのでは時間がかかりすぎます。

ここでは篩を使った解法を紹介します。
同じアイデアは CodeIQ の過去の問題（数学の問題をプログラミングで解こう！「タンジェント・フラクション」問題解説｜CodeIQ MAGAZINE）でも紹介しました。

2 つの配列を用意しましょう。一方は m²＋1 で初期化し、もう一方は 1 で初期化します。
f:id:riverplus:20180320235142p:plain

上側の配列を左端から見ていきます。1 でない値に遭遇したとき、その値 p は実は必ず素数です。
そしてそこから右に p おきに見た値は、すべて p の倍数となっています。
下図では、m = 1 の値である 2 おきに、上側の配列の値を 2 で割り尽くしていきます。
そして下側の配列に 2 をかけておきます。
f:id:riverplus:20180321000138p:plain

さらに右隣を見ます。隣の値の 5 についても、同様に、配列の値を 5 で割り尽くしていきます。
f:id:riverplus:20180320235147p:plain

その右隣も同様です。
f:id:riverplus:20180320235149p:plain

このように配列をどんどん右に見ていくことで、非常に高速に m²＋1 の値を割っていけます。
こうして最終的に、下の配列に F(m²＋1) の値が並びます。

コード例は以下になります。

#include <iostream>
#include <vector>
using namespace std;
typedef int		s32;
typedef long long int	s64;

s64 G(s64 m) {
	// 2つの配列をそれぞれ1とn^2+1で初期化
	vector<s64> v(m + 1, 1), u(m + 1);
	for (s64 i = 1; i <= m; i++) u[i] = i * i + 1;
	for (s32 i = 1; i <= m; i++) {
		s64 p = u[i];
		if (p == 1) continue;
		for (s64 k = i; k <= m; k += p) {
			v[k] *= p;
			while (u[k] % p == 0) {
				u[k] /= p;
			}
		}
	}
	s64 answer = 0;
	for (s32 i = 1; i <= m; i++) answer += v[i];
	return answer;
}

int main(int argc, char *argv[])
{
	s64 m;
	cin >> m;
	cout << G(m) << endl;
}

これで ideone で実行時間 0.4 秒です。
あまり m の上限を緩くすると、愚直解＋コンパイル言語でちょっと待てばクリアできてしまうので、ギリギリ大きめの値にしました。

2018-03-19

第３問「ヒストリー和」問題解答

「ヒストリー和」問題の解法です。
問題文はこちら。
他の問題の解法は、第１問、第２問、第４問、第５問。

漸化式にしたがって F_k(n) の項の値を求める問題です。
求めるべき答えは、(k, n) = (10, 10¹⁰) と (10⁴, 10⁶) に対する値の和です。
やけに奇妙な組み合わせだな、と思われたかもしれません。
実は本問は、2 つの場合でそれぞれ異なる解法が必要です。

(k, n) = (10, 10¹⁰) のとき

F_k(n) の値は、直前の 10 項の値によって決まります。
行列で書くとこういうことです。
$\left( \begin{array}{l} f_{10}(n+1)\\ f_{10}(n+2)\\ f_{10}(n+3)\\ f_{10}(n+4)\\ f_{10}(n+5)\\ f_{10}(n+6)\\ f_{10}(n+7)\\ f_{10}(n+8)\\ f_{10}(n+9)\\ f_{10}(n+10) \end{array} \right) = \left( \begin{array}{cc} 0&1&0&0&0&0&0&0&0&0\\ 0&0&1&0&0&0&0&0&0&0\\ 0&0&0&1&0&0&0&0&0&0\\ 0&0&0&0&1&0&0&0&0&0\\ 0&0&0&0&0&1&0&0&0&0\\ 0&0&0&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&0&0&1&0\\ 0&0&0&0&0&0&0&0&0&1\\ 1&2&3&4&5&6&7&8&9&10 \end{array} \right) \left( \begin{array}{c} f_{10}(n)\\ f_{10}(n+1)\\ f_{10}(n+2)\\ f_{10}(n+3)\\ f_{10}(n+4)\\ f_{10}(n+5)\\ f_{10}(n+6)\\ f_{10}(n+7)\\ f_{10}(n+8)\\ f_{10}(n+9) \end{array} \right)$

行列をべき乗する典型的な問題に帰着できます。
計算量は O(k³ log n) となります。

(k, n) = (10⁴, 10⁶) のとき

上記のべき乗の方法でこちらの組を解こうとすると、k が大きすぎてうまくいきません。
まったく別のアプローチが必要になります。

F_k(n) だけでなく、直近 k 項の和として G_k(n) を新たに定義することがポイントです。

$\begin{eqnarray} G_k(n)=\sum_{i=1}^{k} G_k(n-i) \\ \end{eqnarray}$

G_k(n) を使うと、例えば k = 5 の場合、F_k(n) の漸化式は次のように表せます。

$\begin{eqnarray} F_5(n) &=& 5F_5(n-1)+4F_5(n-2)+3F_5(n-3)+2F_5(n-4)+F_5(n-5)\\ &=& 5F_5(n-1) + \left( 5F_5(n-2)+4F_5(n-3)+3F_5(n-4)+2F_5(n-5)+F_5(n-6)\right) - \left( F_5(n-2)+F_5(n-3)+F_5(n-4)+F_5(n-5)+F_5(n-6)\right)\\ &=& 5F_5(n-1) + F_5(n-1) - G_5(n-1)\\ &=& 6F_5(n-1) - G_5(n-1) \end{eqnarray}$

一般的に書くとこういうことです。

$\begin{eqnarray} F_k(n) = (k+1)F_k(n-1) - G_k(n-1) \end{eqnarray}$

G_k(n) の値を求めるには、まず長さ k の配列を用意します。
これをリングバッファとして用いて、直近 k 個の値を保持しておきます。
さらにこの和を別に変数として持っておいて、都度更新するようにすれば、n 回程度の計算でG_k(n) の値が求められます。

計算量は O(n) です。

以上の 2 通りの解法を切り替えれば、最終的な答えが出せます。コード例は以下になります。

D = 100000000

def pow_mod(x, m):
    if m == 1: return x
    if m % 2 == 0:
        p = pow_mod(x, m / 2)
        return mul(p, p)
    else:
        p = pow_mod(x, m - 1)
        return mul(x, p)

def mul(x, y):
    z = [[0 for i in range(len(y[0]))] for j in range(len(x))]
    for i in range(len(x)):
        for j in range(len(y[0])):
            for k in range(len(x[0])):
                z[i][j] = (z[i][j] + x[i][k] * y[k][j]) % D
    return z

def f_matrix(n, k):
    m = [[] for i in range(k)]
    for i in range(k-1):
        for j in range(k):
            m[i].append(1 if i+1==j else 0)
    for j in range(k):
        m[k-1].append(j+1)
    p = pow_mod(m, n+1)
    return p[k-1][0]

def f_buffer(n, k):
    h, f, g = [0]*k, k, 1
    p = 0
    h[0], h[1] = 1, k
    for i in range(2, n+1):
        f_new = f*(k+1) - g
        g_new = g - p + f
        p = h[i%k]
        h[i%k] = f_new
        f, g = f_new % D, g_new % D
    return f

print (f_matrix(10000000000, 10) + f_buffer(1000000, 10000)) % D

2018-03-19

第２問「奇数和平方数」問題解答

「奇数和平方数」問題の解法です。
問題文はこちら。
他の問題の解法は、第１問、第３問、第４問、第５問。

連続する奇数の和が平方数になるものを探す問題です。
1 から m 番目の奇数の和については、有名な公式が存在します。イコール m² です。
ただ今回は、1 以外の奇数から始まるものについても考慮しなければなりません。
始点の一つ前の奇数を p 番目の奇数、終点を q 番目の奇数とするならば、この間にある奇数の和は q² - p² となります。
これが平方数になるということは・・要するに、ある自然数 r が存在して、

$\begin{eqnarray} p^2 - q^2 = r^2 \\ \Leftrightarrow q^2 + p^2 = p^2 \end{eqnarray}$

となる (q, r, p) の組、すなわちピタゴラス数を求める問題に帰着されます。

ピタゴラス数の列挙といえば定番ネタです。
ピタゴラスの定理 - Wikipedia の記事に、原始ピタゴラス数（q, r, p が互いに素であるピタゴラス数）を列挙するやり方が載ってます。互いに素な自然数の組 (s, t) で、かつ s＞t、s－t が奇数であるものについて、

$\begin{eqnarray} (s^2-t^2, 2st, s^2+t^2) または (2st, s^2-t^2, s^2+t^2) \end{eqnarray}$

は必ず原始ピタゴラス数になります。これでまずは原始ピタゴラス数を列挙しましょう。

題意を満たすのは原始ピタゴラス数だけではありません。実際には、これの定数倍もすべてが題意を満たします。そうしたものの個数は、p 番目の奇数 2p-1 が n 以下であることに注意して、(n+1)/p/2 であると分かります。

コード例は上の通りです。q と r の入れ替わりを考慮してループ内で足す数を 2 倍している点と、最後に p = 0 の場合を足しこんでいる点に注意して下さい。

def gcd(x, y):
	return gcd(y % x, x) if x > 0 else y

def f(n):
    answer = 0
    for s in range(1, n):
        if s*s > n: break
        for t in range(1, s):
            p = s*s + t*t
            if p > n: break
            if (s - t) % 2 == 0: continue
            if gcd(s, t) > 1: continue
            answer += 2*((n+1)/p/2)
    answer += (n+1)/2
    return answer

print f(10000000)

このコードで ideone で 3 秒程度です。

ちなみに原始ピタゴラス数を列挙するには、魔法の行列を使う、もっと高速な方法があります。
ピタゴラスの数を生む行列を参照して下さい。

こちらは C++ で書きました。C++ の速さも相まって、ideone で 0.1 秒切りますね。

#include <iostream>
using namespace std;
typedef long long int	s64;

s64 Do(s64 q, s64 r, s64 p, s64 n) {
	s64 a = 2*((n+1)/p/2);
	if (a == 0) return 0;
	a += Do(q-2*r+2*p, 2*q-r+2*p, 2*q-2*r+3*p, n);
	a += Do(q+2*r+2*p, 2*q+r+2*p, 2*q+2*r+3*p, n);
	a += Do(-q+2*r+2*p, -2*q+r+2*p, -2*q+2*r+3*p, n);
	return a;
}

int main() {
	s64 n = 10000000;
	s64 r = Do(3, 4, 5, n) + (n+1)/2;
	cout << r << endl;
}

2018-03-19

第１問「メダルツアー」問題解答

「メダルツアー」問題の解法です。
問題文はこちら。
他の問題の解法は、第２問、第３問、第４問、第５問。

拾うメダルの個数の期待値を求める問題です。
w, h, m が大した大きさではないので、いろんな方法で解くことができます。

一つは、「和の期待値は期待値の和」のアプローチです。
この言葉の意味については、和の期待値は期待値の和 | 高校数学の美しい物語などを参照して下さい。
メダルの個数の和の期待値を求めるには、それぞれのメダルについて、そのメダルを拾う確率を求めて、和をとればよいのです。

自然数 i（1 ≦ i ≦ h）について、下から i 番目のメダルを考えましょう。
このメダルを拾う確率はいくらでしょうか。
下のような図を書いて考えましょう。

f:id:riverplus:20180319005019p:plain

スタート（A）からゴール（D）までの行き方は全部で ${}_{w+h} C _w$ 通りです。
このうち、A→B→C→D の順で通れば、i 番目のメダルを拾うことができます。
そのような行き方は、 ${}_{m+i-2} C _{m-1} \times {}_{w+h-m-i+1} C _{h-i}$ 通りです。
この 2 つの値の比が、i 番目のメダルを求める確率です。

i についてループを書けばOKです。
コード例は次の通りです。分母の部分は i には依存しないため、最後に割っています。

def bi(a, b):
    r = 1
    for i in range(b):
        r *= a-i
        r /= i+1
    return r

def f(w, h, m):
    answer = 0
    for i in range(1, h+1):
        answer += bi(m+i-2, m-1) * bi(w+h-m-i+1, h-i)
    return float(answer) / bi(w+h, w)

print f(10, 10, 3)

　
　
もうひとつ、別解を紹介します。

def f(w, h, m):
    return float(h) / (w+1)

print f(10, 10, 3)

・・実はこの問題、答えは h/(w+1) なんです。
O(1) です。
m の値に依存しません。
ビックリしますよね。

出題者もビックリしました。
この記事で証明を書こうとしたのですが、結局ギブアップしました。ごめんなさい。

追記

その後、まさにズバリな解法を教えてもらいました！ Vandermonde の恒等式という打ってつけがあるのですね！

第１問「メダルツアー」問題解答 - riverplus のブログ https://t.co/dggf3MaGdN
By using Vandermonde's identity : \sum_i {}_{m+i-2} C_{m-1} \times {}_{w+h-m-i+1} C_{h-i} = {}_{w+h_1} C_{h-1} which yields, divided by {}_{w+h} C_h, the surprizing formula \frac{h}{w+1}
Wow!
— TJ Takei (@karoyakani) 2018年3月19日

そしてこちらは↑と→の並び替えで捉え、これまたズバリで期待値を得る考え方です。

書きました。https://t.co/TdFDDm7Sv3 https://t.co/h4Zww8lgJK
— haruya (@haruya1212) 2018年3月19日

ありがとうございました！

2018-02-25

数学プログラミング・チャレンジ（問題）

数学プログラミングの問題を５問出題します。
いずれもCodeIQで出題しようと思っていたネタです。
各問とも正解者用のひとこと掲示板を用意してます。一番乗りを目指して下さい！

ルールなど

答えの数字をそのままフォームに入力して正解チェックする Project Euler 形式の出題です。
難易度は Project Euler の Easy～Medium ぐらいです。
うまくコードを書くと、普通の性能の PC で 1 分以下の実行時間で答えが出せるよう設計しています。出題者は C++ または Python で ideone で 1 秒を切れることを確認しています。
本問については、ブログ記事など自由に書いて下さってＯＫです。ただとりあえずの解禁日ということで 3/18(日) AM 8 時以降に公開して頂けると嬉しいです。
問題文についてのお問い合わせはこの記事のコメント欄か @riverplus までお願いします。

それではどうぞチャレンジを！

第１問：「メダルツアー」問題

自然数 w, h に対し、横と縦の長さがそれぞれ w, h となる格子状の道を考えます。
さらに、自然数 m に対し、左から m 番目の縦の道の、各交差点の中間にメダルが 1 個ずつ置いてあります。

例えば下図は、w, h, m = (4, 3, 2) のときの様子です。

この図の左下の点からスタートして、右上のゴールまで最短距離で向かいます。
途中で通過したメダルはすべて拾います。

取り得る経路を等しい確率で選ぶとき、拾うメダルの個数の期待値を F(w, h, m) と定義します。

例えば F(3, 2, 2) = 0.5 です。
取り得る経路は下図の 10 通りで、そのうちメダルを 2 個拾うものが 1 通り、1 個拾うものが 3 通り、0 個拾うものが 6 通りです。
したがって、期待値は 2×(1/10)＋1×(3/10)＋0×(6/10) ＝ 0.5 と計算できます。

同様に、F(1, 1, 1)＝0.5, F(2, 1, 3)＝1/3, F(4, 3, 2)＝0.6 となることが確かめられます。

小数点以下が 7 桁になるよう四捨五入した F(10, 10, 3) の値を求めて下さい。

正解チェックはこちらから

第２問：「奇数和平方数」問題

自然数 n に対し、連続する n 以下の奇数の和が平方数（自然数の 2 乗で表される数）となるものを探します。

例えば n = 10 の場合、以下の 7 通りです。
項が 1 つだけのものもカウントしている点に注意して下さい。

　　1 ＝ 1²
　　1＋3 ＝ 2²
　　1＋3＋5 ＝ 3²
　　9 ＝ 3²
　　1＋3＋5＋7 ＝ 4²
　　7＋9 ＝ 4²
　　1＋3＋5＋7＋9 ＝ 5²

このような表し方の数を F(n)と定義します。例えば F(10) = 7 です。
同様に、F(20)＝14, F(30)＝23, F(1000)＝1272 となることが確かめられます。

F(10⁷) の値を求めて下さい。

正解チェックはこちらから

第３問：「ヒストリー和」問題

自然数 k に対して、数列 F_k(n) を次のように定義します。
$\begin{eqnarray} F_k(n)=\left\{ \begin{array}{ll} 0 & (n<0) \\ 1 & (n=0) \\ \sum_{i=1}^{k} (k-i+1) F_k(n-i) & (n>0) \\ \end{array} \right. \end{eqnarray}$

例えば k = 3 なら次のようになります。
$\begin{eqnarray} F_3(n)=\left\{ \begin{array}{ll} 0 & (n<0) \\ 1 & (n=0) \\ 3F_3(n-1)+2F_3(n-2)+F_3(n-3) & (n>0) \\ \end{array} \right. \end{eqnarray}$

例えば F₃(1) = 3, F₃(2) = 11, F₃(3) = 40, F₃(4) = 145, F₅(10) = 36709715 です。

F₁₀(10¹⁰) + F_10⁴(10⁶) の値の、下 8 桁を求めて下さい。

正解チェックはこちらから

第４問：「ラディカル和」問題

自然数 n に対し、n の互いに異なる素因数の積を F(n) と定義します。

例えば F(24) = 6 です。24 を素因数分解すると 24 = 2³・3¹ で、互いに異なる素因数は 2 と 3 だからです。
同様に、F(5) = 5、F(60) = 30、F(1024) = 2 となることが確かめられます。

また、自然数 m に対し、1 から m までの自然数 n に対する F(n² + 1) の和を G(m) と定義します。

例えば G(3) = F(1² + 1) + F(2² + 1) + F(3² + 1) = F(2) + F(5) + F(10) = 2 + 5 + 10 = 17 です。
同様に、G(7) = 107、G(20) = 2590、G(100) = 299434、G(1000) = 303647368 となることが確かめられます。

G(2×10⁶) の値を求めて下さい。

正解チェックはこちらから