ミラーラビン法

 投稿者:しばっち  投稿日:2016年 6月13日(月)22時08分2秒
  以前投稿したミラーラビン法をC++で書いてみました。
#3894

多倍長計算にboostライブラリーを使用しました。

無限精度整数型(cpp_int型)を使用しています。
第2引数(NUM)は判定回数(精度)です。素数リストを1000個用意しましたので
2~1000までの数値を指定してください。

また、ビルドしたdllファイルは想定していたより大き過ぎたので
今回は掲示板への投稿は控えさせて頂きます。
しかしながら、C++開発環境を持たない方のため、VC++2015(x86)にてビルド(コンパイル)
したファイルをアップローダーにupしておきました。(miller rabin.zip)

ダウンロードパス:shibacchi
なお、上記URLには有効期限があり、現時点より1ヶ月間となります。

OPTION ARITHMETIC RATIONAL
LET S=10^500+1
LET E=S+10000
FOR I=S TO E STEP 2
   LET A$=STR$(I)
   IF ISPRIME(A$,100)<>0 THEN
      PRINT A$
   END IF
NEXT  I
END

EXTERNAL FUNCTION ISPRIME(N$,NUM)
OPTION ARITHMETIC RATIONAL
ASSIGN "miller.dll","isprime"
END FUNCTION

EXTERNAL FUNCTION ISPRIME2(N$,NUM)
OPTION ARITHMETIC RATIONAL
ASSIGN "isprime.dll","isprime"
END FUNCTION
---------------------------------------------------------------------------------
                                              miller.cpp

#include <boost/multiprecision/cpp_int.hpp>

using namespace boost::multiprecision;
using namespace std;

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

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

extern "C" __declspec(dllexport) int isprime(char *x,int num)
{
cpp_int d,r,n=0;
int i,j,isp,s=0;
int a[]={
  2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103
,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211
,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331
,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449
,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587
,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709
,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853
,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991
,997,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,1093
,1097,1103,1109,1117,1123,1129,1151,1153,1163,1171,1181,1187,1193,1201,1213,1217
,1223,1229,1231,1237,1249,1259,1277,1279,1283,1289,1291,1297,1301,1303,1307,1319
,1321,1327,1361,1367,1373,1381,1399,1409,1423,1427,1429,1433,1439,1447,1451,1453
,1459,1471,1481,1483,1487,1489,1493,1499,1511,1523,1531,1543,1549,1553,1559,1567
,1571,1579,1583,1597,1601,1607,1609,1613,1619,1621,1627,1637,1657,1663,1667,1669
,1693,1697,1699,1709,1721,1723,1733,1741,1747,1753,1759,1777,1783,1787,1789,1801
,1811,1823,1831,1847,1861,1867,1871,1873,1877,1879,1889,1901,1907,1913,1931,1933
,1949,1951,1973,1979,1987,1993,1997,1999,2003,2011,2017,2027,2029,2039,2053,2063
,2069,2081,2083,2087,2089,2099,2111,2113,2129,2131,2137,2141,2143,2153,2161,2179
,2203,2207,2213,2221,2237,2239,2243,2251,2267,2269,2273,2281,2287,2293,2297,2309
,2311,2333,2339,2341,2347,2351,2357,2371,2377,2381,2383,2389,2393,2399,2411,2417
,2423,2437,2441,2447,2459,2467,2473,2477,2503,2521,2531,2539,2543,2549,2551,2557
,2579,2591,2593,2609,2617,2621,2633,2647,2657,2659,2663,2671,2677,2683,2687,2689
,2693,2699,2707,2711,2713,2719,2729,2731,2741,2749,2753,2767,2777,2789,2791,2797
,2801,2803,2819,2833,2837,2843,2851,2857,2861,2879,2887,2897,2903,2909,2917,2927
,2939,2953,2957,2963,2969,2971,2999,3001,3011,3019,3023,3037,3041,3049,3061,3067
,3079,3083,3089,3109,3119,3121,3137,3163,3167,3169,3181,3187,3191,3203,3209,3217
,3221,3229,3251,3253,3257,3259,3271,3299,3301,3307,3313,3319,3323,3329,3331,3343
,3347,3359,3361,3371,3373,3389,3391,3407,3413,3433,3449,3457,3461,3463,3467,3469
,3491,3499,3511,3517,3527,3529,3533,3539,3541,3547,3557,3559,3571,3581,3583,3593
,3607,3613,3617,3623,3631,3637,3643,3659,3671,3673,3677,3691,3697,3701,3709,3719
,3727,3733,3739,3761,3767,3769,3779,3793,3797,3803,3821,3823,3833,3847,3851,3853
,3863,3877,3881,3889,3907,3911,3917,3919,3923,3929,3931,3943,3947,3967,3989,4001
,4003,4007,4013,4019,4021,4027,4049,4051,4057,4073,4079,4091,4093,4099,4111,4127
,4129,4133,4139,4153,4157,4159,4177,4201,4211,4217,4219,4229,4231,4241,4243,4253
,4259,4261,4271,4273,4283,4289,4297,4327,4337,4339,4349,4357,4363,4373,4391,4397
,4409,4421,4423,4441,4447,4451,4457,4463,4481,4483,4493,4507,4513,4517,4519,4523
,4547,4549,4561,4567,4583,4591,4597,4603,4621,4637,4639,4643,4649,4651,4657,4663
,4673,4679,4691,4703,4721,4723,4729,4733,4751,4759,4783,4787,4789,4793,4799,4801
,4813,4817,4831,4861,4871,4877,4889,4903,4909,4919,4931,4933,4937,4943,4951,4957
,4967,4969,4973,4987,4993,4999,5003,5009,5011,5021,5023,5039,5051,5059,5077,5081
,5087,5099,5101,5107,5113,5119,5147,5153,5167,5171,5179,5189,5197,5209,5227,5231
,5233,5237,5261,5273,5279,5281,5297,5303,5309,5323,5333,5347,5351,5381,5387,5393
,5399,5407,5413,5417,5419,5431,5437,5441,5443,5449,5471,5477,5479,5483,5501,5503
,5507,5519,5521,5527,5531,5557,5563,5569,5573,5581,5591,5623,5639,5641,5647,5651
,5653,5657,5659,5669,5683,5689,5693,5701,5711,5717,5737,5741,5743,5749,5779,5783
,5791,5801,5807,5813,5821,5827,5839,5843,5849,5851,5857,5861,5867,5869,5879,5881
,5897,5903,5923,5927,5939,5953,5981,5987,6007,6011,6029,6037,6043,6047,6053,6067
,6073,6079,6089,6091,6101,6113,6121,6131,6133,6143,6151,6163,6173,6197,6199,6203
,6211,6217,6221,6229,6247,6257,6263,6269,6271,6277,6287,6299,6301,6311,6317,6323
,6329,6337,6343,6353,6359,6361,6367,6373,6379,6389,6397,6421,6427,6449,6451,6469
,6473,6481,6491,6521,6529,6547,6551,6553,6563,6569,6571,6577,6581,6599,6607,6619
,6637,6653,6659,6661,6673,6679,6689,6691,6701,6703,6709,6719,6733,6737,6761,6763
,6779,6781,6791,6793,6803,6823,6827,6829,6833,6841,6857,6863,6869,6871,6883,6899
,6907,6911,6917,6947,6949,6959,6961,6967,6971,6977,6983,6991,6997,7001,7013,7019
,7027,7039,7043,7057,7069,7079,7103,7109,7121,7127,7129,7151,7159,7177,7187,7193
,7207,7211,7213,7219,7229,7237,7243,7247,7253,7283,7297,7307,7309,7321,7331,7333
,7349,7351,7369,7393,7411,7417,7433,7451,7457,7459,7477,7481,7487,7489,7499,7507
,7517,7523,7529,7537,7541,7547,7549,7559,7561,7573,7577,7583,7589,7591,7603,7607
,7621,7639,7643,7649,7669,7673,7681,7687,7691,7699,7703,7717,7723,7727,7741,7753
,7757,7759,7789,7793,7817,7823,7829,7841,7853,7867,7873,7877,7879,7883,7901,7907
,7919 };
if (num>1000) num=1000;
if (num<2) num=2;
n=atocpp(x);
if (n<=1) return 0;
for (i=0;i<=num-1;i++)
if (a[i]==n) return 1;
else if (n % a[i]==0) return 0;
d=(n-1)/2;
while (d % 2==0) {
    d=d/2;
    s++;
                        }
for(i=0;i<=num-1;i++){
    isp=0;
    r=powmod(a[i],d,n);
    if (r==1 || r==n-1) isp=1;
    r=powmod(r,2,n);
    for(j=0;j<=s-1;j++){
        if (r==n-1) isp=1;
        r=powmod(r,2,n);
                       }
    if (isp==0) return 0;
                    }
return 1;
  }

また、boostライブラリーには下記ソースのようにmiller-rabin法が定義されています。
こちらは素数リストではなく、乱数を使用しているようです。
第2引数(NUM)は判定回数を指定してください。
------------------------------------------------------------------------------------
                                            isprime.cpp

#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>

using namespace boost::multiprecision;
using namespace std;

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

extern "C" __declspec(dllexport) int isprime(char *x,int num)
{
cpp_int n;
n=atocpp(x);
if (n==2) return 1;
    bool is_prime = miller_rabin_test(n, num);
    return (is_prime ? 1:0);
}
 

Re: ミラーラビン法

 投稿者:たろさ  投稿日:2016年 6月16日(木)09時45分48秒
  > No.4085[元記事へ]

しばっちさんへのお返事です。

> 以前投稿したミラーラビン法をC++で書いてみました。

素晴らしいプログラムを、ありがとうございます。

動作報告です。

OPTION ARITHMETIC RATIONAL     !有理数モード
LET t0=TIME
LET Q=7 !1020379 (80000th prime)【変数】
LET k9=80000
DIM Pz(k9)

OPEN #2:NAME "E:\prime_10_8x1v.txt",ACCESS INPUT !素数リスト
FOR i=1 TO k9
   INPUT #2: Pz(i)
NEXT i
CLOSE #2

LET m1=1
LET m2=1
LET m3=1
LET m4=1
LET m5=1
LET m6=1
LET m7=1
LET m8=1
LET m9=1
LET m10=1
LET m11=1
LET m12=1
LET m13=1
LET m14=1
LET m15=1
LET m16=1
LET m17=1
LET m18=1
LET m19=1
LET m20=1
LET m21=1
LET m22=1
LET m23=1
LET m24=1
LET m25=1
LET m26=1
LET m27=1
LET m28=1
LET m29=1
LET m30=1
LET m31=1
LET m32=1
LET m33=1
LET m34=1
LET m35=1
LET m36=1
LET m37=1
LET m38=1
LET m39=1
LET m40=1
LET m41=1
LET m42=1
LET m43=1
LET m44=1
LET m45=1
LET m46=1
LET m47=1
LET m48=1
LET m49=1
LET m50=1
LET m51=1
LET m52=1
LET m53=1
LET m54=1
LET m55=1
LET m56=1
LET m57=1
LET m58=1
LET m59=1
LET m60=1
LET m61=1
LET m62=1
LET m63=1
LET m64=1
LET m65=1
LET m66=1
LET m67=1
LET m68=1
LET m69=1
LET m70=1
LET m71=1
LET m72=1
LET m73=1
LET m74=1
LET m75=1
LET m76=1
LET m77=1
LET m78=1
LET m79=1
LET m80=1

FOR iz=1 TO k9
   LET z=MOD(iz,40)
   SELECT CASE z
   CASE 1
      LET m1=m1*Pz(iz)
   CASE 2
      LET m2=m2*Pz(iz)
   CASE 3
      LET m3=m3*Pz(iz)
   CASE 4
      LET m4=m4*Pz(iz)
   CASE 5
      LET m5=m5*Pz(iz)
   CASE 6
      LET m6=m6*Pz(iz)
   CASE 7
      LET m7=m7*Pz(iz)
   CASE 8
      LET m8=m8*Pz(iz)
   CASE 9
      LET m9=m9*Pz(iz)
   CASE 10
      LET m10=m10*Pz(iz)
   CASE 11
      LET m11=m11*Pz(iz)
   CASE 12
      LET m12=m12*Pz(iz)
   CASE 13
      LET m13=m13*Pz(iz)
   CASE 14
      LET m14=m14*Pz(iz)
   CASE 15
      LET m15=m15*Pz(iz)
   CASE 16
      LET m16=m16*Pz(iz)
   CASE 17
      LET m17=m17*Pz(iz)
   CASE 18
      LET m18=m18*Pz(iz)
   CASE 19
      LET m19=m19*Pz(iz)
   CASE 20
      LET m20=m20*Pz(iz)
   CASE 21
      LET m21=m21*Pz(iz)
   CASE 22
      LET m22=m22*Pz(iz)
   CASE 23
      LET m23=m23*Pz(iz)
   CASE 24
      LET m24=m24*Pz(iz)
   CASE 25
      LET m25=m25*Pz(iz)
   CASE 26
      LET m26=m26*Pz(iz)
   CASE 27
      LET m27=m27*Pz(iz)
   CASE 28
      LET m28=m28*Pz(iz)
   CASE 29
      LET m29=m29*Pz(iz)
   CASE 30
      LET m30=m30*Pz(iz)
   CASE 31
      LET m31=m31*Pz(iz)
   CASE 32
      LET m32=m32*Pz(iz)
   CASE 33
      LET m33=m33*Pz(iz)
   CASE 34
      LET m34=m34*Pz(iz)
   CASE 35
      LET m35=m35*Pz(iz)
   CASE 36
      LET m36=m36*Pz(iz)
   CASE 37
      LET m37=m37*Pz(iz)
   CASE 38
      LET m38=m38*Pz(iz)
   CASE 39
      LET m39=m39*Pz(iz)
   CASE 40
      LET m40=m40*Pz(iz)
   CASE 41
      LET m41=m41*Pz(iz)
   CASE 42
      LET m42=m42*Pz(iz)
   CASE 43
      LET m43=m43*Pz(iz)
   CASE 44
      LET m44=m44*Pz(iz)
   CASE 45
      LET m45=m45*Pz(iz)
   CASE 46
      LET m46=m46*Pz(iz)
   CASE 47
      LET m47=m47*Pz(iz)
   CASE 48
      LET m48=m48*Pz(iz)
   CASE 49
      LET m49=m49*Pz(iz)
   CASE 50
      LET m50=m50*Pz(iz)
   CASE 51
      LET m51=m51*Pz(iz)
   CASE 52
      LET m52=m52*Pz(iz)
   CASE 53
      LET m53=m53*Pz(iz)
   CASE 54
      LET m54=m54*Pz(iz)
   CASE 55
      LET m55=m55*Pz(iz)
   CASE 56
      LET m56=m56*Pz(iz)
   CASE 57
      LET m57=m57*Pz(iz)
   CASE 58
      LET m58=m58*Pz(iz)
   CASE 59
      LET m59=m59*Pz(iz)
   CASE 60
      LET m60=m60*Pz(iz)
   CASE 61
      LET m61=m61*Pz(iz)
   CASE 62
      LET m62=m62*Pz(iz)
   CASE 63
      LET m63=m63*Pz(iz)
   CASE 64
      LET m64=m64*Pz(iz)
   CASE 65
      LET m65=m65*Pz(iz)
   CASE 66
      LET m66=m66*Pz(iz)
   CASE 67
      LET m67=m67*Pz(iz)
   CASE 68
      LET m68=m68*Pz(iz)
   CASE 69
      LET m69=m69*Pz(iz)
   CASE 70
      LET m70=m70*Pz(iz)
   CASE 71
      LET m71=m71*Pz(iz)
   CASE 72
      LET m72=m72*Pz(iz)
   CASE 73
      LET m73=m73*Pz(iz)
   CASE 74
      LET m74=m74*Pz(iz)
   CASE 75
      LET m75=m75*Pz(iz)
   CASE 76
      LET m76=m76*Pz(iz)
   CASE 77
      LET m77=m77*Pz(iz)
   CASE 78
      LET m78=m78*Pz(iz)
   CASE 79
      LET m79=m79*Pz(iz)
   CASE 0
      LET m80=m80*Pz(iz)
   END SELECT
NEXT iz
LET C=3
FOR i=Q TO Q+100 STEP 1
   PRINT "【";STR$(i)&"】";
   LET N=(2*10^i+1)/3
   IF n=NUMER(n/m1) AND n=NUMER(n/m2) AND n=NUMER(n/m3) AND n=NUMER(n/m4) AND n=NUMER(n/m5) AND n=NUMER(n/m6) AND n=NUMER(n/m7) AND n=NUMER(n/m8) AND n=NUMER(n/m9) AND n=NUMER(n/m10) AND n=NUMER(n/m11) AND n=NUMER(n/m12) AND n=NUMER(n/m13) AND n=NUMER(n/m14) AND n=NUMER(n/m15) AND n=NUMER(n/m16) AND n=NUMER(n/m17) AND n=NUMER(n/m18) AND n=NUMER(n/m19) AND n=NUMER(n/m20) AND n=NUMER(n/m21) AND n=NUMER(n/m22) AND n=NUMER(n/m23) AND n=NUMER(n/m24) AND n=NUMER(n/m25) AND n=NUMER(n/m26) AND n=NUMER(n/m27) AND n=NUMER(n/m28) AND n=NUMER(n/m29) AND n=NUMER(n/m30) AND n=NUMER(n/m31) AND n=NUMER(n/m32) AND n=NUMER(n/m33) AND n=NUMER(n/m34) AND n=NUMER(n/m35) AND n=NUMER(n/m36) AND n=NUMER(n/m37) AND n=NUMER(n/m38) AND n=NUMER(n/m39) AND n=NUMER(n/m40) AND n=NUMER(n/m41) AND n=NUMER(n/m42) AND n=NUMER(n/m43) AND n=NUMER(n/m44) AND n=NUMER(n/m45) AND n=NUMER(n/m46) AND n=NUMER(n/m47) AND n=NUMER(n/m48) AND n=NUMER(n/m49) AND n=NUMER(n/m50) AND n=NUMER(n/m51) AND n=NUMER(n/m52) AND n=NUMER(n/m53) AND n=NUMER(n/m54) AND n=NUMER(n/m55) AND n=NUMER(n/m56) AND n=NUMER(n/m57) AND n=NUMER(n/m58) AND n=NUMER(n/m59) AND n=NUMER(n/m60) AND n=NUMER(n/m61) AND n=NUMER(n/m62) AND n=NUMER(n/m63) AND n=NUMER(n/m64) AND n=NUMER(n/m65) AND n=NUMER(n/m66) AND n=NUMER(n/m67) AND n=NUMER(n/m68) AND n=NUMER(n/m69) AND n=NUMER(n/m70) AND n=NUMER(n/m71) AND n=NUMER(n/m72) AND n=NUMER(n/m73) AND n=NUMER(n/m74) AND n=NUMER(n/m75) AND n=NUMER(n/m76) AND n=NUMER(n/m77) AND n=NUMER(n/m78) AND n=NUMER(n/m79) AND n=NUMER(n/m80) THEN
      LET A$=STR$(n)
      IF ISPRIME(A$,100)<>0 THEN
         PRINT A$;
         LET C=C+1
         PRINT ":";c;"個"
         GOTO 100
      ELSE
         LET Cx=Cx+1
         PRINT "Miller-Rabin法";":";Cx;"回"
         GOTO 100
      END IF
   ELSE
      LET Cy=Cy+1
      PRINT "合成数";":";Cy;"回"
      GOTO 100
   END IF
   PRINT
100 NEXT i

    LET TM=TIME-t0
    PRINT USING"####." & REPEAT$("#",2):TM;
    PRINT "秒"
END

EXTERNAL FUNCTION ISPRIME(N$,NUM)
    OPTION ARITHMETIC RATIONAL
    ASSIGN "miller.dll","isprime"
END FUNCTION

-------------------------------------------------

OPTION ARITHMETIC RATIONAL     !有理数モード
LET t0=TIME
FOR i=1 TO 17
   PRINT STR$(i)&"=";
   LET n=(2*10^i+1)/3
   LET A$=STR$(n)
   IF ISPRIME(A$,100)<>0 THEN
      PRINT A$
   ELSE
      DO WHILE MOD(n,2)=0
         PRINT STR$(2)&"*";
         LET n=n/2
      LOOP
      LET f=3
      LET s=INT(SQR(n))
      DO UNTIL f>s
         DO WHILE MOD(n,f)=0
            PRINT STR$(f)&"*";
            LET n=n/f
            LET A$=STR$(n)
            IF ISPRIME(A$,100)<>0 THEN
               PRINT A$
               GOTO 100
            END IF
            LET s=INT(SQR(n))
         LOOP
         LET f=f+2
      LOOP
      IF n>1 THEN PRINT STR$(n)
   END IF
100 NEXT i
    LET TM=TIME-t0
    PRINT USING"####." & REPEAT$("#",2):TM;
    PRINT "秒"

END

EXTERNAL FUNCTION ISPRIME(N$,NUM)
    OPTION ARITHMETIC RATIONAL
    ASSIGN "miller.dll","isprime"
END FUNCTION

----------------------------------------------
素数判定高速です。

n=a*b*c
a,bが11桁を超えると低速になります。

http://blogs.yahoo.co.jp/donald_stinger

 

Re: ミラーラビン法

 投稿者:たろさ  投稿日:2017年 4月 5日(水)20時44分3秒
  > No.4085[元記事へ]

しばっちさんへのお返事です。


!ミラーラビン法 #4085
OPTION ARITHMETIC RATIONAL
LET k9=1000
DIM B(k9)
OPEN #2:NAME "E:\カーマイケル.TXT",ACCESS INPUT
FOR i=1 TO k9
   INPUT #2: B(i)
NEXT i
CLOSE #2
FOR I=1 TO k9
   LET A$=STR$(B(I))
   IF ISPRIME2(A$,100)=1 THEN
      PRINT A$
   END IF
NEXT  I
END

EXTERNAL FUNCTION ISPRIME(N$,NUM)
OPTION ARITHMETIC RATIONAL
ASSIGN "miller.dll","isprime"
END FUNCTION

EXTERNAL FUNCTION ISPRIME2(N$,NUM)
OPTION ARITHMETIC RATIONAL
ASSIGN "isprime.dll","isprime"
END FUNCTION

-----------------------------------------
計算結果

214852609因数分解229*457*2053
405739681因数分解229*457*3877
702683101因数分解421*701*2381
739444021因数分解277*1381*1933
775368901因数分解373*1117*1861
1213619761因数分解433*1297*2161
1312332001因数分解241*1361*4001
1648076041因数分解277*829*7177
1726372441因数分解307*1531*3673
1750412161因数分解241*2161*3361
2295209281因数分解313*1093*6709
2323147201因数分解277*1381*6073
2561945401因数分解421*701*8681

----------------------------------
IF ISPRIME2(A$,1000)=1 THEN
試しましたが、結果は同様でした。

IF ISPRIME(A$,100)=1 THEN
こちらは、問題ナシです。

http://blogs.yahoo.co.jp/donald_stinger

 

戻る