大整数的快速质因子分解，用到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 @
2009-04-03 20:06 sdfond 阅读(693) |
评论 (0) |
编辑 收藏
题目大意是给定一个数n，问约数个数为n的最小的数k是多少。其中1 <= n <= 10000, k <= 10 ^ 15。
这是一个经典问题了，我一直以为会有经典算法，开始的时候一直往贪心上想，结果owen给出了反例。后来经过吉大牛点拨，因为k <= 10 ^ 15，可以根据这个定界，最差情况k的素因子也不会超过13，这样就可以搜索了！
实现的时候我也犯了几个小错，一个是把10 ^ 15少打了一个0，还有一个剪枝必须加：如果当前结果的约数个数为f，那么如果n % f不为0，则剪掉，因为约数个数是以乘积的关系累加的。
1 #include <cstdio>
2 const int M = 14;
3 const long long max = 1000000000000000LL;
4
5 int p[M] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41}, k;
6 long long ans;
7 void solve(long long v, int factor, int pos)
8 {
9 if (factor >= k)
10 {
11 if (factor == k) ans <?= v;
12 return;
13 }
14 if (k % factor) return;
15 if (pos == M) return;
16 for (int i = 1; i <= 50; i++)
17 {
18 v *= p[pos];
19 if (v > max) break;
20 solve(v, factor * (i + 1), pos + 1);
21 }
22 }
23
24 int main()
25 {
26 while (scanf("%d", &k) == 1)
27 {
28 ans = max + 1;
29 solve(1, 1, 0);
30 if (ans > max) printf("-1\n");
31 else printf("%lld\n", ans);
32 }
33
34 return 0;
35 }
36
posted @
2009-03-30 21:44 sdfond 阅读(277) |
评论 (0) |
编辑 收藏
给定一个数N（N <= 10 ^ 1000），如何快速求得N！的最末位非零数是一个经典的问题。一直以来都被这个问题困扰，今天仔细想了下，终于给想通了，尽管可能有些笨拙，现把想法记录于此。
在N很小的情况下，有一个简便的方法：求出1到N之间每个数的2的因子数和5的因子数，记为F(2)和F(5)，显然F(2) >= F(5)。由于在末尾只有2和5相乘才能产生0，如果我们把2和5抛去，那么肯定不会有0，这样就可以一边乘一边模10，防止溢出。剩下的一堆2和5如何处理呢？因为2肯定比5多，因此最末位肯定是偶数（0的阶乘和1的阶乘除外）。而一个偶数不停地乘2，最末位的规律是：2 -> 4 -> 8 -> 6 -> 2 -> ...出现了4位1循环，这样我们先用F(2) - F(5)，使得一部分2和5匹配上，2 * 5 = 10，对末尾不产生影响，剩下的2就模一下4，剩几再乘几次2就可以了。
但是这个方法在N非常大的时候肯定就不行了，但是可以利用找循环这个思想继续做。如果算阶乘的时候跳过5的倍数，记G(n)为跳过5的倍数的时候，从1乘到n的最末非零位，也就是把5的倍数当1乘。可以发现：
G(1) = 1, G(2) = 2, G(3) = 6, G(4) = 4, G(5) = 4, G(6) = 4, G(7) = 8, G(8) = 4, G(9) = 6, G(10) = 6, G(11) = 6, G(12) = 2, G(13) = 6...
又出现了循环，每10个数循环一次。如何计算G(n)就变的很简单，求出n的最末位，就知道对应的G(n)是多少了，当然需要特判n = 1的情况。由于我们把5的倍数的数都提出来了，提出来的这些数（5、10、15、20、25、30...）每个除以5后又组成了一个阶乘序列！除完5一共提出了n / 5个5，根据之前的分析，每个5都可以拿出一个2和它配对然后把它消去，这样一个5就相当于少一个2，我们就要把原来的数乘以3个2（模四循环）。这样一来5的个数其实也可以模四，模完四之后剩k的话，就可以乘以k个8，就把所有的5消去了。现在总结一下：对一个数n的阶乘，计算它的末尾非零位，先计算G(n)，相当于非5的倍数的数的乘积最末非零位先算好了，然后乘以n / 5 % 4个8，处理了提出的n / 5个5，这样之后还剩下n / 5的阶乘没有算。递归的求解n / 5的阶乘的最末位非零数，再乘上去就得到结果了。
这个做法的复杂度就很低了，达到O(log n)，对于10 ^ 1000的数据，利用高精度做就行了。利用这种循环的思想，算排列数P(n, k)的最末非零数也就可以做到了。
附HOJ 1013代码：
1
#include <cstdio>
2#include <cstring>
3const int N = 1024;
4
5int hash[10] = 
{6, 6, 2, 6, 4, 4, 4, 8, 4, 6};
6int one_digit_hash[10] = {1, 1, 2, 6, 4, 2, 2, 4, 2, 8};
7
8int last_digit(char str[N], int st, int to)
9{
10
int i, tmp = 0, ret, num_of_five = 0;
11
12 if (st == to)
13 return one_digit_hash[str[st]-'0'];
14
15 ret = hash[str[to]-'0'];
16 for (i = st; i <= to; i++)
17
{
18 tmp = tmp * 10 + str[i] - '0';
19 str[i] = tmp / 5 + '0';
20 tmp %= 5;
21 num_of_five = (num_of_five * 10 + str[i] - '0') % 4;
22
}
23 if (str[st] == '0') st++;
24 ret = last_digit(str, st, to) * ret % 10;
25 while (num_of_five--) ret = ret * 8 % 10; //mul one 5 equals mul one 8
26
27 return ret;
28
}
29
30int main()
31{
32 char str[N];
33
34 while (scanf("%s", str) == 1)
35 printf("%d\n", last_digit(str, 0, strlen(str) - 1));
36
37 return 0;
38}
39
posted @
2009-03-29 21:00 sdfond 阅读(647) |
评论 (2) |
编辑 收藏
摘要: 一道思路很简单的计算几何题目，就是先判是不是“凸多边形”，然后计算点到直线的最短距离。但是我错了很多次。一个问题就是有可能peg不在多边形内，这要单独判断一下；还有一个问题找了很久才发现，原来题目中说满足条件的多边形不是纯粹的凸多边形，题目中的多边形是“任意内部两点连线不会和多边形的边相交”，这样如果多边形的多个顶点存在共线的情况，其实也是可以的，但...
阅读全文
posted @
2009-03-29 11:58 sdfond 阅读(291) |
评论 (0) |
编辑 收藏最近又研究了线性筛素数方法的扩展，的确非常强大。
经典的应用是线性时间筛欧拉函数和求约数个数。我想了一下线性时间筛约数和(积性函数)，也是可行的。
对于一个大于1的数n，可以写成p1 ^ a1 * p2 ^ a2 * ... * pn ^ an，那么n的约数和就是(p1 ^ 0 + p1 ^ 1 + ... p1 ^ a1) * ... * (pn ^ 0 + ... + pn ^ an)，由此就可以有递推关系了：
设f[i]表示i的约数和，e[i]表示i的最小素因子个数，t[i]表示(p1 ^ 0 + .. + p1 ^ a1)，p1是t的最小素因子，a1是p1的幂次，这样对于i * p[j]，如果p[j]不是i的因子，那么根据积性条件，f[i*p[j]] = f[i] * (1 + p[j])，e[i] = 1，t[i] = 1 + p[j]；如果p[j]是i的因子，那么相当于t[i]多了一项p1 ^ (a1 + 1)，首先e[i]++，然后tmp = t[i]，t[i] += p[j] ^ e[i]，f[i*p[j]] = f[i] / tmp * t[i]。这样也就做到了O(1)的时间计算出了f[i*p[j]]同时也计算出了附加信息。
这种方法还可以继续推广，例如可以记录i的最小素因子，这样就可以做到O(log n)时间的素因子分解。
posted @
2009-03-28 15:02 sdfond 阅读(183) |
评论 (0) |
编辑 收藏