大整数的快速质因子分解,用到pollard-rho启发式算法。
算法导论上介绍pollard-rho介绍得比较详细,由于小因子的循环节长度很小,通过倍增步长,pollard-rho能够很快的找到一个大整数的一个较小的素因子p,书中说复杂度在O(sqrt(p))内,用的什么概率的分析方法,不懂。实际中pollard-rho的速度还是很快的,当然可能出现死循环。
这个题目要利用pollard-rho找到一个数的最小素因子,因此还需要Miller-Rabin测试来辅助。原来写的那个Miller-Rabin很快挂掉了,因为没有用到二次探测,判不出来Carmichael数。如果x ^ 2 = 1 (mod n),如果n是质数,那么x只能是1和n - 1;二次探测就是利用这个定理来进行检测。
POJ的论坛里面更有牛人列出了N多Carmichael数,真不知道他怎么找到的。最初怎么也不知道二次探测加在哪里好,后来参考网上一位大牛的代码,它的方法是计算a ^ b % n的时候,先将b折半到一个奇数b'为止,计算a ^ b',然后倍增b',同时进行二次检测,想想觉得很有道理,因为如果x ^ 2 = 1 mod n成立的话,那么(x ^ 2) ^ 2 = 1 mod n也成立。
这个题目还有一个trick就是模能达到2 ^ 54,如果这样计算一个数平方的时候,即使long long也会溢出。后来发现可以用快速幂取模的思想弄个"快速积取模",同样将b表示成二进制的形式,倍增的同时加到结果上就行了。这样每次运算的数范围都在2 ^ 54以内,并且都是加操作,不会溢出了。
POJ上还有一个用pollard-rho做的题是PKU 2429,这个比Prime Test还恶,因为这个题目是彻彻底底进行factorization,而且之后还要枚举一下找最优解,总之我的代码非常的长,而且这个题目数据范围2 ^ 63,必须用unsigned long long才能过。我的代码中间出现了一些减操作,都要特殊处理一下。最后还犯了个低级错误,函数返回值写错了,找了好几遍才找出来。
附PKU 1811代码:
PKU 1811
1 #include <cstdio>
2 #include <stdlib.h>
3 #include <time.h>
4 const int MAX = 8, N = 1000010;
5
6 long long ans, p[N/5];
7 int tag[N];
8
9 long long abs(long long t)
10 {
11 return t < 0 ? -t : t;
12 }
13 void prime()
14 {
15 int cnt = 0;
16
17 for (int i = 2; i < N; i++)
18 {
19 if (!tag[i]) p[cnt++] = i;
20 for (int j = 0; j < cnt && p[j] * i < N; j++)
21 {
22 tag[i*p[j]] = p[j];
23 if (i % p[j] == 0)
24 break;
25 }
26 }
27 }
28 long long gcd(long long a, long long b)
29 {
30 return !b ? a : gcd(b, a % b);
31 }
32 long long MultiMod(long long a, long long b, long long k)
33 {
34 long long ret = 0, tmp = a % k;
35
36 while (b)
37 {
38 if (b & 1)
39 ret = (ret + tmp) % k;
40 tmp = (tmp << 1) % k;
41 b >>= 1;
42 }
43
44 return ret;
45 }
46 long long PowerMod(long long a, long long b, long long k)
47 {
48 long long ret = 1, f = a % k;
49
50 while (b)
51 {
52 if (b & 1)
53 ret = MultiMod(ret, f, k);
54 f = MultiMod(f, f, k);
55 b >>= 1;
56 }
57 return ret;
58 }
59 int TwiceDetect(long long a, long long b, long long k)
60 {
61 int t = 0;
62 long long x, y;
63
64 while ((b & 1) == 0)
65 {
66 b >>= 1;
67 t++;
68 }
69 y = x = PowerMod(a, b, k);
70 while (t--)
71 {
72 y = MultiMod(x, x, k);
73 if (y == 1 && x != 1 && x != k - 1)
74 return 0;
75 x = y;
76 }
77 return y != 1 ? 0 : 1;
78 }
79 bool MillerRabin(long long n)
80 {
81 int i;
82 long long tmp;
83
84 srand(time(0));
85 for (i = 0; i < MAX; i++)
86 {
87 tmp = rand() % (n - 1) + 1;
88 if (TwiceDetect(tmp, n - 1, n) != 1)
89 break;
90 }
91 return (i == MAX);
92 }
93 long long pollard_rho(long long n)
94 {
95 srand(time(0));
96 int i = 1, k = 2;
97 long long d, x = rand() % (n - 1) + 1, y = x;
98 while (true)
99 {
100 i++;
101 x = (MultiMod(x, x, n) - 1) % n;
102 d = abs(gcd(y - x, n));
103 if (d != 1 && d != n)
104 return d;
105 if (i == k)
106 {
107 y = x;
108 k <<= 1;
109 }
110 }
111 }
112 void factorization(long long n)
113 {
114 long long tmp;
115
116 if (n < N)
117 {
118 ans <?= !tag[n] ? n : tag[n];
119 return;
120 }
121 if (MillerRabin(n))
122 {
123 ans <?= n;
124 return;
125 }
126 tmp = pollard_rho(n);
127 factorization(tmp);
128 factorization(n / tmp);
129 }
130
131 int main()
132 {
133 int T;
134 long long n;
135
136 scanf("%d", &T);
137 prime();
138 while (T--)
139 {
140 scanf("%I64d", &n);
141 if (n < N)
142 {
143 if (tag[n] == 0)
144 puts("Prime");
145 else
146 printf("%d\n", tag[n]);
147 continue;
148 }
149 if (MillerRabin(n))
150 puts("Prime");
151 else
152 {
153 ans = 0xfffffffffffffffLL;
154 factorization(n);
155 printf("%I64d\n", ans);
156 }
157 }
158
159 return 0;
160 }
161
posted on 2009-04-03 20:06
sdfond 阅读(710)
评论(0) 编辑 收藏 引用 所属分类:
Algorithm - Number Theory