ガジェット通信 GetNews

見たことのないものを見に行こう
ワンダーウーマン

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

DATE:
  • ガジェット通信 GetNewsを≫

問題のおさらい

一辺の長さが 1 の正方形を考えます。これを P0 と呼びます。
各頂点を時計回りに A, B, C, D とします。
P0 に対し、次の手順で新しい正方形を描きます。

辺 BC と同じ辺を持つ正三角形 BCE を P0 の外側に描きます。
D と E を線分で結びます。
そして線分 DE を一辺とする正方形を右側に描きます。この正方形を P1 と呼びます。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

同じ作業を P1 に対し行い、正方形 P2 を描きます。
以降同様に、P3, P4, … と次々と正方形を描くことができます。
P3 までを描いた様子を下図に示します。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

自然数 n に対し、Pn の面積の値を小数点以下で切り捨てた値を F(n) と定義します。

例えば F(1)=3 です。P1 の面積はおよそ 3.732 だからです。
同様に、F(2)=13, F(3)=51 となることが確かめられます。

標準入力から、自然数 n (1 ≦ n ≦ 106)が与えられます。
標準出力に F(n) を 106 で割った余りを出力するプログラムを書いてください。

方針:まずは定式化

平面幾何に関する問題です。本問は、まずは F(n) の公式を求めることから始まります。正方形 Pn から正方形 Pn+1 を作る様子に注目します。点 E から辺 AD に下した垂線の足を F とし、三角形 DEF に着目しましょう。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

Pn の一片の長さを Ln とします。DF の長さは、正方形の一片の長さの半分ですから、Ln/2 です。また、EF の長さは、一片の長さが Ln の正三角形の高さと、正方形の一片の長さの和ですから、つまり (1+(√3)/2)・Ln です。したがって、三平方の定理より、以下の関係が成り立ちます。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

さて、正方形 Pn+1 の面積は DE の長さの 2 乗です。ですから、F(n) と F(n+1) の間に次の関係が言えます。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

要は、正方形が一つ大きくなるごとに、面積は (2+√3) 倍になるわけです。P0 の面積は 1 ですから、次のようになります。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

つまり本問は、(2+√3)n の整数部分を求めて、それを 106 で割った余りを求める問題だと言い換えられます。

素直な実装では精度に困難あり

さて、まずは素直に実装してみましょう。次のようなコードになると思います。

import math

n = int(input())
f = math.pow(2 + math.sqrt(3.0), n)
print int(f) % 1000000

この実装は、n が比較的小さい場合はうまくいきます。しかし、n が大きくなるとエラーになってしまいます。値が大きくなりすぎて、ふつうの浮動小数点型の変数で扱えなくなるためです。

どうすればよいでしょうか? ひとつの考えは、高精度に値を保持できるライブラリを使うことです。が、本問のように、n が最大で 106 にもなる状況では、(2+√3)n の値を十分な精度で保存するのは、メモリ容量的に厳しいでしょう。

解法の手がかり―共役な値とセットで考えよう

そこで次のように考えます。まず、整数 pn, qn を用いて (2+√3)n を次のように表すことが本問の一つ目のポイントです。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

pn, qn を求めるには、次のように漸化式を立てます。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

行列で表すと次のようになります。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

行列を 1 回作用させるごとに次の n になるわけですから、列ベクトル (p0, q0) = (1, 0) から始めて、行列の n 乗を作用させればOKです。後述しますが、これにはよく知られた計算方法が存在します。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

ただ、pn や qnを求めるだけでは本問は解けません。2つ目のポイントは、共役を考える ことです。つまり、(2-√3)n もセットに考えます。(2+√3)n のときと同じ計算をやってみましょう。すると面白いことに、(2+√3)n のときと同じ pn, qn を使って、(2-√3)n は次のように表せるのです。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

二つの式の左辺と右辺同士を足します。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

(2-√3)n を右辺に移行します。

数学の問題をプログラミングで解こう!「ルーム・アンド・ルーフ」問題解説

これで全ての準備が整いました。(2+√3)n の整数部分はどのように求められるでしょうか? 右辺をじっと見ましょう。(2-√3)n のところです。この部分、n が大きくなると、きわめて 0 に近い値になることに気が付いたでしょうか。2-√3 ≒ 0.268 です。これを何乗かしていくと、どんどん小さくなります。とはいえ決して 0 にはなりません。

(2+√3)n とは、2 pn(これは整数)から、そのような微小な正の値を引いたものなんですね。つまり、(2+√3)n の整数部分とは、2 pn から 1 を引いた値だということです。

実装

これで実装ができます。まずは pn を計算しましょう。これは漸化式をふつうに計算すればOKです。あるいは下記のコード例のように、行列のべき乗を効率よく行うアルゴリズムを使ってもよいでしょう。

最後に、行列の (1, 1) 要素を 2 倍して、1 を引けば、それが求めるべき答えです。本問では 106 での余りさえ分かればよいので、適宜余りの計算を挟んでいます。また (1, 1) 要素が 0 になる可能性を考慮して、1 を引いてから 106 を足した後に余りの計算をしています。コード例は以下です。

D = 1000000

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

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

matrix = [
[2, 3],
[1, 2]
]

n = int(input())
p = pow(matrix, n)
print (2*p[0][0] – 1 + D) % D

みんなのコード

本問の挑戦者で、コードを公開して頂いた方のコードを Togetter でまとめました。(随時追加していきます。)こちらもぜひチェックしてください。

CodeIQ「ルーム・アンド・ルーフ」 みんなのコード

CodeIQ運営事務局より

Kawazoeさん、ありがとうございました!
現在、Kawazoeさんの最新問題が出題中です。
ぜひ挑戦してみてくださいね!

カテゴリー : デジタル・IT タグ :
CodeIQ MAGAZINEの記事一覧をみる ▶
  • 誤字を発見した方はこちらからご連絡ください。
  • ガジェット通信編集部への情報提供はこちらから
  • 記事内の筆者見解は明示のない限りガジェット通信を代表するものではありません。