GNU MP(GMP)を用いたプログラミングについて

理動 · #1

先日因数分解のプログラムに関する質問をさせていただきました。
その際GMPを使うと大きな桁数も扱えるとアドバイスをいただきましたので、
それを用いてソースコードを書き直しました。

今回はポラードのρ法を用いてレピュニットの因数分解をしようとしています。
レピュニットとはすべての桁が1である自然数のことです。
今回は1が53個並んだR(53)を因数分解するプログラムを組もうとしています。

ここまで正しく動作したとコメントアウトしてある箇所の後のループから抜け出せなくなっています。
原因は該当ループの2回目の V の値が期待するものと違うため、
そのあとの最大公約数が正しく求められず正常な動作をしないと考えています。

何かおかしな操作をしているか、足りていないものはあるか、
自力で見つけられなかったのでこちらに質問させていただきました。
ソースコードがかなり読みにくいとは思いますが、どうかお力添えください。

OS:Win7
コンパイラ:gcc
ライブラリ:GMP(The GNU MP Bidnum Library)
今回はCygwinを使用しています。

出力結果は
11111111111111111111111111111111111111111111111111111
= 107 * 1659431 * 1325815267337711173 * 47198858799491425660200071
となるのが正しいです。

以下にアルゴリズムとソースコードを示します。

nはユーザが入力する大きな桁数の整数です。
1[乱数種の選択]
　ランダムにa∈[1,n-3]を選ぶ;
　ランダムにs∈[0,n-1]を選ぶ;
　U=V=s;
　F(x)=(x * x + a) mod nと定める;
2[因子の探索]
　U=F(U);
　V=F(V);
　V=F(V); //2回適用するのは意図的
　g=gcd(U-V,n);　//最大公約数を求める
　if(g==1) goto [因子の探索];
3[乱数種が悪かった場合]
　if(g==n) goto [乱数種の選択];
4[成功]
　return g;　//非自明な因子が見つかった

コード:

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

#define BASE 10
 
int main()
{

  int i,j,k;
  mpz_t a;
  mpz_t s;
  mpz_t d;
  mpz_t n;
  mpz_t z;
  mpz_t U;
  mpz_t V;
  mpz_t rnd;
  mpz_t temp;
  gmp_randstate_t rstate;
  
  mpz_init(a);
  mpz_init(s);
  mpz_init(d);
  mpz_init(n);
  mpz_init(z);
  mpz_init(U);
  mpz_init(V);
  mpz_init(rnd);
  mpz_init(temp);
  gmp_randinit_default(rstate);

  mpz_set_str(n,"11111111111111111111111111111111111111111111111111111",BASE);
  //Repunit R(53)
  gmp_printf("%Zd \n= ",n);

    mpz_sub_ui(n,n,3);         //n=n-3
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n-3)-1
    mpz_add_ui(a,rnd,1);       //a=rnd+1 == 1<a<n-3
    mpz_add_ui(n,n,4);         //n=n+4
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n+4)-1
    mpz_set(s,rnd);            //s = rnd == 0<s<n
    mpz_set(U,s);              //U=s
    mpz_set(V,s);              //V=s
    mpz_sub_ui(n,n,1);         //(n=n-1)

    while(1){//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_mul(temp,temp,a); //temp = temp * a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_mul(temp,temp,a); //temp = temp * a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_mul(temp,temp,a); //temp = temp * a 
      mpz_cdiv_r(V,temp,n); //   V = temp % n

      mpz_sub(temp,U,V);    //   temp = U - V
      
      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);

      i = mpz_cmp_ui(d,1);  //
      if(i<=0)
        continue;           //loop until d > 1
      break;
    }//2

  gmp_printf(" %Zd ",d);

  while(1){//3
    j=mpz_cmp(n,d);
    if(j<=0) break;            //if(n > d) break; 
    mpz_cdiv_q(n,n,d);         //n = n / d
    mpz_sub_ui(n,n,3);         //n = n - 3
    mpz_urandomm(rnd,rstate,n);
    mpz_add_ui(a,rnd,1);       //a = rnd + 1
    mpz_add_ui(n,n,4);         //n = n + 4
    mpz_urandomm(rnd,rstate,n);
    mpz_set(s,rnd);            //s = rnd
    mpz_set(U,s);              //U = s
    mpz_set(V,s);              //V = s
    mpz_sub_ui(n,n,1);         //n = n - 1
//ここまでは正しく動作した
    while(1){//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_mul(temp,temp,a); //temp = temp * a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_mul(temp,temp,a); //temp = temp * a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_mul(temp,temp,a); //temp = temp * a
      mpz_cdiv_r(V,temp,n); //   V = temp % n

      mpz_sub(temp,U,V);    //   temp = U - V

      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);
      
      k = mpz_cmp_ui(d,1);  //

      if(k<=0)
        continue;           //loop until d > 1
      break;
    }//2

    gmp_printf("* %Zd ",d);

  }//3

  mpz_clear(a);
  mpz_clear(s);
  mpz_clear(d);
  mpz_clear(n);
  mpz_clear(z);
  mpz_clear(U);
  mpz_clear(V);
  mpz_clear(rnd);
  mpz_clear(temp);

  return 0;

}

#2

理動さんが書きました:F(x)=(x * x + a) mod nと定める;

と書いてあるのに、ソースコードではF(x)=(x * x * a) mod nになっているようですが、どっちが正しいのでしょうか？

sleep · #3

・1 < d < n
・d == n
の場合の考慮は 73, 74行目の以下でやってあるようですが、

コード:

    j=mpz_cmp(n,d);
    if(j<=0) break;            //if(n > d) break;

103行目で d の値が１だったときは、どうするご予定でしょうか？

コード:

      k = mpz_cmp_ui(d,1);  //

今のところ、それが原因で無限ループに入ってます。

理動 · #4

みけCATさん
返信遅れてしまい申し訳ございません。
ご指摘ありがとうございます。
アルゴリズムの方が正しく、僕のソースコードが間違っています。

理動 · #5

sleepさん
返信遅れてしまい申し訳ございません。
ご指摘ありがとうございます。
d == 1の場合互いに素数ですが、
その場合に対処する部分がなかったので追加しておきます。

sleep · #6

理動さん、逆に質問させていただいてよろしいでしょうか？

ポラードのρ法って、私は使ったことがなくて
Wikiの解説読みながらちょっと自分でも書いてみたんです。
そうしたら・・・
普通の値（128など）で素因数分解してみると素数ではない値も
抽出されてしまうこともありまして、これって、最初の乱数の値によっては
ちゃんとした素数って求まらない場合もありますか？（汗） ←（これ質問です（笑））
というか、かなりブレるんですが・・・・

128 = 16 * 8
128 = 4 * 4 * 4 * 2
128 = 8 * 8 * 2
128 = 32 * 4

なんだろう、レピュニットの場合だけ抽出する値がかなり定まっている様に見えるんですが
こういうもんなんでしょうか・・・・
（むしろ、私が書いたプログラムが正しくないかもしれませんけど）

50桁だと私のPCだと物理コア４つフルで使ってもまるで終わる気がしなかったので
少し少な目の桁数でそれぞれ100回程度ずつテストはしてみました。
（理動さんどんなCPU使ってるんだろう・・・・）

コード:

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#define BASE 10

int main()
{
  int ret;
  mpz_t i;
  mpz_t s;
  mpz_t d;
  mpz_t n;
  mpz_t U;
  mpz_t V;
  mpz_t rnd;
  mpz_t upperbound;
  mpz_t quot;
  mpz_t rem;
  mpz_t temp;
  gmp_randstate_t rstate;

  mpz_init(i);
  mpz_init(s);
  mpz_init(d);
  mpz_init(n);
  mpz_init(U);
  mpz_init(V);
  mpz_init(rnd);
  mpz_init(upperbound);
  mpz_init(quot);
  mpz_init(rem);
  mpz_init(temp);
  gmp_randinit_default(rstate);

  //mpz_set_str(n,"11111111111111111111111111111111111111111111111111111",BASE);
  mpz_set_str(n,"111111111111111111111",BASE);
  //Repunit R(53)
  gmp_printf("%Zd \n=",n);

  mpz_urandomm(rnd,rstate,n);     //rnd
  mpz_set(s,rnd);                 //s = rnd
  mpz_set(U,s);                   //U = s
  mpz_mul(temp,s,s);              //
  mpz_add_ui(temp,temp,1);        //
  mpz_cdiv_r(V,temp,n);           //V = (s^2+1)%n
  mpz_sqrt(temp, n);              //
  mpz_add_ui(upperbound,temp,1);  //n = sqrt(n)+1

  mpz_set_ui(i,0);                // i = 0
  while(1)                        //for(int i=0; i<upperbound; i++)
  {
    mpz_add_ui(i,i,1);
    ret = mpz_cmp(i,upperbound);
    if(ret>=0) break;

    mpz_sub(temp,V,U);                      //
    mpz_gcd(temp,temp,n);                   //
    mpz_abs(d, temp);                       //d = abs(gcd((V-U),n));

    if((ret = mpz_cmp(d,n)) == 0) ;         //d == n
    else if((ret = mpz_cmp_ui(d,1)) != 0){  //d != 1
      while(1){
        mpz_cdiv_q(quot,n,d);               //quot = n/d
        mpz_cdiv_r(rem,n,d);                //rem = n%d
        ret = mpz_cmp_ui(rem,0);
        if(ret != 0) break;
        gmp_printf(" %Zd *",d);
        mpz_set(n,quot);                    //n = quot
      }
      mpz_sqrt(temp, n);                    //
      mpz_add_ui(upperbound,temp,1);        //upperbound = sqrt(n)+1
    }
    mpz_mul(temp,U,U);            //
    mpz_add_ui(temp,temp,1);      //
    mpz_cdiv_r(U,temp,n);         //U = (U^2+1)%n

    mpz_mul(temp,U,U);            //
    mpz_add_ui(temp,temp,1);      //
    mpz_mul(temp,temp,temp);      //
    mpz_add_ui(temp,temp,1);      //
    mpz_cdiv_r(V,temp,n);         //V = ((U^2+1)^2+1)%n
  }

  gmp_printf(" %Zd ",n);

  mpz_clear(i);
  mpz_clear(s);
  mpz_clear(d);
  mpz_clear(n);
  mpz_clear(U);
  mpz_clear(V);
  mpz_clear(rnd);
  mpz_clear(upperbound);
  mpz_clear(quot);
  mpz_clear(rem);
  mpz_clear(temp);

  return 0;
}

乱数使ってる時点でどうだったら正しいんだかわけがわからない・・・
何度か繰り返して統計取らないとダメなんですかね。

つまり、ρ法や楕円曲線法の場合

理動さんが書きました: 出力結果は
11111111111111111111111111111111111111111111111111111
= 107 * 1659431 * 1325815267337711173 * 47198858799491425660200071
となるのが正しいです。

これで安定して出ないような気がしてならないんですが・・・
（レピュニットに限り、１つずつ必ずバラバラの素数で構成されている、とかそういう法則があるんでしょうか？）

理動 · #7

sleepさん
夜遅くに返信ありがとうございます。
大切なお時間を割いていただき感謝します。

ポラードのρ法の注意点は
・プログラムで得られる各因数 d は必ずしも素数とは限らない。
・もとの数が出力されたからといって、それが素数とは限らない。
・求まる因数は必ずしも小さい順とは限らない
の3つがあるようです。

小さい数にρ法を適用した時にどのような動作をするかは、完全に僕の調査不足ですのでわかりません・・・

僕もR(53)以外の数字を試していたのですが、途中までですが安定した結果がでるのがR(53)でしたので
この数字を因数分解するコードで相談しようとトピックを立てました。
先ほど他のレピュニット数をいくつか試してみたところ、どれも一番小さい因数が最初に発見できていました。
何か法則がありそうですね。

さて、修正したコードを載せます。

コード:

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
 
#define BASE 10
 
int main()
{
 
  int i,j,k;
  mpz_t a;
  mpz_t s;
  mpz_t d;
  mpz_t n;
  mpz_t z;
  mpz_t U;
  mpz_t V;
  mpz_t rnd;
  mpz_t temp;
  gmp_randstate_t rstate;
  
  mpz_init(a);
  mpz_init(s);
  mpz_init(d);
  mpz_init(n);
  mpz_init(z);
  mpz_init(U);
  mpz_init(V);
  mpz_init(rnd);
  mpz_init(temp);
  gmp_randinit_default(rstate);
 
  mpz_set_str(n,"11111111111111111111111111111111111111111111111111111",BASE);
  //Repunit R(53)
  gmp_printf("%Zd \n= ",n);
 
    mpz_sub_ui(n,n,3);         //n=n-3
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n-3)-1
    mpz_add_ui(a,rnd,1);       //a=rnd+1 == 1<a<n-3
    mpz_add_ui(n,n,4);         //n=n+4
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n+4)-1
    mpz_set(s,rnd);            //s = rnd == 0<s<n
    mpz_set(U,s);              //U=s
    mpz_set(V,s);              //V=s
    mpz_sub_ui(n,n,1);         //(n=n-1)
 
    while(1){//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a 
      mpz_cdiv_r(V,temp,n); //   V = temp % n
 
      mpz_sub(temp,U,V);    //   temp = U - V
      
      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);
 
      i = mpz_cmp_ui(d,1);  //
      if(i<=0)
        continue;           //loop until d > 1
      break;
    }//2
 
  gmp_printf(" %Zd ",d);
 
  while(1){//3
    j=mpz_cmp(n,d);
    if(j<=0) break;            //if(n <= d) break; 
    mpz_cdiv_q(n,n,d);         //n = n / d
    mpz_sub_ui(n,n,3);         //n = n - 3
    mpz_urandomm(rnd,rstate,n);
    mpz_add_ui(a,rnd,1);       //a = rnd + 1
    mpz_add_ui(n,n,4);         //n = n + 4
    mpz_urandomm(rnd,rstate,n);
    mpz_set(s,rnd);            //s = rnd
    mpz_set(U,s);              //U = s
    mpz_set(V,s);              //V = s
    mpz_sub_ui(n,n,1);         //n = n - 1

    while(1){//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
 
      mpz_sub(temp,U,V);    //   temp = U - V
 
      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);
      
      k = mpz_cmp_ui(d,1);  //
 
      if(k<0)
        continue;           //loop until d > 1
      break;                 //if(d>=1) break
    }//2

    if(k!=0){
      gmp_printf(" * %Zd\n",d);
      continue;
    }
    //if(k==0) -> output n
    gmp_printf("* %Zd ",n);
    break;
 
  }//3
 
  mpz_clear(a);
  mpz_clear(s);
  mpz_clear(d);
  mpz_clear(n);
  mpz_clear(z);
  mpz_clear(U);
  mpz_clear(V);
  mpz_clear(rnd);
  mpz_clear(temp);
 
  return 0;
 
}

おそらく因数分解はできていると考えます。
ですが、

理動さんが書きました: 出力結果は
11111111111111111111111111111111111111111111111111111
= 107 * 1659431 * 1325815267337711173 * 47198858799491425660200071
となるのが正しいです。

とあるように、
一度しか因数を探さないのではなく、何度も探索したいのですがうまく動作しません。
何を修正したらいいでしょうか。

sleep · #8

理動さん、返答ありがとうございます。

理動さんが書きました: ポラードのρ法の注意点は
・プログラムで得られる各因数 d は必ずしも素数とは限らない。
・もとの数が出力されたからといって、それが素数とは限らない。
・求まる因数は必ずしも小さい順とは限らない
の3つがあるようです。

なるほど、勉強になります。

理動さんが書きました: 僕もR(53)以外の数字を試していたのですが、途中までですが安定した結果がでるのがR(53)でしたので
この数字を因数分解するコードで相談しようとトピックを立てました。

やはりそうなのですね。
安定している数字もあるのですね。

理動さんが書きました: 一度しか因数を探さないのではなく、何度も探索したいのですがうまく動作しません。
何を修正したらいいでしょうか。

今回の修正版のコードの場合、105～107行目のコードを

コード:

      if(k<0)
        continue;           //loop until d > 1
      break;                 //if(d>=1) break
    }//2

以下のように修正すると探索してくれるようになります。

コード:

      if(k>0) break;      //if(d>1) break
    }//2

ただ、

理動さんが書きました: おそらく因数分解はできていると考えます。

私の環境で試してみてましたが、最初の３つ目くらいまでは
安定して素数を求められているようなんですが、それ以降が
なんか途中から計算がおかしくなってますね。

コード:

  mpz_set_str(n,"111111111111111111111",BASE);

以下が理動さんのコードで実行した結果で

コード:

  = 3  * 37 * 43 * 239 * 4649 * 1933 * 10838689

以下が私のコードで実行した結果ですね。

コード:

  = 3 * 43 * 37 * 1933 * 12042986573479

12042986573479 が素数なんですけど、それを分解しようとしてるんですよね。
しかも、掛け直してみても元の数字にはならない、という。

n の求め方だけ私と大分違うので、多分 n じゃないかなぁ、とは思っているんですが
乱数を取り込みながらされてるようなので、ちょっとトレースが大変そうだったので諦めました。
探索は上記の修正でできるようになるので、確認してみていただけますか？

sleep · #9

12042986573479 が素数なんですけど、それを分解しようとしてるんですよね。
しかも、掛け直してみても元の数字にはならない、という。

いや、違いますね。
よくよく見ると、私が計算するときに１部コードが抜けてました。(37が抜けてました)
逆に、私のコードの方が因数分解できてないっぽいだけみたいですね。
申し訳ないです。

コード:

  //=  3  * 37 * 43 * 239 * 4649 * 1933 * 10838689
  mpz_set_ui(temp,3);
  mpz_mul_ui(temp,temp,37);       // 抜けてました
  mpz_mul_ui(temp,temp,43);
  mpz_mul_ui(temp,temp,239);
  mpz_mul_ui(temp,temp,4649);
  mpz_mul_ui(temp,temp,1933);
  mpz_mul_ui(temp,temp,10838689);
  gmp_printf("temp==%Zd ",temp);

  printf(\n);

  //= 3 * 43 * 37 * 1933 * 12042986573479
  mpz_set_ui(temp,3);
  mpz_mul_ui(temp,temp,43);
  mpz_mul_ui(temp,temp,37);
  mpz_mul_ui(temp,temp,1933);
  mpz_mul_ui(temp,temp,12042986573479);
  gmp_printf("temp==%Zd ",temp);

大丈夫そうですね。

#10

オフトピック

理動さんが書きました:

コード:

    while(1){//2
      // 省略

      i = mpz_cmp_ui(d,1);  //
      if(i<=0)
        continue;           //loop until d > 1
      break;
    }//2

この部分は

コード:

    do{//2
      // 省略

    }while((i=mpz_cmp_ui(d,1))<=0);//2

ではいけなかったのっでしょうか？

理動 · #11

sleepさん
返信ありがとうございます。

R(53)の場合、d=107の次にd=1となってしまい複数回の探索ができなくなっているようです。
2回目のdはd=1659431となってほしいのですが、うまくいきません。
ソースコードと出力結果を載せます。

コード:

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
 
#define BASE 10
 
int main()
{
 
  int i,j,k;
  mpz_t a;
  mpz_t s;
  mpz_t d;
  mpz_t n;
  mpz_t z;
  mpz_t U;
  mpz_t V;
  mpz_t rnd;
  mpz_t temp;
  gmp_randstate_t rstate;
  
  mpz_init(a);
  mpz_init(s);
  mpz_init(d);
  mpz_init(n);
  mpz_init(z);
  mpz_init(U);
  mpz_init(V);
  mpz_init(rnd);
  mpz_init(temp);
  gmp_randinit_default(rstate);
 
  mpz_set_str(n,"11111111111111111111111111111111111111111111111111111",BASE);
  //Repunit R(53)
  gmp_printf("%Zd \n= ",n);
 
    mpz_sub_ui(n,n,3);         //n=n-3
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n-3)-1
    mpz_add_ui(a,rnd,1);       //a=rnd+1 == 1<a<n-3
    mpz_add_ui(n,n,4);         //n=n+4
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n+4)-1
    mpz_set(s,rnd);            //s = rnd == 0<s<n
    mpz_set(U,s);              //U=s
    mpz_set(V,s);              //V=s
    mpz_sub_ui(n,n,1);         //(n=n-1)
 
    do{//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a 
      mpz_cdiv_r(V,temp,n); //   V = temp % n
 
      mpz_sub(temp,U,V);    //   temp = U - V
      
      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);
 
    }while((i=mpz_cmp_ui(d,1))<=0);//2 loop until d > 1
 
  gmp_printf(" %Zd ",d);
 
  while((j=mpz_cmp(n,d))>0){//3 if(n<=d) break
    mpz_cdiv_q(n,n,d);         //n = n / d
    mpz_sub_ui(n,n,3);         //n = n - 3
    mpz_urandomm(rnd,rstate,n);
    mpz_add_ui(a,rnd,1);       //a = rnd + 1
    mpz_add_ui(n,n,4);         //n = n + 4
    mpz_urandomm(rnd,rstate,n);
    mpz_set(s,rnd);            //s = rnd
    mpz_set(U,s);              //U = s
    mpz_set(V,s);              //V = s
    mpz_sub_ui(n,n,1);         //n = n - 1

    do{//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
 
      mpz_sub(temp,U,V);    //   temp = U - V
 
      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);
      
    }while((k=mpz_cmp_ui(d,1))<=0);//2 loop until d > 1

    if(k=0){
      gmp_printf(" * %Zd\n",n);
      break;
    }
    gmp_printf("* %Zd ",d);
    break;
  }//3
 
  mpz_clear(a);
  mpz_clear(s);
  mpz_clear(d);
  mpz_clear(n);
  mpz_clear(z);
  mpz_clear(U);
  mpz_clear(V);
  mpz_clear(rnd);
  mpz_clear(temp);
 
  return 0;
 
}

$ ./pollard
111111111111111111111 (R21)
= 3 * 37

$ ./pollard
11111111111111111111111111111111111111111111111111111 (R53)
= 107 * 1659431

69行目と78行目のnの値は一致しているのを確認しましたので、
tempの内容が正しい値になっていないのかと考えています。

理動 · #12

みけCATさん
ご指摘ありがとうございます。

似たような箇所すべて修正しました。
多少見やすいコードになったかと思います。
ありがとうございます。

sleep · #13

理動さん、私はそのような形の修正方法を提示してません。
私は以下のように修正してください、と書きましたよね。

sleep さんが書きました: 今回の修正版のコードの場合、105～107行目のコードを
コード:
      if(k<0)
        continue;           //loop until d > 1
      break;                 //if(d>=1) break
    }//2
以下のように修正すると探索してくれるようになります。
コード:
      if(k>0) break;      //if(d>1) break
    }//2

理動 · #14

sleepさん
返信ありがとうございます。

大変申し訳ございません、sleepさんの指示通りの修正で正しく動作するものになっていました。
僕が処理に時間がかかっているものを無限ループから抜け出せなくなっていると勘違いしていたため、
コードをさらに改変してしまいました。

これでポラードのρ法による因数分解のプログラムが完成しました。
ご協力いただきありがとうございました。

またわからないことがありましたらトピックを立てて質問いたしますので、
そのときはまたよろしくお願いします。

最後に、少々見にくいですがソースコードを載せておきます。
お役立ていただければ幸いです。

コード:

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
 
#define BASE 10
 
int main()
{
 
  int i,j,k;
  mpz_t a;
  mpz_t s;
  mpz_t d;
  mpz_t n;
  mpz_t z;
  mpz_t U;
  mpz_t V;
  mpz_t rnd;
  mpz_t temp;
  gmp_randstate_t rstate;
  
  mpz_init(a);
  mpz_init(s);
  mpz_init(d);
  mpz_init(n);
  mpz_init(z);
  mpz_init(U);
  mpz_init(V);
  mpz_init(rnd);
  mpz_init(temp);
  gmp_randinit_default(rstate);
 
  mpz_set_str(n,"11111111111111111111111111111111111111111111111111111",BASE);
  //Repunit R(53)
  gmp_printf("%Zd \n= ",n);
 
    mpz_sub_ui(n,n,3);         //n=n-3
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n-3)-1
    mpz_add_ui(a,rnd,1);       //a=rnd+1 == 1<a<n-3
    mpz_add_ui(n,n,4);         //n=n+4
    mpz_urandomm(rnd,rstate,n);//0<rnd<(n+4)-1
    mpz_set(s,rnd);            //s = rnd == 0<s<n
    mpz_set(U,s);              //U=s
    mpz_set(V,s);              //V=s
    mpz_sub_ui(n,n,1);         //(n=n-1)
 
    while(1){//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a 
      mpz_cdiv_r(V,temp,n); //   V = temp % n
 
      mpz_sub(temp,U,V);    //   temp = U - V
      
      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);
 
      i = mpz_cmp_ui(d,1);  //
      if(i<=0)
        continue;           //loop until d > 1
      break;
    }//2
 
  gmp_printf(" %Zd ",d);
 
  while(1){//3
    j=mpz_cmp(n,d);
    if(j<=0) break;            //if(n <= d) break; 
    mpz_cdiv_q(n,n,d);         //n = n / d
    mpz_sub_ui(n,n,3);         //n = n - 3
    mpz_urandomm(rnd,rstate,n);
    mpz_add_ui(a,rnd,1);       //a = rnd + 1
    mpz_add_ui(n,n,4);         //n = n + 4
    mpz_urandomm(rnd,rstate,n);
    mpz_set(s,rnd);            //s = rnd
    mpz_set(U,s);              //U = s
    mpz_set(V,s);              //V = s
    mpz_sub_ui(n,n,1);         //n = n - 1
 
    while(1){//2
      mpz_mul(temp,U,U);    //temp = U * U
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(U,temp,n); //   U = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
      
      mpz_mul(temp,V,V);    //temp = V * V
      mpz_add(temp,temp,a); //temp = temp + a
      mpz_cdiv_r(V,temp,n); //   V = temp % n
 
      mpz_sub(temp,U,V);    //   temp = U - V
 
      mpz_gcd(d,temp,n);    //   d = gcd(temp,n);
      
      k = mpz_cmp_ui(d,1);  //
 
      if(k>0) break;        //if(d>1) break
    }//2
 
    if(k!=0){
      gmp_printf(" * %Zd\n",d);
      continue;
    }
    //if(k==0) -> output n
    gmp_printf("* %Zd ",n);
    break;
 
  }//3
 
  mpz_clear(a);
  mpz_clear(s);
  mpz_clear(d);
  mpz_clear(n);
  mpz_clear(z);
  mpz_clear(U);
  mpz_clear(V);
  mpz_clear(rnd);
  mpz_clear(temp);
 
  return 0;
 
}

GNU MP(GMP)を用いたプログラミングについて

GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

Re: GNU MP(GMP)を用いたプログラミングについて

		7月 2025
日	月	火	水	木	金	土
		1	2	3	4	5
6	7	8	9	10	11	12
13	14	15	16	17	18	19
20	21	22	23	24	25	26
27	28	29	30	31