这是我第一次使用gmp库，我想通过Fermat测试和Rabin Miller测试数字是否为质数。该程序可以很好地处理20位数的数字，但是我想测试300和500位数甚至更大的非常大的数字，并且测试非常大的数字会出错，我不知道为什么！
这里是代码
#include <gmp.h>
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
void S_and_M(mpz_t a,mpz_t n,mpz_t h, mpz_t r) // square and multiply
{
char * bin = mpz_get_str(NULL,2,h);
int i;
mpz_set(r,a);
for(i = 1; i < strlen(bin);i++)
{
mpz_mul(r,r,r);
mpz_mod(r,r,n);
if(bin[i] == '1')
{
mpz_mul(r,r,a);
mpz_mod(r,r,n);
}
}
}
void M2Pow(mpz_t s,mpz_t S_pow)
{
int i = 0;
mpz_set_ui(S_pow,1);
while(mpz_cmp_si(s,i) > 0 )
{
mpz_mul_ui(S_pow,S_pow,2);
i++;
}
}
void Decomp(mpz_t x,mpz_t s,mpz_t t) // Function to decompose x in 2^s * t
{
mpz_t y,S_pow;
mpz_init(y);
mpz_set_ui(y,0);
mpz_init(S_pow);
while(mpz_cmp(y,x)!=0) // while we don't find 2^s * t = x
{
mpz_set_ui(t,1); //we restart with t = 1
mpz_mul(y,S_pow,t);// y = 2^s * t (we test the values)
while(mpz_cmp(y,x) < 0)// we stop when 2^s * t > x
{
mpz_add_ui(t,t,2);
M2Pow(s,S_pow);
mpz_mul(y,S_pow,t);
}
mpz_add_ui(s,s,1);
}
mpz_sub_ui(s,s,1);
mpz_clear(y);
mpz_clear(S_pow);
}
void testFermat(mpz_t n, mpz_t rep)
{
gmp_randstate_t state;
gmp_randinit_mt(state);
gmp_randseed_ui(state, time(NULL));
mpz_t i;
mpz_init(i);
mpz_set_si(i,1);
mpz_t n2;
mpz_init(n2);
mpz_sub_ui(n2,n,1);
mpz_t a;
mpz_init(a);
mpz_t r;
mpz_init(r);
mpz_t n3;
mpz_init(n3);
mpz_sub_ui(n3,n,3);
while(mpz_cmp(i,rep)<=0 && mpz_cmp_si(n,2)!= 0 && mpz_cmp_si(n,3)!=0)
{
mpz_urandomm(a,state,n3);
mpz_add_ui(a,a,2);
S_and_M(a,n,n2,r);
if(mpz_cmp_si(r,1)!=0)
{
printf("The number is composite \n");
mpz_clear(i);mpz_clear(n2);
mpz_clear(a);mpz_clear(r);
mpz_clear(n3);gmp_randclear(state);
return ;
}
mpz_add_ui(i,i,1);
}
printf("The number is prime \n");
mpz_clear(i);mpz_clear(n2);
mpz_clear(a);mpz_clear(r);
mpz_clear(n3);gmp_randclear(state);
}
void Miller_Rabin(mpz_t n, mpz_t rep)
{
if(mpz_get_ui(n) % 2 == 0)
{
if(mpz_cmp_ui(n,2) == 0)
printf("The number is prime \n");
else
printf("The number is composite \n");
return;
}
int i=1;
mpz_t a,y,s,t,n1,n2,deux;
gmp_randstate_t state;
gmp_randinit_mt(state);
gmp_randseed_ui(state, time(NULL));
mpz_init(a);
mpz_init(y);
mpz_init(s);
mpz_set_ui(s,1);
mpz_init(deux);
mpz_set_ui(deux,2);
mpz_init(t);
mpz_init(n1);
mpz_sub_ui(n1,n,1);
mpz_init(n2);
mpz_sub_ui(n2,n,2);
Decomp(n1,s,t);
mpz_sub_ui(s,s,1);
while(mpz_cmp_ui(rep,i)>=0)
{
mpz_urandomm(a,state,n1);
mpz_add_ui(a,a,1);
S_and_M(a,n,t,y);
if(mpz_cmp_si(y,1)!=0 && mpz_cmp(y,n1)!=0)
{
for(int j=1;mpz_cmp_ui(s,j)>=0;j++)
{
mpz_set(n2,y);
S_and_M(y,n,deux,y);
if(mpz_cmp_si(y,1)==0)
{
printf("The number is composite\n");
mpz_clear(a);mpz_clear(y);
mpz_clear(s);mpz_clear(t);
mpz_clear(n1);mpz_clear(n2);
mpz_clear(deux);gmp_randclear(state);
return;
}
if(mpz_cmp(y,n1)==0) //Si y congrue a -1 mod n on sort de la boucle
break;
}
if(mpz_cmp(y,n1)!=0)
{
printf("The number is composiste \n");
mpz_clear(a);mpz_clear(y);
mpz_clear(s);mpz_clear(t);
mpz_clear(n1);mpz_clear(n2);
mpz_clear(deux);gmp_randclear(state);
return;
}
}
i++;
}
printf("The number is prime \n");
mpz_clear(a);mpz_clear(y);
mpz_clear(s);mpz_clear(t);
mpz_clear(n1);mpz_clear(n2);
mpz_clear(deux);gmp_randclear(state);
}
int main()
{
int t=1;
mpz_t n;
mpz_init(n);
mpz_t rep;
mpz_init(rep);
printf("########## Primality test ##########\n");
while(t==1)
{
printf("\n");
printf("Choose an integer to test : ");
gmp_scanf("%Zd", &n);
if(mpz_cmp_ui(n,1)<=0)
{
printf("\n choose an integer bigger than 1 !");
}
else
{
printf("Choose the number of repetitions : ");
gmp_scanf("%Zd", &rep);
printf("#########################################################\n");
printf("Miller_Rabin : ");
Miller_Rabin(n,rep);
printf("\n");
printf("Fermat : ");
testFermat(n,rep);
printf("#########################################################\n");
printf("tape 1 for an other test !");
scanf("%d",&t);
}
}
mpz_clear(n);
mpz_clear(rep);
return 0;
}
请帮助我。
首先；您可以/应该在其中进行额外的打印（例如printf(" Round %d\n", i);中的“ Miller_Rabin()”），以确定代码是否仍在工作以及它在做什么。
一些性能提示（假设您的代码花费了很长时间，以至于您假设它没有锁定时就锁定了它：]
将花费大量时间来尝试找到模数（例如，在S_and_M()中）。因为除数总是相同的，所以可以通过一次找到倒数（以足够的精度）并执行“ modulo = X - (X * reciprocal); if(modulo => divisor) { modulo -= divisor}; return modulo;”（而不是使用mpz_mod()）来加快除数。注意：仅使用整数，您实际上想执行“ shifted_reciprocal = (1 << bits) / divisor”来获得已向左移位bits位的倒数；然后执行“ modulo = X - ((X * shifted_reciprocal) >> bits);”。
当乘以大数时，您最终将第一个数字的每个“数字”（字，32位，...）与第二个数字的每个“数字”相乘。对于平方，两个数字具有相同的数字，因此可以对（乘以数字的）中间结果进行回收，从而使这些中间乘法的成本减半。这意味着设计为平方的代码几乎可以快于通用乘法（例如mpz_mul()）的两倍。
我不知道您的Decomp()在做什么（我怀疑您写它是想s需要处理大于计算机RAM位数的数字，这是荒谬的）。通常，您要计算“最低有效清除位”的数量（以找到s），然后进行一次右移以找到t（例如“ t = X >> s”）。为此，请查看GMP库的文档（我认为整个Decomp()可以/应该只是“ s = mpz_scan1(x, 0); mpz_tdiv_q_2exp(t, x, s);”）。
使用更多的审判庭。对于大数，如果除数小于“数字”（字，32位，...），则可以相对快速地找到模数。同样，您可以将多个小除数相乘，将它们相乘成一个“大数模”，然后使用结果的较小（整数）模。例如，可以不用“ if( (mpz_mod_ui(dummy, big_number, 3) == 0) || (mpz_mod_ui(dummy, big_number, 5) == 0) || (mpz_mod_ui(dummy, big_number, 7) == 0) || (mpz_mod_ui(dummy, big_number, 11) == 0)”，而可以用“ temp = mpz_mod_ui(dummy, big_number, 3*5*7*11); if( (temp%3 == 0) || (temp%5 == 0) || (temp%7 == 0) || (temp%11 == 0) )”。这明显更快。通过具有（预计算的）“除数除法器组”表，可以进一步利用此功能，以查找除法器的最大数目，这些除数相乘后将适合“数字”（字，32位，... ）。为了最大程度地提高效率（在产品不能放入“数字”之前最大化表的条目数，并且我必须减少表条目的除数的数量），我使用“大小”方法生成这些表（例如19和137较小，421和587大，并且19 * 137 * 421 * 587 = 0x26578B9D =一个足够小的数字以适合32位）。
您不需要大整数（mpz）即可跟踪Miller Rabin轮数（while(mpz_cmp_ui(rep,i)>=0)循环）。在开始循环之前，将rep转换为int，然后执行while(i <= int_rep。请注意（对于查找RSA-4096的质数来说，）不必要的256轮Miller Rabin就是多余的。]]
米勒·拉宾尴尬地平行。如果您有8个CPU，那么您想使用8个线程，并希望并行执行（最多）8轮Miller Rabin。注意：如果没有进行足够的试验划分，则Miller Rabin的第一次迭代将经常说“复合”，在这种情况下，您可能希望使用单个线程（仅使用多个线程）进行Miller Rabin的第一次迭代。第一次迭代后）。