ネイピア数eの小数点以下の桁数をできるだけ多く求めて表示する
小数点以下の桁数が例えば500桁まで求めるこができること
多倍長計算として比較的簡単に実現できます。私の手元では数十行程度のプログラムで小数点以下1000桁を求めることができました。
1/n!のn!が配列で表現されているとき、
1/n!をどのようにして計算するかを教えていただきたいです。
具体的にはn! = A[3]A[2]A[1]であるとき、
1/n! = 1 / A[3]A[2]A[1] の計算の仕方を知りたいです。
「n! = A[3]A[2]A[1]である」の意味がわかりませんが、それはともかく・・・
ネイピア数
階乗を使うのは、このページの「微分積分学の基本的な関数を使った定義」式を使おうということだと思います。その式をそのままプログラムすれば次のようなコードになるでしょう。
C
1 double napier = 1.0;
2 for (int n = 1; 十分な回数; n++) {
3 napier += 1.0 / factrial(n);
4 }
5 printf("napier = %g\n", napier);
factrial(i) は i の階乗を計算する関数を想定したものです。しかしプログラムを作るには、階乗そのものを計算しないほうが良いと考えます。
上のプログラムは次のように書き変えることができます。
C
1#include <stdio.h>
2#include <float.h>
3int main(void)
4{
5 double napier = 1.0, differ = 1.0;
6 for (int n = 1; differ > DBL_EPSILON; n++) {
7 differ /= n;
8 napier += differ;
9 }
10 printf("napier = %.20g\n", napier);
数学的な意味は同じですが、こう書き変えれば
- 直接、階乗を計算しないですむ
- 多倍長計算にも向いている
ことから、計算機的に良い方法だと考えています。まずは、
- このコードはどう書き変えたのか
- 実際に動作結果を見る(桁数の少ない近似値が表示されるはず)
をご確認ください。その上で、この計算方法をお伝えすべきかどうかを確認する必要があると思っています。なぜならお気づきと思いますが、要件と思われる「階乗の配列」を使わない方法だからです。
なお、階乗の配列を使い500桁計算する方法であれば、私の手には負えませんので、他の方におたずねください。
気づいたら反応が無いまま日が経ちましたが、せっかくなのでコードをアップします。
この問題は多倍長演算としては初歩的な部類に入ると思います。私でもゼロから書けた位。
多倍長演算はある数を配列で表すことで多くの桁数を保持しますが、アルゴリズム等によって配列の使い方はいくつかあるようです。
今回は、配列の一要素に0〜99999の値を、即ち5桁分を格納することにしました。たとえば int napier[200]; とした場合、10 進数で最大 1000 桁を保持します(ただし、最下位2桁位は誤差が残るので、小数点以下 1000 桁を求める場合 1005 桁まで求める)。
こうした理由はふたつあります。
- 表示が簡単、"%05d"で済む。上記 wikipedia に挙げられた、5桁ずつの表示に合わせるのに都合がよい。
- 私のやり方なら小さい整数で割り算ができればよく、筆算の要領で簡単にプログラムできる。
今回作った足し算と割り算を解説します。どちらもそんなに難しくありません。
add() 関数は多倍長 + 多倍長の足し算です。2つの数の、下の桁から、同じ桁同士を足していきます。100000 を超えた場合、99999 以下を残し、次の桁に桁上がり(carry)を送る・・・筆算のやり方です。
div() 関数は多倍長 ÷ 整数の割り算です。今度は上の桁から整数で割り、余りを次の桁への繰り下がり(borrow)として割り算を繰り返します。これも筆算の応用です。今回は割る数が比較的小さい整数に収まるので、簡単なコードで済みます。
differ[] 配列を除数 divisor で割った結果、配列全てが0になった時点で計算は終了です。小数点以下 500 桁を求める場合は 256 で割った時、1000 桁を求める場合は、452 で割った時、それぞれ終了しました。
このことは 256 の階乗も 500 桁程度になることを意味します。質問の
1/n! = 1 / A[3]A[2]A[1] の計算の仕方を知りたい
は多倍長 ÷ 多倍長の割り算、即ち多倍長整数での割り算を求めているようにも受け取れます。しかし、階乗を求めるのも、割り算するのも、大変です。計算の負荷をムダに重くするばかりで、もはや挑戦しようという気になりません。
C
1#include <stdio.h>
2#include <stdbool.h>
3
4#ifdef DEBUG
5# define DIGITS (40) // with -D DEBUG
6#else
7# define DIGITS (1000) // 小数点以下に桁数--
8#endif
9#define ARRSIZE ((DIGITS / 5) + 2) // 必要な配列サイズ
10#define UNIT (100000) // 配列の一要素は00000〜99999を格納
11
12void display(char *msg, int xx[]); // 配列 xx[] を表示
13void add(int sum[], int yy[]); // sum += yy
14bool div(int ddd[], int divisor); // ddd /= divisor
15
16int main(void)
17{
18 int napier[ARRSIZE] = {1}; // napier = 1.0000...
19 int differ[ARRSIZE] = {1}; // differ = 1.0000...
20
21 for (int n = 1; ; n++) {
22 bool zero = div(differ, n); // differ /= n
23 if (zero) {
24 printf("DEBUG: exit on %d as divisor.\n", n);
25 break; // 商が0なら終了
26 }
27 add(napier, differ); // napier += differ
28
29#ifdef DEBUG // 途中経過を表示してみる
30 printf("DEBUG: divided by %d\n", n);
31 display("DEBUG: differ", differ);
32 display("DEBUG: napier", napier);
33#endif
34 }
35
36 display(" napier", napier); // 最終結果
37 return 0;
38}
39
40void display(char *msg, int xx[]) // 配列 x[] を表示
41{
42 printf("%s = %d.", msg, xx[0]); // 整数部と小数点
43 for (int k = 1; k < ARRSIZE; k++)
44 printf("%05d ", xx[k]); // 5桁ずつ表示
45 printf("\n");
46}
47
48void add(int sum[], int yy[]) // sum += y
49{
50 int carry = 0; // 桁上がり
51 for (int k = ARRSIZE - 1; k >= 0; --k) {
52 // 同じ桁位置の2数+下位からの繰上がり
53 int num = sum[k] + yy[k] + carry;
54
55 if (num >= UNIT) { // 桁あふれ?
56 num -= UNIT; // 99999 以下を残す
57 carry = 1; // 桁上がりあり
58 } else {
59 carry = 0; // 桁上がりなし
60 }
61 sum[k] = num; // その桁の足し算結果
62 }
63}
64
65bool div(int ddd[], int divisor) // ddd /= divisor
66{
67 bool allzero = true; // 商の全ての桁が0か?
68 int borrow = 0; // 次の位置に送る余り=桁下がり
69
70 for (int k = 0; k < ARRSIZE; k++) {
71 // 「上の桁から桁下がり * 100000 + この桁位置の値」を割る
72 int num = borrow * UNIT + ddd[k];
73
74 borrow = num % divisor; // 余りを次の位置に送る
75 ddd[k] = num / divisor; // その位置の商
76 if (ddd[k] > 0) allzero = false; // 結果は0じゃない
77 }
78 return allzero; // 商が0になったら true を返す
79}
80