素因数分解

 投稿者:しばっち  投稿日:2016年 6月18日(土)20時44分50秒
  試し割り法 + ρ法で素因数分解をします。
ISPRIME(X)はミラー・ラビン法などの素数判定ルーチンです。ルーチンは省略しています。

ポラード・ロー法(ρ法)をC++で書いています。
また、boostライブラリーを使い、任意精度整数型(cpp_int型)を使用しています。

ビルドしたdllファイルは大き過ぎるので、今回も掲示板への投稿は控えさせて頂きます。
C++開発環境を持たない方のため、VC++2015(x86)にてビルドしたファイルをアップローダーへ
upしておきました。(factor.zip) 下記URLは現時点より1ヶ月間有効です。

ダウンロードパス:shibacchi

試しに2^N-1 として Nを100~120(2^100-1 ~ 2^120-1)までを素因数分解してみました。

OPTION ARITHMETIC RATIONAL
LET L=1000000
DIM S(L)
FOR I=2 TO L !'エラトステネスの篩
   IF S(I)=0 THEN
      FOR J=I*I TO L STEP I
         LET S(J)=1
      NEXT J
   END IF
NEXT I
FOR K=100 TO 120
   LET N=2^K-1
   PRINT "2 ^";K;"- 1"
   PRINT N;"-->";
   IF ISPRIME(N)=0 THEN !'ミラー・ラビン法など。合成数なら試し割り
      LET NN=N
      FOR I=2 TO L
         IF S(I)=0 THEN
            DO WHILE MOD(N,I)=0
               LET N=N/I
               PRINT I;
            LOOP
         END IF
      NEXT I
      IF NN<>N AND N>1 THEN PRINT N ELSE PRINT
      CALL FACTOR(N)
   ELSE
      PRINT " 素数"
   END IF
   PRINT
NEXT K
END

EXTERNAL SUB FACTOR(N)
OPTION ARITHMETIC RATIONAL
IF N>1 AND ISPRIME(N)=0 THEN !'ミラー・ラビン法などで確認。合成数なら
   PRINT N;"-->";
   LET D=POLLARDRHO(N)
   PRINT D;N/D
   CALL FACTOR(D) !'再帰呼び出し
   CALL FACTOR(N/D)
END IF
END SUB

EXTERNAL  FUNCTION POLLARDRHO(X) !'ポラード・ロー素因数分解法(ρ法)
OPTION ARITHMETIC RATIONAL
LET N$=STR$(X)
CALL POLLARD(N$)
FOR I=1 TO LEN(N$) !' POLLARDRHO=VAL(N$) とできない
   LET P$=N$(I:I)
   IF ORD(P$)>=48 AND ORD(P$)<=57 THEN LET B$=B$ & P$
NEXT I
LET POLLARDRHO=VAL(B$)
END FUNCTION

EXTERNAL SUB POLLARD(X$) !'素数を入れないこと。無限ループのおそれ。
OPTION ARITHMETIC RATIONAL
ASSIGN "pollard.dll","pollard"
END SUB
----------------------------------------------------------------------------
                                  pollard.cpp
#include <boost/multiprecision/cpp_int.hpp>

using namespace boost::multiprecision;
using namespace std;

cpp_int atocpp(char str[])
{
cpp_int result = 0;
int i=0;
while (str[i]>='0' && str[i]<='9'){
    result=result*10+(str[i++] - '0');
                  }
return result;
}

cpp_int mulmod(cpp_int a,cpp_int b,cpp_int m)
{
cpp_int result=0;
a %= m;
while (b>0) {
    if  (b % 2 == 1) result = (result+a) % m;
    a = (a*2) % m;
    b = b / 2;
                 }
return result;
}

cpp_int pollardrho(cpp_int num,unsigned long c){
  cpp_int x=1,y,p=1,n=2,i=1;
  while (p == 1){
    ++i;
    if (i == n){
      y = x;
      n *= 2;
    }
    x = (mulmod(x,x,num) + c) % num;
    if (y == x) return pollardrho(num,c+1);
    p = gcd(num,abs(x-y));
  }
  return p;
}

extern "C" __declspec(dllexport) void pollard(char x[])
{
cpp_int num,n;
int i=0,tmp,len;
num=atocpp(x);
len=strlen(x);
for(int j=0;j<len;j++) x[j]='\0';
n=pollardrho(num,1u);
while (n != 0) {
    x[i] = (char)(n % 10) + '0';
    n /= 10;
    i++;
            }
x[i] = '\0';
i--;
for(int j = 0;j <= i/2;j++) {
    tmp = x[j];
    x[j] = x[i-j];
    x[i-j] = tmp;
                       }

}

実行結果
2 ^ 100 - 1
1267650600228229401496703205375 --> 3  5  5  5  11  31  41  101  251  601  1801  4051  8101  268501

2 ^ 101 - 1
2535301200456458802993406410751 -->
2535301200456458802993406410751 --> 7432339208719  341117531003194129

2 ^ 102 - 1
5070602400912917605986812821503 --> 3  3  7  103  307  2143  2857  6529  11119  43691  131071

2 ^ 103 - 1
10141204801825835211973625643007 -->
10141204801825835211973625643007 --> 2550183799  3976656429941438590393

2 ^ 104 - 1
20282409603651670423947251286015 --> 3  5  17  53  157  1613  2731  8191  858001  308761441

2 ^ 105 - 1
40564819207303340847894502572031 --> 7  7  31  71  127  151  337  29191  106681  122921  152041

2 ^ 106 - 1
81129638414606681695789005144063 --> 3  107  6361  69431  572263032673174337633
572263032673174337633 --> 20394401  28059810762433

2 ^ 107 - 1
162259276829213363391578010288127 --> 素数

2 ^ 108 - 1
324518553658426726783156020576255 --> 3  3  3  3  5  7  13  19  37  73  109  87211  246241  262657  279073

2 ^ 109 - 1
649037107316853453566312041152511 -->
649037107316853453566312041152511 --> 745988807  870035986098720987332873

2 ^ 110 - 1
1298074214633706907132624082305023 --> 3  11  11  23  31  89  683  881  2971  3191  201961  48912491

2 ^ 111 - 1
2596148429267413814265248164610047 --> 7  223  321679  5170159074941697420729913
5170159074941697420729913 --> 26295457  196617958567584409
196617958567584409 --> 616318177  319020217

2 ^ 112 - 1
5192296858534827628530496329220095 --> 3  5  17  29  43  113  127  257  5153  859166727965929937
859166727965929937 --> 15790321  54410972897

2 ^ 113 - 1
10384593717069655257060992658440191 --> 3391  23279  65993  1993423291715412685783
1993423291715412685783 --> 1868569  1066818132868207

2 ^ 114 - 1
20769187434139310514121985316880383 --> 3  3  7  571  32377  174763  524287  194620086937183
194620086937183 --> 1212847  160465489

2 ^ 115 - 1
41538374868278621028243970633760767 --> 31  47  14951  178481  10683848415441845139401
10683848415441845139401 --> 4036961  2646507710984041

2 ^ 116 - 1
83076749736557242056487941267521535 --> 3  5  59  233  1103  2089  174850288360342272373981
174850288360342272373981 --> 3033169  57646075230342349
57646075230342349 --> 107367629  536903681

2 ^ 117 - 1
166153499473114484112975882535043071 --> 7  73  79  937  6553  8191  86113  121369  7830118297

2 ^ 118 - 1
332306998946228968225951765070086143 --> 3  2833  37171  179951  5845385390147915655817
5845385390147915655817 --> 1824726041  3203431780337

2 ^ 119 - 1
664613997892457936451903530140172287 --> 127  239  20231  131071  8257410955834335790279
8257410955834335790279 --> 131105292137  62983048367

2 ^ 120 - 1
1329227995784915872903807060280344575 --> 3  3  5  5  7  11  13  17  31  41  61  151  241  331  1321  61681  4562284561

 

戻る