【C++】大きな数の階乗計算が途中で止まってしまう

前提・実現したいこと

nを入力として受け取り、n! を 1000000007（10^9 + 7）で割った余りを表示するプログラムを作ろうとしています。nの大きさとしては10万以下を想定しているのですが、nが数万程度以上になると急にプログラムが正常に動作しなくなるという問題が発生しています。原因を突き止めて正常に動作させたいと思っています。

発生している問題・エラーメッセージ

n = 26033 までは正常に動作するのですが、n が26034以上になると急にプログラムが正常に動作しなくなります。具体的には出力まで達することなくプログラムが終了するようになります。

(8/27追記)
皆様ご回答ありがとうございます。お返事が遅くなってしまい大変申し訳ありません。皆様のアドバイスをもとにループを使う方式と、末尾再帰を用いる方式の2パターンで試してみました。末尾再帰を用いる方式ではn が100万を超えるような大きさの問題でも途中で止まることなく計算ができるようです。一方でループを用いる方式では、n が20万程度のところでループがストップしてしまうようです。どちらにせよ今回はn が10万以下の想定なので問題はありませんが、後学のためにこのような差異が表れた原因についてお教えいただけますと幸いです。

最初に作成した末尾再帰のコードは配列に階乗の値を格納することができなかったため、配列に値を格納できるよう末尾再帰を書き直しました。しかしこのコードもまたn が10万強のところでストップしてしまうようです。

大変初歩的なミスでお恥ずかしい限りなのですが、回答者様からご指摘をいただき MAX_N を書き換えていなかったことが原因と分かりました。そのために配列へのアクセスのある関数は途中で止まってしまったようです。MAX_N をもっと大きな値で書き換えた場合、MAX_N 以下の n では末尾再帰、ループともに正常に動作することを確認いたしました。

該当のソースコード

C++
1#pragma GCC target("avx2")
2#pragma GCC optimize("O3")
3#pragma GCC optimize("unroll-loops")
4
5#include <algorithm>
6#include <iostream>
7
8const long long MOD = 1E+09 + 7;
9const int MAX_N = 100000;
10
11using namespace std;
12
13long long modfac[MAX_N + 1];
14long long modfacinv[MAX_N + 1];
15
16long long modPow(long long x, long long a);
17long long modfactor(long long n);
18
19int main()
20{
21    ios::sync_with_stdio(false);
22    cin.tie(nullptr);
23
24    long long n;
25
26    cin >> n;
27    cout << "(" << n << "!) % " << MOD << " = " << modfactor(n) << "\n";
28
29    return 0;
30}
31
32// x^a　をMODで割った余りを求める
33long long modPow(long long x, long long a)
34{
35    if (a == 1)
36        return x % MOD;
37    if (a % 2 == 0) {
38        long long t = modPow(x, a / 2);
39        return (t * t) % MOD;
40    } else {
41        return ((x % MOD) * modPow(x, a - 1)) % MOD;
42    }
43}
44
45
46// n! をMODで割った余りを求める
47long long modfactor(long long n)
48{
49    if (n == 0) {
50        modfacinv[n] = 1;
51        return modfac[n] = 1;
52    }
53    modfac[n] = (n * modfactor(n - 1)) % MOD;
54    modfacinv[n] = (modPow(n, MOD - 2) * modfacinv[n - 1]) % MOD;
55    return modfac[n];
56}
57

試したこと

関数 modfactor の中の modfacinv[n] = (modPow(n, MOD - 2) * modfacinv[n - 1]) % MOD; の記述を消したところ n = 65139 まで正常に計算できるようになりました。このプログラムでは modfacinv は使っていませんが、元々コンビネーションの計算用に作っていたプログラムなのでこちらも計算できるようにしたいと思っています。

###(8/27追記)
末尾再帰を用いる方式

C++
1long long modfactor(long long n, long long f)
2{
3    if (n == 1) {
4        return f;
5    }
6    return modfactor(n - 1, (n * f) % MOD);
7}

末尾再帰＋配列格納(modFactor(n, 0, 1) を呼ぶとn! までの階乗とその逆元が全て求まるようにしました)

C++
1long long modFactor(long long n, long long i, long long f)
2{
3    modfac[i] = f;
4    if (n == i) {
5        modfacinv[n] = modPow(modfac[n], MOD - 2);
6        modFactorInv(n, modfacinv[n]);
7        return f;
8    }
9    return modFactor(n, i + 1, ((i + 1) * f) % MOD);
10}
11
12long long modFactorInv(long long i, long long f)
13{
14    if (i == 0) {
15        modfacinv[i] = f;
16        return f;
17    }
18    modfacinv[i] = f;
19    return modFactorInv(i - 1, (i * f) % MOD);
20}

ループを用いる方式((M - 1)! の逆元がM(M!)^-1 で求められることを用いています)

C++
1// n!を求める
2long long modfactor(long long n)
3{
4    modfac[0] = 1;
5    long long p = 1;
6    while (p <= n) {
7        //cout << p << "\n";
8        modfac[p] = (p * modfac[p - 1]) % MOD;
9        p++;
10    }
11    modfacinv[n] = modPow(modfac[n], MOD - 2);
12    modfactorInv(n);
13    return modfac[n];
14}
15
16// (n!)^-1 を求める
17void modfactorInv(long long n)
18{
19    long long p = n - 1;
20    while (p >= 0) {
21        modfacinv[p] = ((p + 1) * modfacinv[p + 1]) % MOD;
22        p--;
23    }
24}

行動規範の内容に同意します

回答3件

スタックオーバーフローしているということはありませんか？

再帰呼び出しを繰り返していくとスタックを消費しますので、再帰できる回数には限度があります。

投稿2020/08/21 14:04

maisumakun

総合スコア145930

maisumakun

2020/08/21 14:05

あとは、最適化系のプラグマが書かれていますが、これらを止めてみることも有効かもしれません。

行動規範の内容に同意します

ベストアンサー

maisumakunさんも指摘されているように、再帰呼び出しによるスタックオーバーフローでしょう。
再帰呼び出しをループで書き直せば動くはずです。
関数呼び出しが減るので高速化にもなります。

C++
1// x^a　をMODで割った余りを求める
2long long modPow(long long x, long long a)
3{
4    long long ret = 1;
5    x %= MOD;
6    while(0 < a)
7    {
8        if (a % 2 == 0) {
9            x = (x * x) % MOD;
10            a /= 2;
11        }
12        else {
13            ret = (ret * x) % MOD;
14            --a;
15        }
16    }
17    return ret;
18}
19
20// n! をMODで割った余りを求める
21long long modfactor(long long n)
22{
23    modfacinv[0] = 1;
24    modfac[0] = 1;
25
26    for (long long cnt = 1; cnt <= n; cnt++)
27    {
28        modfac[cnt] = (cnt * modfac[cnt - 1]) % MOD;
29        modfacinv[cnt] = (modPow(cnt, MOD - 2) * modfacinv[cnt - 1]) % MOD;
30    }
31    return modfac[n];
32}

投稿2020/08/22 00:38

編集2020/08/22 01:35

SHOMI

総合スコア4079

shukrin

2020/08/27 10:35

ご回答ありがとうございます。お返事が遅くなってしまい大変申し訳ありません。アドバイスをもとにループを用いた階乗計算コードを作成したのですが（コードを追記いたしましたのでご覧ください）ループ回数が20万程度になると止まってしまうようです。この原因も再帰計算と同じくスタックオーバーフローによるものなのでしょうか。

SHOMI

2020/08/27 12:07 編集

ループでスタックは消費しません。私の貼ったコードでMAX_Nを200000以上にすれば問題なく実行できます。

shukrin

2020/08/27 12:30 編集

MAX_N を書き換えるのを忘れていました。書き換えたところ問題なく動作しました。大変失礼いたしました。

行動規範の内容に同意します

ご指摘がありましたの回答を削除しました。

投稿2020/08/21 19:03

編集2020/08/21 21:17

mjk

総合スコア303

raccy

2020/08/21 20:35

質問の課題は1000000007で割った余りを求めるというものであるため、剰余をきちんと扱えればlong longで計算しても、最終結果や計算途中でLLONG_MAXを越えることはありません。回答者のコードについて、その部分は問題ないように思えます。質問の重要な部分を見逃して、勘違いな回答をしているため、低評価させていただきます。

mjk

2020/08/21 20:41

ご指摘ありがとうございます。再度確認してみます。

raccy

2020/08/21 20:50

回答にある再帰処理のコードですが、末尾再帰になっていますので、一般的なコンパイラ(GCCやcl.exe等)で標準的な最適化を行っていれば、末尾再帰最適化によってスタックを消費しない形でコンパイルされます。通常はスタックオーバーフローが発生しないような場合を例示しているのはいささかおかしいと思います。なお、デバッグモードでは最適化されないので、デバッグを有効にした状態で確認することはできません。また、質問者のコードは末尾再帰ではありません。

mjk

2020/08/21 20:59

すみません知識が及ばずraccyさんのコメントの単語（末尾再帰など）すら理解出来ないので恐縮です。身の程もわきまえず知識ないのに回答してごめんなさい。後ほど削除しておきます。コメント欄で恐縮なのですが後学の為に１つだけ教えて下さい。質問のコードの関数にあるmodfac[n] = (n * modfactor(n - 1)) % MOD; 階乗計算をしている過程で20桁を超えているはずなのですが内部的には何桁になろうが正確な計算が出来ているということでよろしいですか？自分の浅はかな考えでは計算途中に%MODする前に20桁超えているので計算結果がおかしくなっていると思っていました。

raccy

2020/08/21 21:04

nは10万以下想定なので6桁以下、modfactor()が返す値は常に10桁以下、なので、n * modfactor(n - 1)は16桁以下です。

mjk

2020/08/21 21:08

何度もすみません。途中でn!階乗を求めていますが21!などの20桁超えてる階乗の計算は20桁を内部的に超えても大丈夫ということでしょうか？

mjk

2020/08/21 21:17 編集

＞＞n! を 1000000007（10^9 + 7）で割った余りを表示するプログラムごめんなさい。質問にもこう書いてあったので勘違いしたのかもしれません。もしかするとコード内ではnの階乗してないということですかね。

raccy

2020/08/21 21:19

n!階乗を1000000007で割った余りを求めています。階乗の結果そのものを求めているのではありません。余りのみを求める場合は、桁あふれが起きないようにする方法が存在します。競技プログラミング等において基本のテクニック(アルゴリズム)の一つになっていますので、詳しいことはアルゴリズムや競技プログラミングの本やサイト等をを参考にしてください。

mjk

2020/08/21 21:21

＞＞余りのみを求める場合は、桁あふれが起きないようにする方法が存在します。貴重な勉強になりました！早速調べてみます。またコメント欄の質問にも関わらず回答して頂きありがとうございました！

行動規範の内容に同意します

あなたの回答