回答編集履歴
2
コードを追加
answer
CHANGED
|
@@ -14,8 +14,8 @@
|
|
|
14
14
|
階乗を使うのは、このページの「微分積分学の基本的な関数を使った定義」式を使おうということだと思います。その式をそのままプログラムすれば次のようなコードになるでしょう。
|
|
15
15
|
```C
|
|
16
16
|
double napier = 1.0;
|
|
17
|
-
for (int
|
|
17
|
+
for (int n = 1; 十分な回数; n++) {
|
|
18
|
-
napier += 1.0 / factrial(
|
|
18
|
+
napier += 1.0 / factrial(n);
|
|
19
19
|
}
|
|
20
20
|
printf("napier = %g\n", napier);
|
|
21
21
|
```
|
|
@@ -27,8 +27,8 @@
|
|
|
27
27
|
int main(void)
|
|
28
28
|
{
|
|
29
29
|
double napier = 1.0, differ = 1.0;
|
|
30
|
-
for (int
|
|
30
|
+
for (int n = 1; differ > DBL_EPSILON; n++) {
|
|
31
|
-
differ /=
|
|
31
|
+
differ /= n;
|
|
32
32
|
napier += differ;
|
|
33
33
|
}
|
|
34
34
|
printf("napier = %.20g\n", napier);
|
|
@@ -50,7 +50,7 @@
|
|
|
50
50
|
|
|
51
51
|
この問題は多倍長演算としては初歩的な部類に入ると思います。私でもゼロから書けた位。
|
|
52
52
|
多倍長演算はある数を配列で表すことで多くの桁数を保持しますが、アルゴリズム等によって配列の使い方はいくつかあるようです。
|
|
53
|
-
今回は、配列の一要素に0〜99999の値を、即ち５桁分を格納することにしました。たとえば int napier[200]; とした場合、10 進数で最大 1000 桁を保持します（ただし、最下位２桁位は誤差が残るので、小数点以下 1000 桁を求める場合 1005 桁まで求める）。
|
|
53
|
+
今回は、**配列の一要素に0〜99999の値を、即ち５桁分を格納**することにしました。たとえば int napier[200]; とした場合、10 進数で最大 1000 桁を保持します（ただし、最下位２桁位は誤差が残るので、小数点以下 1000 桁を求める場合 1005 桁まで求める）。
|
|
54
54
|
|
|
55
55
|
こうした理由はふたつあります。
|
|
56
56
|
- 表示が簡単、"%05d"で済む。上記 wikipedia に挙げられた、5桁ずつの表示に合わせるのに都合がよい。
|
|
@@ -67,4 +67,87 @@
|
|
|
67
67
|
|
|
68
68
|
> 1/n! = 1 / A[3]A[2]A[1] の計算の仕方を知りたい
|
|
69
69
|
|
|
70
|
-
は多倍長 ÷ 多倍長の割り算、即ち多倍長整数での割り算を求めているようにも受け取れます。しかし、階乗を求めるのも、割り算するのも、大変です。計算の負荷をムダに重くするばかりで、もはや挑戦しようという気になりません。
|
|
70
|
+
は多倍長 ÷ 多倍長の割り算、即ち多倍長整数での割り算を求めているようにも受け取れます。しかし、階乗を求めるのも、割り算するのも、大変です。計算の負荷をムダに重くするばかりで、もはや挑戦しようという気になりません。
|
|
71
|
+
|
|
72
|
+
```C
|
|
73
|
+
#include <stdio.h>
|
|
74
|
+
#include <stdbool.h>
|
|
75
|
+
|
|
76
|
+
#ifdef DEBUG
|
|
77
|
+
# define DIGITS (40) // with -D DEBUG
|
|
78
|
+
#else
|
|
79
|
+
# define DIGITS (1000) // 小数点以下に桁数--
|
|
80
|
+
#endif
|
|
81
|
+
#define ARRSIZE ((DIGITS / 5) + 2) // 必要な配列サイズ
|
|
82
|
+
#define UNIT (100000) // 配列の一要素は00000〜99999を格納
|
|
83
|
+
|
|
84
|
+
void display(char *msg, int xx[]); // 配列 xx[] を表示
|
|
85
|
+
void add(int sum[], int yy[]); // sum += yy
|
|
86
|
+
bool div(int ddd[], int divisor); // ddd /= divisor
|
|
87
|
+
|
|
88
|
+
int main(void)
|
|
89
|
+
{
|
|
90
|
+
int napier[ARRSIZE] = {1}; // napier = 1.0000...
|
|
91
|
+
int differ[ARRSIZE] = {1}; // differ = 1.0000...
|
|
92
|
+
|
|
93
|
+
for (int n = 1; ; n++) {
|
|
94
|
+
bool zero = div(differ, n); // differ /= n
|
|
95
|
+
if (zero) {
|
|
96
|
+
printf("DEBUG: exit on %d as divisor.\n", n);
|
|
97
|
+
break; // 商が０なら終了
|
|
98
|
+
}
|
|
99
|
+
add(napier, differ); // napier += differ
|
|
100
|
+
|
|
101
|
+
#ifdef DEBUG // 途中経過を表示してみる
|
|
102
|
+
printf("DEBUG: divided by %d\n", n);
|
|
103
|
+
display("DEBUG: differ", differ);
|
|
104
|
+
display("DEBUG: napier", napier);
|
|
105
|
+
#endif
|
|
106
|
+
}
|
|
107
|
+
|
|
108
|
+
display(" napier", napier); // 最終結果
|
|
109
|
+
return 0;
|
|
110
|
+
}
|
|
111
|
+
|
|
112
|
+
void display(char *msg, int xx[]) // 配列 x[] を表示
|
|
113
|
+
{
|
|
114
|
+
printf("%s = %d.", msg, xx[0]); // 整数部と小数点
|
|
115
|
+
for (int k = 1; k < ARRSIZE; k++)
|
|
116
|
+
printf("%05d ", xx[k]); // 5桁ずつ表示
|
|
117
|
+
printf("\n");
|
|
118
|
+
}
|
|
119
|
+
|
|
120
|
+
void add(int sum[], int yy[]) // sum += y
|
|
121
|
+
{
|
|
122
|
+
int carry = 0; // 桁上がり
|
|
123
|
+
for (int k = ARRSIZE - 1; k >= 0; --k) {
|
|
124
|
+
// 同じ桁位置の２数＋下位からの繰上がり
|
|
125
|
+
int num = sum[k] + yy[k] + carry;
|
|
126
|
+
|
|
127
|
+
if (num >= UNIT) { // 桁あふれ？
|
|
128
|
+
num -= UNIT; // 99999 以下を残す
|
|
129
|
+
carry = 1; // 桁上がりあり
|
|
130
|
+
} else {
|
|
131
|
+
carry = 0; // 桁上がりなし
|
|
132
|
+
}
|
|
133
|
+
sum[k] = num; // その桁の足し算結果
|
|
134
|
+
}
|
|
135
|
+
}
|
|
136
|
+
|
|
137
|
+
bool div(int ddd[], int divisor) // ddd /= divisor
|
|
138
|
+
{
|
|
139
|
+
bool allzero = true; // 商の全ての桁が０か？
|
|
140
|
+
int borrow = 0; // 次の位置に送る余り=桁下がり
|
|
141
|
+
|
|
142
|
+
for (int k = 0; k < ARRSIZE; k++) {
|
|
143
|
+
// 「上の桁から桁下がり * 100000 + この桁位置の値」を割る
|
|
144
|
+
int num = borrow * UNIT + ddd[k];
|
|
145
|
+
|
|
146
|
+
borrow = num % divisor; // 余りを次の位置に送る
|
|
147
|
+
ddd[k] = num / divisor; // その位置の商
|
|
148
|
+
if (ddd[k] > 0) allzero = false; // 結果は０じゃない
|
|
149
|
+
}
|
|
150
|
+
return allzero; // 商が０になったら true を返す
|
|
151
|
+
}
|
|
152
|
+
|
|
153
|
+
```
|
1
コードと解説を追加
answer
CHANGED
|
@@ -42,4 +42,29 @@
|
|
|
42
42
|
- 実際に動作結果を見る（桁数の少ない近似値が表示されるはず）
|
|
43
43
|
|
|
44
44
|
をご確認ください。その上で、この計算方法をお伝えすべきかどうかを確認する必要があると思っています。なぜならお気づきと思いますが、要件と思われる「階乗の配列」を使わない方法だからです。
|
|
45
|
-
なお、階乗の配列を使い500桁計算する方法であれば、私の手には負えませんので、他の方におたずねください。
|
|
45
|
+
なお、階乗の配列を使い500桁計算する方法であれば、私の手には負えませんので、他の方におたずねください。
|
|
46
|
+
|
|
47
|
+
----
|
|
48
|
+
|
|
49
|
+
気づいたら反応が無いまま日が経ちましたが、せっかくなのでコードをアップします。
|
|
50
|
+
|
|
51
|
+
この問題は多倍長演算としては初歩的な部類に入ると思います。私でもゼロから書けた位。
|
|
52
|
+
多倍長演算はある数を配列で表すことで多くの桁数を保持しますが、アルゴリズム等によって配列の使い方はいくつかあるようです。
|
|
53
|
+
今回は、配列の一要素に0〜99999の値を、即ち５桁分を格納することにしました。たとえば int napier[200]; とした場合、10 進数で最大 1000 桁を保持します（ただし、最下位２桁位は誤差が残るので、小数点以下 1000 桁を求める場合 1005 桁まで求める）。
|
|
54
|
+
|
|
55
|
+
こうした理由はふたつあります。
|
|
56
|
+
- 表示が簡単、"%05d"で済む。上記 wikipedia に挙げられた、5桁ずつの表示に合わせるのに都合がよい。
|
|
57
|
+
- 私のやり方なら小さい整数で割り算ができればよく、筆算の要領で簡単にプログラムできる。
|
|
58
|
+
|
|
59
|
+
今回作った足し算と割り算を解説します。どちらもそんなに難しくありません。
|
|
60
|
+
|
|
61
|
+
add() 関数は**多倍長 + 多倍長の足し算**です。２つの数の、下の桁から、同じ桁同士を足していきます。100000 を超えた場合、99999 以下を残し、次の桁に桁上がり（carry）を送る・・・筆算のやり方です。
|
|
62
|
+
|
|
63
|
+
div() 関数は**多倍長 ÷ 整数の割り算**です。今度は上の桁から整数で割り、余りを次の桁への繰り下がり（borrow）として割り算を繰り返します。これも筆算の応用です。今回は割る数が比較的小さい整数に収まるので、簡単なコードで済みます。
|
|
64
|
+
|
|
65
|
+
differ[] 配列を除数 divisor で割った結果、配列全てが０になった時点で計算は終了です。小数点以下 500 桁を求める場合は 256 で割った時、1000 桁を求める場合は、452 で割った時、それぞれ終了しました。
|
|
66
|
+
このことは 256 の階乗も 500 桁程度になることを意味します。質問の
|
|
67
|
+
|
|
68
|
+
> 1/n! = 1 / A[3]A[2]A[1] の計算の仕方を知りたい
|
|
69
|
+
|
|
70
|
+
は多倍長 ÷ 多倍長の割り算、即ち多倍長整数での割り算を求めているようにも受け取れます。しかし、階乗を求めるのも、割り算するのも、大変です。計算の負荷をムダに重くするばかりで、もはや挑戦しようという気になりません。
|