予選課題: 長い桁数の答えを求める

東京工業大学 SuperCon'98 問題策定委員会
1998.6.1

問題

整数を係数に持つ３次方程式の解(実数根)を、小数点以下400桁まで 求める。問題の方程式を

x3 - 6 x2 + 8 = 0
とすると、この方程式は３つの実数根α, β, γ (α < β < γ) を持つ。ここで、
-2 < α < -1、　 1 < β < 2、 　5 < γ < 6
である。根β を十進数で小数点以下400桁まで標準出力に書き出す プログラムを作成せよ。なお、小数点部分は横50桁単位で改行すること。

ヒント

問題を解く鍵は3つある。1) 根を近似的に求めるためのアルゴリズムの選択、 2) 長い桁数(多倍長)の計算をどうするのか、3)計算機上で効率よく計算するための 手法の活用、である。以下、それぞれについて簡単なヒントを与える。

1)近似解アルゴリズム

f(x) = x3 - 6 x2 + 8
とすると、f(β)=0である。したがって、ある規則で、xの値を変化させながら、 f(x)=0となるxを求めれば、その答えは β となりゆる。xの値の列を x1 , x2 , ... とすると、
xn+1 = xn - f(xn)/A (A ≠ 0, Aは適当な定数)
が知られている(直感的には A は、この関数のある地点での傾きになっている)。 Aおよびx1 が適切に選べれば、これを無限に繰り返して得られる x∞ が求める β の値である(ニュートン法)。

ただし、これでは永久に計算が終わらないので、適当な精度が得られたところで 計算を打ち切る必要がある。打ち切るための目安としては、

|xn - x∞| ≦ K(1-ε)n
を利用する。ε に関しては、正定数であって
ε ≦ f'/A ≦ 2 - ε
f' : f の１次微分(= 3 x2 - 12 x) で、 x=xi のときの値
となるように選ぶことで、x∞が保証される。なお、Kは、1程度の正数。 この場合、A, x1 の選択は、近似の精度と収束の速さに影響する。

Aを、f'で置き換えるアプローチも存在する(ニュートン・ラプソン法)。すなわち、

xn+1 = xn - f(xn)/f'(xn)
f'は直感的には、傾きを表わす。この関係がうまくいくためには、 f'がx∞の 近傍で連続(f'(x∞) ≦ 0)であり、かつ、 |x∞ - x1| が十分小さいことが 必要である。

この手法での収束の速度は、

|xn - x∞| ≦ B exp(C 2n)
で与えられる。ただしB, Cはある定数で、C < 0である。

2) 長い桁数の計算(多倍長計算)

多倍長の足し算 C = A + B を考えてみる。A, B, Cは整数配列とし、それぞれの 値は、配列要素の列で表わすことにする。たとえば、A=123456789123ならば、

A[0]=1234, A[1]=5678, A[2]=9123 (10進数で1要素あたり4桁格納すると仮定)
と表わせる。実際の計算は、
C[i] = A[i] + B[i], 　 i=N, ..., 0 (Nは最大桁数)
を繰り返しながら、桁上がりが生じるたびに、必要な処理を行えばよい。

引き算は、基本的には足し算と同じである。

掛け算も、通常の掛け算と同じ考え方が使える。ただし、各配列要素同士で 掛け算した結果が、配列要素で表現できる数の範囲を超えてはいけない。 また、要素同士の掛け算の結果が決められた桁数(今の場合4桁)を超えた場合の 処理も必要である。

割り算は、除数が4桁以下の数である時は、直前の要素の計算の余りを考慮することで、 通常の割り算と同じ考え方が使える。しかしながら、多倍長同士の割り算の場合には、 除数の逆数を、被除数に乗じる方が計算速度が速くなる。その場合、正の除数 b を

b = B0(1 - δ)
ただし、B0 は4桁以下の自然数 × 10整数 、 　0 < δ < 1、 　 δ は長桁で表現される数を小数点数以下で表わした数である。
とする。この時、除算 a/bは次のように表わされる。
a/b = a/B0 Πk=0∞ (1 + δ2k )
ここで
Πk=0∞ ...
は続く式の 積を意味する。 (計算の都合上、δ2kは (δ2k-1)2 を利用し、 乗積は必要な精度の項まで用いる)

なお、今の場合、10進表現を仮定しているが、2進表現で多倍長計算を行うことも 可能である。

3) 計算速度向上のためのヒント

本問題のように、繰り返し計算がある場合、アルゴリズム上の繰り返し回数を 減らすことが重要である。その一方で、答えに到達するまでの、収束速度も 重要である。また、簡単な掛け算、割り算は、足し算、引き算に置き換えることも 可能である。さらに、除算を除数の逆数で乗じる場合の公式も活用できる。

いまの問題の場合、xn からxn+1を求める場合、 単純に前述の式から計算する こともできるが、例えば次のような形で求めることもできる。

xn+1 = xn2 /6 + 4/3 xn あるいは、
tn+1 = 3/4 tn - 1/8 tn2, 　 x∞ = 1/t∞ あるいは、
xn+1 = √{(xn3 + 8)/6}
なお、√ の計算には工夫を要する。
また、逆数計算には次の式が利用できる。
1/(1 -δ) = Πk=0∞ (1 + δ2k), 　(|δ| < 1) あるいは、
1/(1 +δ) = (1 - δ)Πk=0∞ (1 + δ2k), 　(|δ| < 1)

参考

小数点以下50桁を示す。

30540728933227860459313349292274081599849729126372

審査のために使用する計算機の int,longのサイズはそれぞれ 32bit,64bitとする。

本選予告

「板詰め問題」 与えられた板を組み合わせて、四角形を埋める。