注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

那一日,泪水打湿雪花冰冷的心

曾经的滋味,回忆在一次次的离合中,回眸时,那一刻,如涟漪般在一刹那融化

 
 
 

日志

 
 

POJ 1811  

2010-04-18 14:17:27|  分类: poj |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
素因数分解,在网上查询了一下,基本上都是用Miller_Rabin算法判断是否是素数,再用Pollard_Rho或squfof进行分解,Pollard_Rho在算法导论上有介绍,而squfof可以查看en.wikipedia.org/wiki/SQUFOF
#include<cstdio>
#include <cstdlib>
typedef unsigned long long u64;
const int MAX = 100;
const int maxt=100;
u64 a,c,x,len, dig ,limit;
int t,tt;
inline long long intsqrt(long long n)
{
long long temp, nHat = 0, b = 0x80000000, bshft = 31;
do
{
if (n >= (temp = (((nHat<<1)+b)<<bshft--)))
{
nHat += b;
n -= temp;
}
} while (b >>= 1);
return nHat;
}
u64 mod(u64 a, u64 b, u64 n) {
     if(! a) return 0;
     return (((a & dig) * b) % n + (mod(a >> len, b, n) << len) % n) % n;
}
u64 by(u64 a, u64 b, u64 n) {
      u64 p;
      p = 8, len = 61;
      while(p < n) {
           p <<= 4;
           len -= 4;
      }
      dig = ((limit / p) << 1) - 1;
      return mod(a, b, n);
}
u64 random() {
      u64 a;
      a = rand();
      a *= rand();
      a *= rand();
      a *= rand();
      return a;
}
u64 square_multiply(u64 x, u64 c, u64 n) {
      u64 z = 1;
      while(c) {
           if(c & 1) z = by(z, x, n);
           x = by(x, x, n);
           c >>= 1;
      }
      return z;
}
bool Miller_Rabin(u64 n) {
      if(n < 2) return false;
      if(n == 2) return true;
      if(! (n & 1)) return false;
      u64 i, j, k, m, a;
      m = n - 1;
      k = 0;
      while(m % 2 == 0) {
           m >>= 1; k ++;
      }
      for(i = 0; i < MAX; i ++) {
           a = square_multiply(random() % (n - 1) + 1, m, n);
           if(a == 1) continue;
           for(j = 0; j < k; j ++) {
                if(a == n - 1) break;
                a = by(a, a, n);
           }
           if(j == k) return false;
      }
      return true;
}
u64 gcd(u64 a, u64 b) {
      if(b == 0) return a;
      return gcd(b, a % b);
}
u64 squfof(u64 n,u64 k){
    tt=0;
    u64 p1=intsqrt(k*n),q1=1,q2=k*n-p1*p1,b0,q0,p0;
    u64 pp=p1,qq;
    do{
      p0=p1,q0=q1,q1=q2;    
      b0=(pp+p0)/q1;
      p1=b0*q1-p0;
      q2=q0+b0*(p0-p1);
      qq=intsqrt(q2);
      tt++;
      if(tt>maxt)
        return 1;
    }while(qq*qq!=q2);
    b0=(pp-p1)/qq;
    p1=b0*qq+p1;
    q1=qq;
    q2=(k*n-p1*p1)/q1;
    do{
      p0=p1,q0=q1,q1=q2;
      b0=(pp+p0)/q1;
      p1=b0*q1-p0;
      q2=q0+b0*(p0-p1);    
      tt++;
      if(tt>maxt)
        return 1;
    }while(p0!=p1);
    return gcd(n,p0);
}      
u64 prime(u64 a) {
  if(Miller_Rabin(a)) return a;
  u64 qq=intsqrt(a);
  if(qq*qq==a)
    return prime(qq);
  u64 k=1,aa,mm;
  while((aa=squfof(a,k))<=1)
    k++;
  mm=prime(aa);
  aa=prime(a/aa);
  if(mm>aa)
    mm=aa;
  return mm;
}
int main(){
    limit=1;
    limit<<=63;
    scanf("%d",&t);
    while(t--){
      scanf("%I64d",&x);  
      if(Miller_Rabin(x))
        printf("Prime\n");
      else printf("%I64d\n",prime(x));
    }
    return 0;
}
这是我写的squfof分解,另附上别人的Pollard_Rho分解
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
using namespace std;
typedef unsigned __int64 u64;
const int MAX = 100;
const int MAXN =100;
u64 len, dig ,limit;
u64 mod(u64 a, u64 b, u64 n) {
     if(! a) return 0;
     return (((a & dig) * b) % n + (mod(a >> len, b, n) << len) % n) % n;
}
u64 by(u64 a, u64 b, u64 n) {
      u64 p;
      p = 8, len = 61;
      while(p < n) {
           p <<= 4;
           len -= 4;
      }
      dig = ((limit / p) << 1) - 1;
      return mod(a, b, n);
}
u64 random() {
      u64 a;
      a = rand();
      a *= rand();
      a *= rand();
      a *= rand();
      return a;
}
u64 square_multiply(u64 x, u64 c, u64 n) {
      u64 z = 1;
      while(c) {
           if(c & 1) z = by(z, x, n);
           x = by(x, x, n);
           c >>= 1;
      }
      return z;
}
bool Miller_Rabin(u64 n) {
      if(n < 2) return false;
      if(n == 2) return true;
      if(! (n & 1)) return false;
      u64 i, j, k, m, a;
      m = n - 1;
      k = 0;
      while(m % 2 == 0) {
           m >>= 1; k ++;
      }
      for(i = 0; i < MAX; i ++) {
           a = square_multiply(random() % (n - 1) + 1, m, n);
           if(a == 1) continue;
           for(j = 0; j < k; j ++) {
                if(a == n - 1) break;
                a = by(a, a, n);
           }
           if(j == k) return false;
      }
      return true;
}
u64 gcd(u64 a, u64 b) {
      if(b == 0) return a;
      return gcd(b, a % b);
}
u64 f(u64 x, u64 n) {
      return (by(x, x, n) + 1) % n;
}
u64 Pollard_Rho(u64 n) {
      if(n <= 2) return 0;
      if(! (n & 1)) return 2;
      u64 i, p, x, xx;
      for(i = 1; i < MAX; i ++) {
           x = random() % n;
           xx = f(x, n);
           p = gcd((xx + n - x) % n, n);
           while(p == 1) {
                x = f(x, n);
                xx = f(f(xx, n), n);
                p = gcd((xx + n - x) % n, n) % n;
           }
           if(p) return p;
      }
      return 0;
}
u64 factor[MAXN], m;
u64 Prime(u64 a) {
      if(Miller_Rabin(a)) return 0;
      u64 t = Pollard_Rho(a);
      u64 p = Prime(t);
      if(p) return p;
      return t;
}
int main() {
#ifdef JUDGE_ONLINE
      freopen("input.txt", "r", stdin);
      //freopen("output.txt", "w", stdout);
#endif
      u64 a, t;
      limit = 1;
      limit = limit << 63;
      int casn;
      scanf("%d", &casn);
      while(casn --) {
           scanf("%I64d", &a);
           m = 0;
           while(a > 1) {
                if(Miller_Rabin(a)) break;
                t = Prime(a);
                factor[m ++] = t;
                a /= t;
           }
           if(a > 0) factor[m ++] = a;
           if(m == 1) printf("Prime\n"); //所有素约数保存于factor[]中,m为其个数
           else {
                int k = 0;
                for(int i = 1; i < m; i ++) {
                     if(factor[k] > factor[i])
                           k = i;
                }
                printf("%I64d\n", factor[k]);
           }
      }
      return(0);
}

  评论这张
 
阅读(1807)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018