p11190.jpg

上式の通り、ごく簡単な、ベキ級数の合計を算出する問題。変数 l, h, k の値はともに1~15000000、hとlとの差は1000以内。

kの値が10以内なら、探せば、それなりの計算式が出てくるが、1500万ではそんな式は存在しないだろう。

力任せで、C言語の数学関数をそのまま使うと、無限大になってしまい、計算不能になる。

考え方
 合計値は級数の後ろの項によってほぼ決まるので、級数の末尾から計算することにしよう。

大事なことは、個々の項の値を、仮数部と指数部に分けておく。

例えば、999100、998100の値を
  999100 仮数部 0.9047921471 指数部 e+300
  998100 仮数部 0.8185668046 指数部 e+300
に分ける。

指数部の値が1つ違うと、項の値として1桁も違ってくる。指数部の値が5も違うと、5桁も右の部分にしか影響しないので、その後の計算を打ち切っても全体の合計値にほとんど影響がないだろうと考えてよい。例えば、
  880100 仮数部 0.2807160311 指数部 e+295
当たりから、その後の合計計算を打ち切る。基数やkの値が大きければ大きいほどすぐに打ち切ることになる。

また、合計の計算では、級数の末尾から項ひとつひとつ足していくが、個々の項の指数部は、末尾項の指数部の値に(つまり、最大指数部)に合わせて、仮数部の小数点をシフトさせる。
 上の例では、級数の末尾項が999100なら、個々の項の指数部の値をe+300に統一して、仮数部の値を合計する。
  970100 仮数部 0.0475525075 指数部 e+300

なお、シフトする前の仮数部、指数部は以下の数学関数でも求めることができる。

 ベキ乗 xk → 仮数部 pow(10, modf(k*log10(x), &e)) 指数部 左の e

最後に、元問題の検証データを載せておく(入力データの各行は左から、l, h, kの値)

10 10 1000000
14999001 15000100 1
14999001 15000000 15000000
0 0 15000000
0 1 15000000
0 0 1
Case 0001: 0.100000e0001000001
Case 0002: 0.164995e0000000011
Case 0003: 0.121628e0107641370
Case 0004: 0.000000e0000000001
Case 0005: 0.100000e0000000001
Case 0006: 0.000000e0000000001

080907.jpg

16進数100万桁表示の円周率計算プログラムについて考える。実行速度は2秒以内としたい。効率のよい計算式と、FFTによる乗算で高速化にチャンレンジしたい。

ここのブログでも紹介している多倍長整数を利用するが、16進数に合わせて、1語の大きさを2進32ビット、つまり、16進8桁に直す。10進数と違って、除算ではビットシフト、剰余計算ではビットマスクだけでできるので、さらなる高速化が期待できそう。

整数一つの配列はこんなものになるだろう。オーバーフローが起きて、64ビット整数にしないといけない可能性もあるが。

  unsigned num[125000];

Chudnovsky)の公式は計算が速いようだが、確かめないとなんともいえない。

p11302.jpg

ネット上を調べたら、全ての桁をいっぺんに計算するのではなく、特定の桁だけを算出する方法もあるらしい。今回の用途に合致するかもしれない。

英文 Algorithm for Calculating Individual Hexadecimal Digits of Pi

なお、16進円周率の値一部はここにて確認できる。

2次多項式の因数分解について考える。

例: 2x2 + 5x – 12 = (x + 4) (2x – 3)

複素数や無理数にまで拡張したら面白くないので、今回は、2次多項式の係数すべては整数ということに制限する。定理により、整数係数の多項式が有理数体で因数分解できるなら、因数の係数もすべて整数にすることができることから、因数の係数すべても整数とする。

前回の記事 「大きさとの戦い」 では結局解き方が判らなくて、そのまま、中途半端に終わってしまった。

その時はまだ簡単に解けそうな問題が結構残っていたのだが、ここに来ると一問一問解くのに益々時間がかかる気がしてきて、再チャレンジしてみた訳だ。

さて、その記事では、最終的に得られた方程式は

   (2t+1)2 – (2s)2 2/(2n+1) = 1  (n = 0, 1, 2, 3, 4, 5, 6, 7)

というものだった。X = 2t+1とし、n に 0~7の値をそれぞれ代入すると、つぎの8つの方程式が得られる。
   X2 – 8s2 = 1
   X2 – 8/3 s2 = 1
   X2 – 8/5 s2 = 1
   X2 – 8/7 s2 = 1
   X2 – 8/9 s2 = 1
   X2 – 8/11 s2 = 1
   X2 – 8/13 s2 = 1
   X2 – 8/15 s2 = 1

分母を払うために、s にそれぞれ、以下のものを代入する。すると見事にペル方程式にたどり着く。

     s = Y        X2 - 8Y2 = 1          Ans = Y2
s = 3Y       X2 - 24Y2 = 1         Ans = 9Y2
s = 5Y       X2 - 40Y2 = 1         Ans = 25Y2
s = 7Y       X2 - 562 = 1          Ans = 49Y2
s = 3Y       X2 - 8Y2 = 1          Ans = 9Y2
s = 9Y       X2 - 72Y2 = 1         Ans = 81Y2
s = 11Y      X2 - 88Y2 = 1         Ans = 121Y2
s = 13Y      X2 - 104Y2 = 1        Ans = 169Y2
s = 15Y      X2 - 120Y2 = 1        Ans = 225Y2

右列の Ans 部分は、ペル方程式の一般解 (X, Y) のうちの Y の値を、問題の解答に対応させるためのもの。

ペル方程式の解き方は既に 昔の記事 に書いたので、詳細はそちらを参照されたい。

結論的にいうと、方程式の数は9個と多いが、それぞれに簡単に最小解が得られる。さらに、多倍長整数とペル方程式の解の漸化式を使えば、一般解もすぐに判る。

たとえば、2番目のペル方程式 X2 – 24Y2 = 1 の最小解は (5, 1)。一般解は

 Xn+1 = 5 Xn + 24 Yn
 Yn+1 = 5 Yn + Xn
 X1 = 5
 Y1 = 1

になり、最小解(初期解)から漸化的に計算していくと、
 (X, Y) = (5, 1), (49, 10), (485, 99), (4801, 980), …
 Ans = 9Y2 = 9, 900, 88209, 8643600, …
が得られる。

解答をすべてここに書くのはまずいが、条件を満たす数は計 112 個。ラストの10個は以下となっている。

73702684612392858376746072080400
151046383493325234090009219613956
663812918895887474609694358375209
1359417451439927106810082976525604
5131130648390546663702275158894481
8758618149740884101499623707907225
35524541072381125235676320801486025
46180175835514919973320476430050329
60153048103383854522005973075150400
65046891745147247744323041849672900

出題側として、難しい問題にするには、自分や仲間達の最新研究成果を使うが、問題を難しくするには、計算すべき数値の範囲(大きさ)や数値の量を大幅に増やす、といった安易なアプローチもある。

例えば、ごく簡単な足し算。数学では、a と b が与えられたら、その和は a+b で計算できるから、それ以上のことを深く考える必要はない。しかし、コンピュータの世界では、10進20桁未満の整数同士なら簡単だが、それ以上の大きさ、例えば、10進200桁同士の整数 a と b を足すにはひと工夫が必要。

近々、2進128桁までの整数はC言語でもふつうに扱えるようになると予想できるが、ある種の計算ではいままで以上に難しくなる。例えば、素数表の生成には、エラトステネスのふるいが効率がもっともよいのでよく利用されるが、128桁になると、あまり役に立たなくなるかもしれない。10秒では計算が終了しないから。

<オリジナル問題>(素数の分布)
 整数dが与えられたら、素数 a とそのつぎの素数 b との間隔 (=b-a) がd以上であるような、最小素数 a, bペアを計算せよ。ただし、a, bの大きさは10進35桁まで。
  例: d = 1 ⇒ a = 2, b = 3
     d = 2 ⇒ a = 3, b = 5
     d = 3 ⇒ a = 7, b = 11
     d = 4 ⇒ a = 7, b = 11
     d = 5 ⇒ a = 23, b = 29
     d = 6 ⇒ a = 23, b = 29
     d = 7 ⇒ a = 89, b = 97
</オリジナル問題>

ということで、プログラミング問題を解いたといっても、数学と違って、いろいろなレベルがある。つまり、計算量を探求するならば、多くの問題はいまだに解けられていない、と強引に主張しても間違いない。しかし、コンピュータの世界では、そのための言い訳というか、責任逃れの逃げ道が理論的にある程度用意されている。NP完全問題がその好例。誰もが解けないから、私ができなくても当然という、おかしな論理が、コンピュータの世界では堂々と成り立つ。

さて、違う問題をここで考えてみる。

問題 #10241
  三角数 (Triangular Number) とはつぎの数列で表される数のこと、つまり、 1, 3, 6, 10, 15, 21, 28 など、k( k+1)/2 を満たす数のことをいう。一方、平方数(四角数ともいう) (Square number) とは 1, 4, 9, 16, 25 のような数のことをいう。
 平方数を S とし、三角数を T とすると、つぎの条件2つを満たすすべての X を表示せよ。
    X = S
    X = (2n+1) * T  (ただし、n = 0, 1, …, 7)
なお、Xの範囲は 16×1034 まで。
 例: X = 1, 9, 36, 196, 225, 324, 900, 1225, 4225, 11025, 41616, 53361, 88209, 108900, …
</問題 #10241>

この問題の難しいところはやはり大きさの一点につきるだろう。10進36桁の整数はいまのプログラミング環境ではかなり厳しい。

<番外編>
 三角数でもあり、平方数でもある 平方三角数 Square Triangular Number というのがあって、そのN番目の値 Y(N) がつぎの漸化式で求まる。
   Y(N) = 34 * Y(N-1) – Y(N-2) + 2
   Y(1) = 1, Y(2) = 36
</番外編>

ここで考えている問題では、n=0の場合が Square Triangular Number にあたる。

問題解くための数式を整理してみる。T = t (t+1)/2、S = s2 とおくと、
   X = s2
   X = (2n+1) t (t+1)/2
なので、
   t (t+1) (2n+1) / 2 = s2
   (t2+t) (2n+1) / 2 = s2
   ((t+1/2)2-1/4) (2n+1) / 2 = s2
   ((2t+1)2-1) (2n+1) / 8 = s2
   (2t+1)2-1 = 8 s2/(2n+1)
最終的に、
   (2t+1)2 – (2s)2 2/(2n+1) = 1
が得られる。

メモ程度に書いておく。

整数のべき乗、320等、多少高速にする計算。
単純計算では、上の例だと、3を20回乗算すればいいわけだが、以下のようにすると5回の乗算で済む。

    (((32)23)2)2

セコイ考えだが、64ビット整数でもべき乗の指数は大きく取れない(2のべき乗でも、最大64以下だから)ので、実質はあまり意味のないことかもしれない。考え方やC言語プログラムを以下に示す(メモ程度)。

<考え>
  変数名の整理: p = ak  (a、k、p はすべて正の整数)
  1. p に対し、 p ← 1 と初期値を代入。
  2. kの2進数ビットパターンを上位ビットから1ビットずつ見ていき、各ビットに対し、
     (2.1) まず p を二乗する p ← p*p
     (2.2) そのビットが1なら、さらに a を掛ける p ← p*a

<例>
  p = 320 を計算。
  1. p = 1
  2. 20の2進数は 10100 であるので、
     一番左のビット1に対し、
       (2.1) p = p*p = 1*1 = 1
       (2.2) p = p*3 = 1*3 = 3      (=31
     左から2ビット目のビット0に対し、
       (2.1) p = p*p = 3*3 = 9      (=32
     左から3ビット目のビット1に対し、
       (2.1) p = p*p = 9*9 = 81     (= 34
       (2.2) p = p*3 = 81*3 = 243    (= 35
     左から4ビット目のビット0に対し、
       (2.1) p = p*p = 243*243 = 59049   (= 310
     左からの5ビット目のビット0に対し、
       (2.1) p = p*p = 59049*59049 = 3486784401  (= 320
     これで終了。

<計算量>
   指数をNとすると、ふつうの計算では、O(N) であるが、この方法だと、明らかに、O(logN) の計算量で済む。

unsigned long calc_pow(unsigned long a, unsigned k)
{
    unsigned long p;
    unsigned bit;
if (a <= 1) return a;
p = 1; for (bit = 0x80; bit > 0; bit >>= 1) { /* 0x80=128。64bit整数でも十分 */ if (p > 1) p *= p; if (k & bit) p *= a; } return p; }

なお、今の Gnu Cでは64ビット整数 unsigned long long も使えるようになったので、unsigned long のところを合わせて直すといいだろう。

自然対数の底eの値を多桁数で算出する関数。

<関数名>
  e —- 自然対数の底eの値を多桁数で算出する

<形式>
  void e(int *data, int keta);

<引数>
  data (出力) eの値が10進4桁ずつ入った整数の配列
  keta (入力) 小数点以下の桁数の指定

<関数値>
  なし

用例
  int data[1000];
  e(data, 1000);

  e の値 (小数点以下1万桁)

<関数本体>
  e.c

<説明>
  マクローリン展開

ex = 1 + x/1! + x2/2! + … +xn/n! + …

に x = 1 とおき、

     

e = 1 + 1/1! + 1/2! + 1/3! + … + 1/n! + …

の式を利用して計算。さらに、右辺の第2項目以降の計算では、計算回数の少ないつぎの式に変形して使う。

     (((…((1/n + 1) 1/(n-1) + 1) 1/(n-2) + 1)…) 1/3 + 1) 1/2 + 1

円周率 π を精度よく計算する関数。

<関数名>
  pi —- 円周率 π の値を多桁数で算出する

<形式>
  int pi(int *data, int keta);

<引数>
  data (出力) 円周率の値が10進4桁ずつ入った整数の配列
  keta (入力) 小数点以下の桁数の指定

<関数値>
  正常に計算できたときは 0、エラーのときは -1。

用例
  int data[1000];
  pi(data, 1000);

  π の値 (小数点以下1万桁)

<関数本体>
  pi.c

<説明>
  πの値を求めるには、いろいろな式があるが、arctan(逆正接関数)を無限級数で表した式がよく使われる。式が単純な形で扱いやすく、ある程度の桁数の値を得るには適している。しかし欠点は、非常に長い桁数を得るには時間がかかりすぎることである。
  ここで使う式は、最もよく知られている、J. Machinの式
     π/4 = 4 arctan (1/5) – arctan (1/239)
  である。彼は1706年にこの公式でπを小数点以下100桁を求めた。同じ公式で、1949年にはENIACにより2037桁が約70時間で、1958年にはIBM704により1万桁が約100分で求められている。
  arctanをによる古典的な公式にほかに以下のものがある。
  ガウスの公式
     π/4 = 12 arctan (1/18) + 8 arctan (1/57) – 5 arctan (1/239)
  クリンジェンシェルナの公式
     π/4 = 8 arctan (1/10) – arctan (1/239) – 4 arctan (1/515)
  シュテルマーの公式
     π/4 = 6 arctan (1/8) + 2 arctan (1/57) + arctan (1/239)
  高野喜久雄の公式
     π/4 = 12 arctan (1/49) + 32 arctan (1/57) – 5 arctan (1/239) + 12 arctan (1/110443)
  柴田昭彦の公式
     π/4 = 17 arctan (1/22) + 3 arctan (1/172) – 2 arctan (1/682) – 7 arctan (1/5357)

  シュテルマーの公式に8の割り算を工夫すれば高速化期待できるし、高野の公式を利用して並列処理をすれば良い成績が期待できよう。

  arctanを使った公式でπを計算すると、桁数を2倍にするのに計算時間は4倍の増加になるが、以下の方法では2.1〜2.2倍で済む。

  ガウス・ルジャンドルの公式

     A = 1、B = (1/2)1/2、T = 1/4、X = 1

  として、AとBの差が必要とする精度より大きい間、次の計算を繰り返し実行する。

     Y = A、A = (A + B)/2、B = (B * Y)1/2
     T = T – X * (Y-A)2、X = 2 * X

  すると、πの値は(A + B)2/(4 * T) となる。ただし、A、B、T、Yの各値は、求めようとする精度以上で計算しておく必要がある。

1 から 99 の間の平方根の値を多桁数で算出する関数。多くの桁数の平方根が必要な時に利用するといいだろう。

平方根の算出は一般的に難しいけど、結果があっているかどうかの検証はものすごく簡単。2乗してみればすぐに判るからだ。円周率 π とは大違い。

では、関数の紹介。

<関数名>
  squareRoot —- 1 から 99 の間の平方根の値を多桁数で算出する。

<形式>
  int squareRoot(int *data, int keta, int n);

<引数>
  data (出力) 平方根の値が10進4桁ずつ入った整数の配列
  keta (入力) 小数点以下の桁数の指定
  n   (入力) 計算しようとする平方根の2乗の値 (1 〜 99)

<関数値>
  正常に計算できたときは、計算結果となる data の長さ。エラーのときは 0。

用例
  2 の平方根を小数点以下1000桁で求める例
  int data[1000];
  squareRoot(data, 1000, 2);

  2 の平方根の値 (小数点以下1万桁)
  3 の平方根の値 (小数点以下1万桁)
  5 の平方根の値 (小数点以下1万桁)
  6 の平方根の値 (小数点以下1万桁)
  7 の平方根の値 (小数点以下1万桁)
  8 の平方根の値 (小数点以下1万桁)
  10 の平方根の値 (小数点以下1万桁)

<関数本体>
  squareRoot.c

<説明>
  必要とする桁数以上に、以下の漸化式

      pi+1 = pi2 – 2
      qi+1 = pi qi

  を計算し、最後に、pi+1 / qi+1 で平方根を求める。

  なお、初期値 p0、q0ペル方程式 p02 – M q02= 4 を満たす整数解。ただし、M は計算しようとする平方根の2乗の値。

  数学的帰納法で証明できるが、上の漸化式で得られた pi+1、qi+1 により、

        pi+1 / qi+1 = (M + 4 / (pi qi)2)1 / 2

  となる。したがって平方根の計算において、本アルゴリズムの方法ではニュートン法等ほかの方法よりも効率がよく、数回の計算で必要とする精度のものが得られる。

階乗 n!
     n! = 1 x 2 x 3…x n

をまともにと計算していくと、かなりのコンピュータでも、すぐにオーバフローしてしまう。ということで、本関数のようにそれなりの工夫が必要。

なお、n! の桁数は スターリングの近似式

n! = nne-n(2nπ)1/2

より、

   

n log10(n/e) + (log10(2nπ)/2+1

で求められる。

<関数名>
  factorial —- 階乗 n! の値を多桁数で算出する

<形式>
  int factorial(int *data, int n);

<引数>
  data (出力) n!の値が10進4桁ずつ入った整数の配列
  n   (入力) 階乗 n の値

<関数値>
  配列dataに計算値の入った要素の数(配列の先頭から)

用例
  int data[1000];
  factorial(data, 1000);

  1000! の値

<関数本体>
  factorial.c