標準入力から与えらる正整数 m を読み取り、m 番目のメルセンヌ数を書き出すプログラムを組みたいです

Question

前提・実現したいこと

標準入力から与えらる正整数 m を読み取り、m 番目のメルセンヌ数を書き出すプログラムを組みたいです。
※メルセンヌ数とは、 2^m-1 です。

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

標準入力に100を入れ実行すると、1267650600228229401496703205376　になります。
正解は1267650600228229401496703205375　なのですが・・・
54以降はなぜか、1足されたものが出力されます。

該当のソースコード

#include <stdio.h>
#include <math.h>

int main(int argc, char *argv[]) {
    double m;
    double a;
    double b;
    
    scanf("%lf", &m);
    
    a = pow(2,m);
    b = (a-1);
    
    printf("%.0f", b);
    printf("\n");
    return 0;
}

Answer 1

doubleは浮動小数点数型です。

浮動小数点数というのはざっくり言うと
(+/-)1.aaaa x 2^bbbb
というような数値表現です。例えば1.01001 x 2^5のように。

そして、doubleの場合、aaaaの部分(仮数部)は52ビット、2進数で52桁と決まっています。
1.の部分も合わせると、53桁になります。

(2^53)-1 は2進数で1111111..11(1が合計53個)ですから、
doubleで誤差なく表すことができるメルセンヌ数はm=53までとなります。

Comments

[Taka787]

m=53以上の値を求めるにはどのようにすればよいのでしょうか？

[nob.]

そういうのを「多倍長整数」といいます。
多倍長整数を扱うライブラリを探してください。
無ければ自分で作りましょう。

[Taka787]

皆さん分かりやすい回答ありがとうございます。
解決できました。

Answer 2

こんにちは

まず初めに

他の方もおっしゃるとおり

倍精度浮動小数点数 任意精度演算

で実装するのがいいと思います。
質問欄にあるソースコードを見て、質問者もその方針で進めていると思います。

しかし、他の方法を思いついたのでその方法を説明したいと思います。
あまり参考にならないかもしれませんが、
一応、求めたいものは求まったので少しでも参考になれば幸いです。

どんな方法で実装するのか?

そもそもintやlongではなくdoubleで求めたかを考えるとそれは
結果がとても大きな自然数であり、intやlongではオーバーフローしてしまうからです。
そこで考えたのは各桁をそれぞれ配列に格納すればいいと思いました。
各桁の数字は0以上9以下なのでintで十分に収まります。

(i)初期化として、配列の0番目を1それ以外は0とする。
(ii)配列に格納されている数字をそれぞれ2倍する。
(iii)配列に格納されている数字が10を超えた場合、次の桁へ繰り上げる。
(iv) (ii),(iii)をそれぞれn回繰り返す。
(v) 一桁目(配列の0番目)を1引く.

以下、私が書いたソースコードを載せておきます。

注意

(i)今回は100番目のメルセンヌ数を求めたいので、配列の要素数は31くらいでいいですが、
もっと大きなメルセンヌ数を求めたい場合を考え,配列の要素数を100にしています。

(ii) n桁目の数字は配列のn - 1番目に格納されていますが、表示するときは、逆順にする必要があります。

例えば、1桁目の数字は配列の0番目に格納されていますが、1桁目は数字の一番右にあります。
しかし、for文で回すときカウントする数字を0から回すと1桁目が1一番左に表示されます。
そのため、、for文で回すときカウントする数字を(桁数 - 1)から回す必要があります。

#include <stdio.h>
#include <stdlib.h>

#define BUFSIZE 256
#define MAX_NUMBER 100

void Mersenne_primes(int);

int main(void) {
    
    char buf[BUFSIZE];
    int n;
    
    fgets(buf,sizeof(buf),stdin);
    n = atoi(buf);
    Mersenne_primes(n);
    
    return 0;
}

void Mersenne_primes(int n)
{
    int answer[MAX_NUMBER];
    int carry = 0; //桁上がり
    int digit = 1; //現在の桁数
    int i,j;
    int num;
    for(i = 0; i < MAX_NUMBER; i++){
        answer[i] = 0;
    }
    answer[0] = 1;
    
    for(i = 0; i < n; i++){
        for(j = 0; j <= digit; j++){
            num = answer[j] * 2 + carry;
            carry = num / 10;
            answer[j] = num % 10;
            if(carry > 0 && j + 1 > digit){
                digit++;
            }
        }
        carry = 0;
    }
    answer[0]--;
 
    for(i = digit; i >= 0; i--){
        printf("%d",answer[i]);
    }
}

<結果>

Comments

[Taka787]

詳しい回答ありがとうございました。
おかげさまで理解できました！

[退会済みユーザー]

質問内容とは全く別の方針だったのですが、役に立ったので良かったです。
ちなみに私は最初,質問のタイトルを見てメルセンヌ素数と勘違いをしました。(笑)

[退会済みユーザー]

ちなみに計算したところ,配列の要素数が100(つまり100桁)で求めることができるの自然数の最大値は2^332 - 1です。

[Taka787]

勘違いさせてすみません( ;∀;)
有難うございました。( ´∀｀ )

[退会済みユーザー]

いえいえ、これは、私の勘違いなので....(おそらく数学のメルセンヌ素数と完全数が好きなのが原因)

Answer 3

double ではなくて long などで試せばおかしな値になってしまうのでわかりやすいのですが、コレは double 計算精度を超える値を使っていることによる丸め誤差ですね。
たとえば、 a が巨大な数になった時、大きな値を表現するためのビットに対して -1 などの小さな値は決められた数のビット数では表現できません。

a = 10000000000000000;
b = a - 1;

もし具体的にこれらの値を計算したい場合は、「任意精度演算」などを参照すると良いと思います。

Comments

[Taka787]

初めて知りました。とても参考になりました。ありがとうございます。

Answer 4

倍精度浮動小数点数

これをみて、２の１００乗という数値がどういう表現になるか考えてみましょう。
そして、そこから１を引く、ということがどういうことかを考えると、おのずと理解できるかと。

Comments

[Taka787]

なるほどです。よく理解できました。有難うございました。