时间戳

首页 > 解决方案 > 在 C 中使用 GMP 库测试素数

c - 在 C 中使用 GMP 库测试素数

问题描述

这是我第一次使用 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;
}

请帮帮我。

标签： cgmp

解决方案

第一的; 您可以/应该在其中添加额外的打印（例如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 -bit piece, ...). 为了最大限度地提高效率（在产品不适合“数字”之前最大限度地增加表中的条目数，并且我必须减少表条目中的除数）我生成这些表使用“smalls and bigs”方法（例如，19 和 137 很小，421 和 587 很大，19*137*421*587 = 0x26578B9D = 一个小到足以容纳 32 位的数字）。

  • 您不需要大整数 (mpz) 来跟踪 Miller Rabin 轮数（while(mpz_cmp_ui(rep,i)>=0)循环）。在开始循环之前将您的转换rep为 aint并while(i <= int_rep改为执行。请注意（对于寻找 RSA-4096 的素数之类的事情）256 轮米勒拉宾是不必要的过度杀伤。

  • 米勒拉宾尴尬地平行。如果您有 8 个 CPU，那么您想使用 8 个线程，并希望并行执行（最多）8 轮 Miller Rabin。注意：如果你没有做足够的试验划分，那么米勒拉宾的第一次迭代会经常说“复合”，在这种情况下你可能想用一个线程做米勒拉宾的第一次迭代（并且只使用多个线程第一次迭代后）。

推荐阅读