c言語の累乗の計算について

Question

初心者です。暗号化のプログラムを作っています。その中で、mのe乗を求めたいのですが正確な値が出ません。それよりも前の部分は正常に動いています。/暗号化/という部分です。どこを直せば良いのでしょうか。

 #include<stdio.h>
#include <stdlib.h> // rand, srand関数
#include <time.h>   // time関数
#include <math.h>

int gcd(int x, int y);  // 最大公約数を求める関数





int main(void){

    unsigned long long int p=0,q=0;
    int i,l;//繰り返し用変数
    int flag1=0,flag2=0;
    int count=0;
    unsigned long long int n;
    unsigned long long int o;//オイラー関数
    unsigned long long int E;
    unsigned long long int num=0;
    int x,y,r;//ユーグリッドの互除法で使う
    int gcd_num=0;//最大公約数
    unsigned long long int d=0;
    unsigned long long int judge , k=0;//dを求める時に必要,dはあまり
    unsigned long long int c;
    int m=0;//平文
    unsigned long long int me=0;//mのe乗
    

    srand((unsigned int)time(NULL)); // 現在時刻の情報で初期化
    
    /*pの値を求める*/
    while(p<1000 || flag1==1 || p==0){//pの条件を満たしてない限り繰り返す
        flag1=0;
        p=rand()%10000;//以降の計算に時間がかかるため4桁に指定
        for( l=2;l<p;++l ) {
            if( p%l==0 ) {//割りきれる場合
            flag1 = 1;
            break;//繰り返し終了
            }
        }
    }
    printf("p:%llu\n",p);//pの値出力

    /*qの値を求める*/
    while(q<1000 || flag2==1 || q==0){//qの条件を満たしてない限り繰り返す
        flag2=0;
        q=rand()%10000;//以降の計算に時間がかかるため4桁に指定
        for( i=2;i<q;++i ) {
            if( q%i==0 ) {//割りきれる場合
                flag2 = 1;
                break;//繰り返し終了
            }
        }
    }
    printf("q:%llu\n",q);//qの値出力

    n=p*q;
    
    o=(p-1)*(q-1);//オイラー関数
    printf("φ:%llu\n",o);

    /*eの値を求める*/   
    while(gcd_num!=1 || E==0 ||E==1){
        E=rand()%100;
        
        if(o<=E){//E,oのうち大きい方をx小さい方をyに設定
            x=E;
            y=o;
        }
        else{
            y=E;
            x=o;
        }
        gcd_num = gcd(x, y);
     }
    printf( "e:%llu\n", E );
    
    
    /*dの値を求める*/ 
    while(judge!=0 || d>o ){
        d=rand()%99960004;//φの範囲:0~9998*9998
        judge=(E*d-1)%o;
    }
    printf( "d:%llu\n", d );
    k=(E*d-1)/o;

    printf("k:");
    if(E*d==k*o+1){//ed=φk+1か確認
        printf("%llu\n", k );
    }
    
    /*暗号化*/
    while(m==0){
        m=rand()%10;//mは2桁
    }
    printf("m:%d\n",m);
    me=m*m;
    for(i=1;i<E-1;i++){
        me=me*m;
    }
    printf("m^e:%llu\n", me );
    c=me%n;
    printf("c:%llu\n",c);
    
    
    
    
    return 0;
}


// 最大公約数を求める関数
int gcd(int x, int y){
    int r;
    

    // ユーグリッドの互除法
    while((r = x % y) != 0) { // 割り切れるまでループ
        x = y;
        y = r;
    }
    return y;
}


    

Comments

[y_waiwai]

どういう値が出るんでしょうか

[退会済みユーザー]

例えば3^7の場合、67,585,198,634,817,523,235,520,443,624,317,923
ですが、このプログラムでは6,572,658,003,009,201,123と出力されます。

[y_waiwai]

階乗じゃないんですか？
それはなんの数字でしょうか

[raccy]

67,585,198,634,817,523,235,520,443,624,317,923は3^73ですよね？

[otn]

3の7乗 は 2187 では？

[退会済みユーザー]

そうです！すみません、間違えました3^73です。

Answer 1

単純に桁溢れだと思います。ちなみに、これは他の部分を読む限り、RSA暗号あたりの暗号/復号のものでしょう。だとすると累乗計算はO(log N)の形にするのが適切です。以下、動作未確認。

//巨大整数を扱う型をTとします。
T power(T x, T n)
{
  T temp;
  if(0 == n){//0乗はxに関わらず1。
    return 1;
  }else if(1 == n){//1乗はxそのもの。
    return x;
  }
  temp = power(x, n / 2);//再帰でx ^ (n / 2)を求めて一時変数へ。
  return temp * temp * (n % 2 ? x : 1);//x ^ nを計算結果とする。
}

このプログラムなら、nを半分しながら反復していくので計算量が対数オーダになります。

Answer 2

一般的なx86_64を含むほとんどの環境においてunsigned long long intの最大値は18446744073709551615です。mとeの値によってはこの最大値を超えてしまうため、その場合は正確な値がでることはありません。最大であり得る9の99乗の値(29512665430652752148753480226197736314359272517043832886063884637676943433478020332709411004889)が十分に入るような整数の型を考えるか、GNU MPのような多倍長整数ライブラリを使う必要があります。

Comments

[退会済みユーザー]

回答ありがとうございます。
GNU MPというものを初めて聞いたのですが、調べてもいまいち理解できません。どのように使うのでしょうか。

[raccy]

どのように使うのかと言われてもマニュアル https://gmplib.org/manual/ を読んでくださいとしか言いようがないです。

Answer 3

細かい状況は分かりませんが、とりあえずバグとして考えられるのは

1. eの値を求める際に低確率でgcdで0剰余演算が発生する(質問とは関係なさそうですが)
2. eの値が低確率で1になるがmeとしてm^2が出力される
3. meが割とオーバーフローしやすい

Eの詳しい算出方法は十分追い切れておりませんが恐らく3.が主な原因かと思います。
unsigned long long intの最大値が2^64-1として、例えばmが9の時eが20もあれば溢れます。