パソコンで階乗を計算

現在、fortran90を使って階乗を計算するプログラムを作っております。

プログラム内容は、(n !を求めえるプログラム)
n=0
do i=1,100
n=n*i
enddo

このプログラムを実行すると、１２！までは予想された値が得られるのですが、１３！以降は電卓で計算した値と遙かに異なる値が得られました。
このプログラムは間違っているとは思えないですが、電卓の計算とパソコンの計算が異なる結果になった理由が分かりません。
どなたか、ヒントや参考情報だけでもいいので教えてください。

ちなみにパソコンによる計算結果は、
i n
1 1
2 2
3 6
4 24
5 120
6 720
7 5040
8 40320
9 362880
10 3628800
11 39916800
12 479001600
13 1932053504
14 1278945280
15 2004310016
16 2004189184
17 -288522240
18 -898433024
19 109641728
20 -2102132736
21 -1195114496
22 -522715136
23 862453760
24 -775946240
25 2076180480
26 -1853882368
27 1484783616
28 -1375731712
29 -1241513984
30 1409286144
31 738197504
32 -2147483648
33 -2147483648
34 0
35 0
36 0
36の階乗以降０です。
計算結果が正となるが、結果が違うモノ（例えば、１３！や３１！）は単精度で約１０桁程度しか有効数字が得られないためであると思われるのですが、負になったり、０になる理由が分かりません。

No.5 回答者： Tacosan 回答日時： 2010/11/19 10:14

「単精度でダメ」なときに, なぜ倍精度ではなく整数に走ったのかがかなり謎.

No.4 回答者： AsanoNagi 回答日時： 2010/11/18 22:19

これは、内部的な整数の表現方法の問題です。

まず、この場合、内部的には 32bit 符号付きの整数表現がなされているようです。

たとえば、最小に負の数が出現する
16 2004189184
17 -288522240
を見てみましょう。
これを内部的な16進表示をすると、

16 -> 2004189184 -> 77758000(16) になります。
※16進数で表現すると、８桁分しかないため、こういう数値になっています。

これに、17 を欠けると、16進数では、
7EECD8000 になります。
ここで、実際には、「８桁」しか有効ではありませんから、下位８桁だけが有効になり、
17! は、（7を無視した）EECD8000 になります。
これを、２進数に書き直すと、
EECD8000 = 1110 1110 1010 1011 1000 0000 0000 0000
になって、32bit の一番上位の桁が '1' になっているのがわかります。
ところが、この、「一番上位の桁」は、符号ビットでもあり、ここが、'1' だと、「負の数」だと認識されます。
このため、17! は、負の数だと表示されるわけです。

もうひとつ、この段階でも、下の方の桁は０が並んでいるのに気づいたと思います。

たとえば、普通の10進数で、10を掛ければ、0 が一つ増えます。
10 * 10 = 100, 100 * 10 = 1000
もしも、10進4桁で計算していれば、
1000 * 10 = 10000 ですが、下位４桁しか残りませんから、
1000 * 10 = 0000 = 0 です。

同じように、２進数で、2 を掛ける度に、0がひとつ増えてゆきます。
そこで、32bit だと、2を32回掛けると（たとえば、4 = 2 * 2 なので、４を掛けることは２を２回かけることになるのに注意）

実際には、34 までで偶数が 17個（２を17回）
４の倍数が、34/4 = 8個（つまり、２をさらに、8個掛けることになる）
８の倍数が、34/8 = 4個
16の倍数が、34/16 = 2個
32の倍数が、1 個
存在するため、34! で、2を 17 + 8 + 4 + 2 + 1 = 32個掛けることになるので、34! で、（２進数で表したときに）0 が 32個並ぶことになり、結果的に「0」になります。

それ以降は、０に対するかけ算なので、0のままです。

No.3 回答者： noname#137556 回答日時： 2010/11/18 22:19

C で計算、16進表記してみました。

整数32 ビット、負数は2の補数表現な環境です。
10進表記では質問の数値と一致してます。

01 00000001
02 00000002
03 00000006
04 00000018
05 00000078
06 000002d0
07 000013b0
08 00009d80
09 00058980
0a 00375f00
0b 02611500
0c 1c8cfc00
0d 7328cc00
0e 4c3b2800
0f 77775800
10 77758000
11 eecd8000
12 ca730000
13 06890000
14 82b40000
15 b8c40000
16 e0d80000
17 33680000
18 d1c00000
19 7bc00000
1a 91800000
1b 58800000
1c ae000000
1d b6000000
1e 54000000
1f 2c000000
20 80000000
21 80000000
22 00000000

No.2 回答者： notnot 回答日時： 2010/11/18 22:11

十進数でなく二進数表記で計算するとわかると思います。

32ビット整数で。

最上位ビットは符号ビットです。数値桁のビットから符号ビットに桁あふれを起こしています。

No.1 回答者： kmee 回答日時： 2010/11/18 22:07

コンピュータは有限の桁しかありません。

整数は2進数を使い、2の補数表現という方法で負の値を使うのがよくある方法です。
nビットの整数だと、2の(n-1)乗-1～ -2の(n-1)乗 まで表現できます。
計算結果が有効な桁を越えた場合は、その越えた分の上の桁が無視されます。
例えば、 4bit整数で最大の 0111(=7)を4倍した場合、 11100(=28)となりますが4桁を越える部分が無視され、 1100 となります。
また、この場合、 1100 は8+4=12ではなく、2の補数表現で負の値 -4 となります。

このあたりは電卓や科学分野、コンピュータでの実数型で使われる「有効数字」とは違います。
有効数字では無視するのは下の桁ですし、正負は別に計算するので変化しません。


出てくる数字を見る限り、32bit符号有り整数での計算をしているようです。
使用できる最大値は2の(n-1)乗-1=2147483647 です。
13! = 6706022400 = 17328cc00(16)
であり、すでに33bit目があるので、32bitに削られ
7328cc00(16) = 1932053504(10)
となっています。


また、 33!までは、素因数分解すると、2が31個以下なので2の31乗以内の桁に0で無い値がありますが、34!以降は2が32個以上になるので、最低でも2の32乗になり、32bitの範囲には0しか無いことになります。