2012-06-02 31 views
1

任何人都可以幫我用pollard rho實現嗎?我已經在C中實現了這個功能,它對於數字高達10位的數字工作正常,但無法處理更大的數字。Pollard Rho因子分解方法在C中的實現

請幫我改進它,以執行數字分解18位數字。我的代碼是this

#include<stdio.h> 
#include<math.h> 
int gcd(int a, int b) 
{ 
    if(b==0) return a ; 
    else 
    return(gcd(b,a%b)) ; 
} 

long long int mod(long long int a , long long int b , long long int n) 
{  
    long long int x=1 , y=a ; 
    while(b>0) 
    { 
     if(b%2==1) x = ((x%n)*(y%n))%n ; 
     y = ((y%n)*(y%n))%n ; 
     b/=2 ; 
    } 
    return x%n ; 
} 

int isprimes(long long int u) 
{ 
    if(u==3) 
    return 1 ; 
    int a = 2 , i ; 
    long long int k , t = 0 , r , p ; 
    k = u-1 ; 
    while(k%2==0) 
    { k/=2 ; t++ ; } 

     while(a<=3)                /*der are no strong pseudoprimes common in base 2 and base 3*/ 
     { 
     r = mod(a,k,u) ; 
      for(i = 1 ; i<=t ; i++) 
      { 
        p = ((r%u)*(r%u))%u ; 
        if((p==1)&&(r!=1)&&(r!=(u-1))) 
        { return 0 ; } 
        r = p ; 
      } 
      if(p!=1) 
      return 0 ; 
     else 
      a++ ; 
     } 

      if(a==4) 
      return 1 ; 

} 

long long int pol(long long int u) 
{ 
    long long int x = 2 , k , i , a , y , c , s; 
    int d = 1 ; 
    k = 2 ; 
    i = 1 ; 
    y = x ; 
    a = u ; 
    if(isprimes(u)==1) 
    { 
    return 1; 
    } 
    c=-1 ; 
    s = 2 ; 
    while(1) 
    { 
    i++; 
    x=((x%u)*(x%u)-1)% u ; 

    d = gcd(abs(y-x),u) ; 

    if(d!=1&&d!=u) 
    { printf("%d ",d); 
     while(a%d==0) { a=a/d; } 

     x = 2 ; 
     k = 2 ; 
     i = 1 ; 
     y = x ; 
     if(a==1) 
     { return 0 ; } 
     if(isprimes(a)!=0) 
     { return a ; } 
     u=a ; 

    } 
    if(i==k) 
    {y = x ; k*=2 ; c = x ;}              /*floyd cycle detection*/ 
     if(c==x)                 
    { x = ++s ; } 
    } 
    return ; 

} 

int main() 
{ 
    long long int t ; 
    long long int i , n , j , k , a , b , u ; 
    while(scanf("%lld",&n)&&n!=0) 
    { u = n ; k = 0 ; 
    while(u%2==0) 
     { u/=2 ; k = 1 ; } 
     if(k==1) printf("2 ") ; 
     if(u!=1) 
     t = pol(u) ; 
     if(u!=1) 
     { 
      if(t==1) 
      { printf("%lld",u) ; } 
      else 
      if(t!=0) 
      { printf("%lld",t) ; } 
     } 
      printf("\n"); 
    } 
    return 0; 
} 

遺憾的長碼.....我是一個新的編碼器。

+1

建議你寫一些單元測試... –

+1

你試過緩存素數? –

+3

該代碼讓我感到噁心...我知道這是有競爭力的編碼風格和所有,但請您在提出問題之前重寫一下它? – nhahtdh

回答

4

當您乘以兩個數字模m時,中間產品可能變得接近m^2。因此,如果使用64位無符號整數類型,則可以處理的最大模數爲2^32,如果模數較大,則可能會發生溢出。當模數僅稍大一些時,這種情況很少見,但這種情況只會不那麼明顯,如果模數允許溢出的可能性,則不能依靠幸運。

uint64_t mod_mul(uint64_t x, uint64_t y, uint64_t m) 
{ 
    int neg = 0; 
    // if x is too large, choose m-x and note that we need one negation for that at the end 
    if (x > m/2) { 
     x = m - x; 
     neg = !neg; 
    } 
    // if y is too large, choose m-y and note that we need one negation for that at the end 
    if (y > m/2) { 
     y = m - y; 
     neg = !neg; 
    } 
    uint64_t prod = (x * y) % m; 
    // if we had negated _one_ factor, and the product isn't 0 (mod m), negate 
    if (neg && prod) { 
     prod = m - prod; 
    } 
    return prod; 
} 

因此,這將允許多達模量:

您可以通過兩個因素,如果你在最m/2或東西相當於選擇剩餘類模m的絕對值的代表獲得更大的範圍以2^33與64位無符號類型。不是一大步。

問題的推薦解決方案是使用大整數庫,例如GMP可作爲大多數(如果不是全部)Linux發行版上的分發包,並且(相對)可在Windows上輕鬆安裝。

如果這不是一個選項(?真的,你確定),你可以得到它的更大的彈性模量的工作(高達2^63爲64位無符號整型)使用俄羅斯農民乘法:

x * y = 2 * (x * (y/2)) + (x * (y % 2)) 

所以爲了計算,您只需要2*(m-1)不會溢出。

uint64_t mod_mult(uint64_t x, uint64_t y, uint64_t m) 
{ 
    if (y == 0) return 0; 
    if (y == 1) return x % m; 
    uint64_t temp = mod_mult(x,y/2,m); 
    temp = (2*temp) % m; 
    if (y % 2 == 1) { 
     temp = (temp + x) % m; 
    } 
    return temp; 
} 

但請注意,該算法需要O(log y)步驟,所以在實踐中它相當慢。對於較小的m可以加快速度,如果2^k*(m-1)沒有溢出,則可以以k位而不是單個位(x*y = ((x * (y >> k)) << k) + (x * (y & ((1 << k)-1))))的步長進行,如果您的模數從未大於48位或56位,則這是一個很好的改進。

使用模乘算法的這種變體,你的算法將適用於更大的數字(但它會明顯更慢)。如果m < 2^32x < (2^64-1)/y,簡單(x * y) % m會做,您也可以嘗試測試模數和/或因素的大小以確定使用哪種方法。

+1

謝謝...這有助於很多:) – SlashGeek