第5問「10/16数」問題 解答
「10/16数」問題の解法です。
問題文はこちら。
他の問題の解法は、第1問、第2問、第3問、第4問。
自然数 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 の各桁を a1a2・・ad と表します。
さらに F(n)÷n の値を r とします。
例えば d = 4,r = 3 のとき、a1, a2, a3, a4 が満たすべき等式は何でしょうか。
それは次の式です。
整理するとこうなります。
a1 の係数が大きいことがポイントです。
例えば a1 = 1 としてみます。第 1 項の値は 1096 です。
a2, a3, a4 は 0~9 の自然数ですから、左辺の第 2 項以降は、どれだけ小さくとも -64・9-14・9-2・9 = -720 にしかなりません。第 1 項の値を打ち消すには足りません。
したがって、具体的な a1, a2, a3, a4 の値を実際に試してみるまでもなく、この等式は決して解を持たないと分かります。
このように、上の桁(a1)から数字を仮決めして、以下の項でどうやっても等式を成立(辺の値をゼロにする)させられないと分かれば、以降の探索を中断する、という枝刈りのアプローチです。
コード例は以下です。ちょっとややこしいです。
始めに桁数 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)
短時間勝負のコンテストでは、この手の全探索+枝刈りの問題はほとんど出題されないかもしれません。
ゆえにか結構難しめだったと思います。
第4問「ラディカル和」問題 解答
「ヒストリー和」問題の解法です。
問題文はこちら。
他の問題の解法は、第1問、第2問、第3問、第5問。
多数の m について、m2+1 の素因数の積を求める問題です。
素朴にm2+1 の素因数をひとつひとつ求めていたのでは時間がかかりすぎます。
ここでは篩を使った解法を紹介します。
同じアイデアは CodeIQ の過去の問題(数学の問題をプログラミングで解こう!「タンジェント・フラクション」問題解説|CodeIQ MAGAZINE)でも紹介しました。
2 つの配列を用意しましょう。一方は m2+1 で初期化し、もう一方は 1 で初期化します。
上側の配列を左端から見ていきます。1 でない値に遭遇したとき、その値 p は実は必ず素数です。
そしてそこから右に p おきに見た値は、すべて p の倍数となっています。
下図では、m = 1 の値である 2 おきに、上側の配列の値を 2 で割り尽くしていきます。
そして下側の配列に 2 をかけておきます。
さらに右隣を見ます。隣の値の 5 についても、同様に、配列の値を 5 で割り尽くしていきます。
その右隣も同様です。
このように配列をどんどん右に見ていくことで、非常に高速に m2+1 の値を割っていけます。
こうして最終的に、下の配列に F(m2+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 の上限を緩くすると、愚直解+コンパイル言語でちょっと待てばクリアできてしまうので、ギリギリ大きめの値にしました。
第3問「ヒストリー和」問題 解答
「ヒストリー和」問題の解法です。
問題文はこちら。
他の問題の解法は、第1問、第2問、第4問、第5問。
漸化式にしたがって Fk(n) の項の値を求める問題です。
求めるべき答えは、(k, n) = (10, 1010) と (104, 106) に対する値の和です。
やけに奇妙な組み合わせだな、と思われたかもしれません。
実は本問は、2 つの場合でそれぞれ異なる解法が必要です。
(k, n) = (10, 1010) のとき
Fk(n) の値は、直前の 10 項の値によって決まります。
行列で書くとこういうことです。
行列をべき乗する典型的な問題に帰着できます。
計算量は O(k3 log n) となります。
(k, n) = (104, 106) のとき
上記のべき乗の方法でこちらの組を解こうとすると、k が大きすぎてうまくいきません。
まったく別のアプローチが必要になります。
Fk(n) だけでなく、直近 k 項の和として Gk(n) を新たに定義することがポイントです。
Gk(n) を使うと、例えば k = 5 の場合、Fk(n) の漸化式は次のように表せます。
一般的に書くとこういうことです。
Gk(n) の値を求めるには、まず長さ k の配列を用意します。
これをリングバッファとして用いて、直近 k 個の値を保持しておきます。
さらにこの和を別に変数として持っておいて、都度更新するようにすれば、n 回程度の計算でGk(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
第2問「奇数和平方数」問題 解答
「奇数和平方数」問題の解法です。
問題文はこちら。
他の問題の解法は、第1問、第3問、第4問、第5問。
連続する奇数の和が平方数になるものを探す問題です。
1 から m 番目の奇数の和については、有名な公式が存在します。イコール m2 です。
ただ今回は、1 以外の奇数から始まるものについても考慮しなければなりません。
始点の一つ前の奇数を p 番目の奇数、終点を q 番目の奇数とするならば、この間にある奇数の和は q2 - p2 となります。
これが平方数になるということは・・要するに、ある自然数 r が存在して、
となる (q, r, p) の組、すなわちピタゴラス数を求める問題に帰着されます。
ピタゴラス数の列挙といえば定番ネタです。
ピタゴラスの定理 - Wikipedia の記事に、原始ピタゴラス数(q, r, p が互いに素であるピタゴラス数)を列挙するやり方が載ってます。互いに素な自然数の組 (s, t) で、かつ s>t、s-t が奇数であるものについて、
は必ず原始ピタゴラス数になります。これでまずは原始ピタゴラス数を列挙しましょう。
題意を満たすのは原始ピタゴラス数だけではありません。実際には、これの定数倍もすべてが題意を満たします。そうしたものの個数は、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; }
第1問「メダルツアー」問題 解答
「メダルツアー」問題の解法です。
問題文はこちら。
他の問題の解法は、第2問、第3問、第4問、第5問。
拾うメダルの個数の期待値を求める問題です。
w, h, m が大した大きさではないので、いろんな方法で解くことができます。
一つは、「和の期待値は期待値の和」のアプローチです。
この言葉の意味については、和の期待値は期待値の和 | 高校数学の美しい物語 などを参照して下さい。
メダルの個数の和の期待値を求めるには、それぞれのメダルについて、そのメダルを拾う確率を求めて、和をとればよいのです。
自然数 i(1 ≦ i ≦ h) について、下から i 番目のメダルを考えましょう。
このメダルを拾う確率はいくらでしょうか。
下のような図を書いて考えましょう。
スタート(A)からゴール(D)までの行き方は全部で 通りです。
このうち、A→B→C→D の順で通れば、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 の恒等式という打ってつけがあるのですね!
第1問「メダルツアー」問題 解答 - riverplus のブログ https://t.co/dggf3MaGdN
— TJ Takei (@karoyakani) 2018年3月19日
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!
そしてこちらは↑と→の並び替えで捉え、これまたズバリで期待値を得る考え方です。
書きました。https://t.co/TdFDDm7Sv3 https://t.co/h4Zww8lgJK
— haruya (@haruya1212) 2018年3月19日
ありがとうございました!
数学プログラミング・チャレンジ(問題)
数学プログラミングの問題を5問出題します。
いずれもCodeIQで出題しようと思っていたネタです。
各問とも正解者用のひとこと掲示板を用意してます。一番乗りを目指して下さい!
ルールなど
- 答えの数字をそのままフォームに入力して正解チェックする Project Euler 形式の出題です。
- 難易度は Project Euler の Easy~Medium ぐらいです。
- うまくコードを書くと、普通の性能の PC で 1 分以下の実行時間で答えが出せるよう設計しています。出題者は C++ または Python で ideone で 1 秒を切れることを確認しています。
- 本問については、ブログ記事など自由に書いて下さってOKです。ただとりあえずの解禁日ということで 3/18(日) AM 8 時以降に公開して頂けると嬉しいです。
- 問題文についてのお問い合わせはこの記事のコメント欄か @riverplus までお願いします。
それではどうぞチャレンジを!
第1問:「メダルツアー」問題
自然数 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) の値を求めて下さい。
正解チェックはこちらから
第2問:「奇数和平方数」問題
自然数 n に対し、連続する n 以下の奇数の和が平方数(自然数の 2 乗で表される数)となるものを探します。
例えば n = 10 の場合、以下の 7 通りです。
項が 1 つだけのものもカウントしている点に注意して下さい。
1 = 12
1+3 = 22
1+3+5 = 32
9 = 32
1+3+5+7 = 42
7+9 = 42
1+3+5+7+9 = 52
このような表し方の数を F(n)と定義します。例えば F(10) = 7 です。
同様に、F(20)=14, F(30)=23, F(1000)=1272 となることが確かめられます。
F(107) の値を求めて下さい。
正解チェックはこちらから
第3問:「ヒストリー和」問題
自然数 k に対して、数列 Fk(n) を次のように定義します。
例えば k = 3 なら次のようになります。
例えば F3(1) = 3, F3(2) = 11, F3(3) = 40, F3(4) = 145, F5(10) = 36709715 です。
F10(1010) + F104(106) の値の、下 8 桁を求めて下さい。
正解チェックはこちらから
第4問:「ラディカル和」問題
自然数 n に対し、n の互いに異なる素因数の積を F(n) と定義します。
例えば F(24) = 6 です。24 を素因数分解すると 24 = 23・31 で、互いに異なる素因数は 2 と 3 だからです。
同様に、F(5) = 5、F(60) = 30、F(1024) = 2 となることが確かめられます。
また、自然数 m に対し、1 から m までの自然数 n に対する F(n2 + 1) の和を G(m) と定義します。
例えば G(3) = F(12 + 1) + F(22 + 1) + F(32 + 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×106) の値を求めて下さい。
正解チェックはこちらから
第5問:「10/16数」問題
10 進数で表された自然数 n に対し、この数字の並びをあたかも 16 進数のように読み、これを 10 進数として表したものを F(n) と定義します。
例えば n = 23 のとき、「23」を 16 進数として読んだものを 10 進数として表すと、2×16+3 = 35 ですから、F(23) = 35 となります。
同様に、F(9) = 9, F(10) = 16, F(11) = 17, F(100) = 256, F(999) = 2457 です。
さて、F(n) が n の整数倍となるような自然数を、「10/16 数」と呼ぶことにしましょう。
例えば 1038 は 10/16 数です。「1038」を 16 進数として読んだものを 10 進数で表すと、1×163+0×162+3×161+8 = 4152 により、F(1038) = 4152 です。たしかに 4152 は 1038 の倍数となっています。
1038 以外に 10/16 数はあるでしょうか。
1 桁の自然数(1~9)は全て 10/16 数です。
その次に小さな 10/16 数は 1038 です。
以降小さい順に、1040, 2078, 2080, 2118, … と続きます。
また 20 番目に小さな 10/16 数は 4238 です。
100 番目に小さな 10/16 数を求めて下さい。
正解チェックはこちらから
「ペア・ドロップ」問題 みんなのコード
「ペア・ドロップ」問題のみんなのコードです。
いつもありがとうございます!!!
CodeIQ「ペア・ドロップ」問題、明日10時で公開終了です。けっこう難しめの確率の問題でした。時間になったら、今回もぜひ皆様のコードを公開して @riverplus まで教えて下さい。Togetterでまとめます! https://t.co/yvHyOrsUKt
— Kawazoe (@riverplus) 2017年12月14日
ペア・ドロップです https://t.co/xETa36Qkd8
— ininsanus (@ininsanus) 2017年12月15日
ペアドロップruby(108)です
— SMZ8110 (@smz_8110) 2017年12月15日
n,k=https://t.co/63Mu9HqRWk &:to_r
r,d=10**6*2**n,n-k
1.upto(n){|i|r/=i+n;i>k||r/=i
i<=d&&i%2<1||r*=i*i-d%2}
p r.to_i
golfじゃない版https://t.co/OPsilnkniU
— SMZ8110 (@smz_8110) 2017年12月15日
(肝心な部分にコメントが付いてません)
身も蓋もない説明しか書いてませんが。https://t.co/XvkDPDZ9XC
— idiotton (@idiotton) 2017年12月15日
コード公開します。https://t.co/1u4VZX4j0t
— haruya (@haruya1212) 2017年12月15日
ペアドロップ解答…。特に捻ったことはしてません…https://t.co/CEucjQVEec
— Nick-IB(和名:高瀬 健) (@Nyagoking) 2017年12月15日
ペア・ドロップのコードネタバレです。身も蓋もないですが、F(n,k)=2^k・(n!)^3/(k!・((n-k)/2)!^2・(2n)!) (n,kの偶奇が一致する場合) をそのまま計算しています。 https://t.co/0m1XPOc5W9
— angel as ㌵㌤の猫 (@angel_p_57) 2017年12月16日
1tweetに収めることもできますね。素直に縮めてRuby(100)です。 n,k=https://t.co/GdVXBGQn0Z &:to_i;p (1000000<<k)*(f=->x{(1..x).reduce 1,&:*})[n]**3/f[k]/f[n-k>>1]**2/f[n*2]
— angel as ㌵㌤の猫 (@angel_p_57) 2017年12月16日