大整數(shù)的快速質(zhì)因子分解，用到pollard-rho啟發(fā)式算法。
算法導(dǎo)論上介紹pollard-rho介紹得比較詳細(xì)，由于小因子的循環(huán)節(jié)長度很小，通過倍增步長，pollard-rho能夠很快的找到一個(gè)大整數(shù)的一個(gè)較小的素因子p，書中說復(fù)雜度在O(sqrt(p))內(nèi)，用的什么概率的分析方法，不懂。實(shí)際中pollard-rho的速度還是很快的，當(dāng)然可能出現(xiàn)死循環(huán)。
這個(gè)題目要利用pollard-rho找到一個(gè)數(shù)的最小素因子，因此還需要Miller-Rabin測試來輔助。原來寫的那個(gè)Miller-Rabin很快掛掉了，因?yàn)闆]有用到二次探測，判不出來Carmichael數(shù)。如果x ^ 2 = 1 (mod n)，如果n是質(zhì)數(shù)，那么x只能是1和n - 1；二次探測就是利用這個(gè)定理來進(jìn)行檢測。
POJ的論壇里面更有牛人列出了N多Carmichael數(shù)，真不知道他怎么找到的。最初怎么也不知道二次探測加在哪里好，后來參考網(wǎng)上一位大牛的代碼，它的方法是計(jì)算a ^ b % n的時(shí)候，先將b折半到一個(gè)奇數(shù)b'為止，計(jì)算a ^ b'，然后倍增b'，同時(shí)進(jìn)行二次檢測，想想覺得很有道理，因?yàn)槿绻鹸 ^ 2 = 1 mod n成立的話，那么(x ^ 2) ^ 2 = 1 mod n也成立。
這個(gè)題目還有一個(gè)trick就是模能達(dá)到2 ^ 54，如果這樣計(jì)算一個(gè)數(shù)平方的時(shí)候，即使long long也會溢出。后來發(fā)現(xiàn)可以用快速冪取模的思想弄個(gè)"快速積取模"，同樣將b表示成二進(jìn)制的形式，倍增的同時(shí)加到結(jié)果上就行了。這樣每次運(yùn)算的數(shù)范圍都在2 ^ 54以內(nèi)，并且都是加操作，不會溢出了。
POJ上還有一個(gè)用pollard-rho做的題是PKU 2429，這個(gè)比Prime Test還惡，因?yàn)檫@個(gè)題目是徹徹底底進(jìn)行factorization，而且之后還要枚舉一下找最優(yōu)解，總之我的代碼非常的長，而且這個(gè)題目數(shù)據(jù)范圍2 ^ 63，必須用unsigned long long才能過。我的代碼中間出現(xiàn)了一些減操作，都要特殊處理一下。最后還犯了個(gè)低級錯(cuò)誤，函數(shù)返回值寫錯(cuò)了，找了好幾遍才找出來。
附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 閱讀(723) |
評論 (0) |
編輯 收藏
題目大意是給定一個(gè)數(shù)n，問約數(shù)個(gè)數(shù)為n的最小的數(shù)k是多少。其中1 <= n <= 10000, k <= 10 ^ 15。
這是一個(gè)經(jīng)典問題了，我一直以為會有經(jīng)典算法，開始的時(shí)候一直往貪心上想，結(jié)果owen給出了反例。后來經(jīng)過吉大牛點(diǎn)撥，因?yàn)閗 <= 10 ^ 15，可以根據(jù)這個(gè)定界，最差情況k的素因子也不會超過13，這樣就可以搜索了！
實(shí)現(xiàn)的時(shí)候我也犯了幾個(gè)小錯(cuò)，一個(gè)是把10 ^ 15少打了一個(gè)0，還有一個(gè)剪枝必須加：如果當(dāng)前結(jié)果的約數(shù)個(gè)數(shù)為f，那么如果n % f不為0，則剪掉，因?yàn)榧s數(shù)個(gè)數(shù)是以乘積的關(guān)系累加的。
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 閱讀(304) |
評論 (0) |
編輯 收藏
給定一個(gè)數(shù)N（N <= 10 ^ 1000），如何快速求得N！的最末位非零數(shù)是一個(gè)經(jīng)典的問題。一直以來都被這個(gè)問題困擾，今天仔細(xì)想了下，終于給想通了，盡管可能有些笨拙，現(xiàn)把想法記錄于此。
在N很小的情況下，有一個(gè)簡便的方法：求出1到N之間每個(gè)數(shù)的2的因子數(shù)和5的因子數(shù)，記為F(2)和F(5)，顯然F(2) >= F(5)。由于在末尾只有2和5相乘才能產(chǎn)生0，如果我們把2和5拋去，那么肯定不會有0，這樣就可以一邊乘一邊模10，防止溢出。剩下的一堆2和5如何處理呢？因?yàn)?肯定比5多，因此最末位肯定是偶數(shù)（0的階乘和1的階乘除外）。而一個(gè)偶數(shù)不停地乘2，最末位的規(guī)律是：2 -> 4 -> 8 -> 6 -> 2 -> ...出現(xiàn)了4位1循環(huán)，這樣我們先用F(2) - F(5)，使得一部分2和5匹配上，2 * 5 = 10，對末尾不產(chǎn)生影響，剩下的2就模一下4，剩幾再乘幾次2就可以了。
但是這個(gè)方法在N非常大的時(shí)候肯定就不行了，但是可以利用找循環(huán)這個(gè)思想繼續(xù)做。如果算階乘的時(shí)候跳過5的倍數(shù)，記G(n)為跳過5的倍數(shù)的時(shí)候，從1乘到n的最末非零位，也就是把5的倍數(shù)當(dāng)1乘?？梢园l(fā)現(xiàn)：
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...
又出現(xiàn)了循環(huán)，每10個(gè)數(shù)循環(huán)一次。如何計(jì)算G(n)就變的很簡單，求出n的最末位，就知道對應(yīng)的G(n)是多少了，當(dāng)然需要特判n = 1的情況。由于我們把5的倍數(shù)的數(shù)都提出來了，提出來的這些數(shù)（5、10、15、20、25、30...）每個(gè)除以5后又組成了一個(gè)階乘序列！除完5一共提出了n / 5個(gè)5，根據(jù)之前的分析，每個(gè)5都可以拿出一個(gè)2和它配對然后把它消去，這樣一個(gè)5就相當(dāng)于少一個(gè)2，我們就要把原來的數(shù)乘以3個(gè)2（模四循環(huán)）。這樣一來5的個(gè)數(shù)其實(shí)也可以模四，模完四之后剩k的話，就可以乘以k個(gè)8，就把所有的5消去了?，F(xiàn)在總結(jié)一下：對一個(gè)數(shù)n的階乘，計(jì)算它的末尾非零位，先計(jì)算G(n)，相當(dāng)于非5的倍數(shù)的數(shù)的乘積最末非零位先算好了，然后乘以n / 5 % 4個(gè)8，處理了提出的n / 5個(gè)5，這樣之后還剩下n / 5的階乘沒有算。遞歸的求解n / 5的階乘的最末位非零數(shù)，再乘上去就得到結(jié)果了。
這個(gè)做法的復(fù)雜度就很低了，達(dá)到O(log n)，對于10 ^ 1000的數(shù)據(jù)，利用高精度做就行了。利用這種循環(huán)的思想，算排列數(shù)P(n, k)的最末非零數(shù)也就可以做到了。
附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 閱讀(682) |
評論 (2) |
編輯 收藏
摘要: 一道思路很簡單的計(jì)算幾何題目，就是先判是不是“凸多邊形”，然后計(jì)算點(diǎn)到直線的最短距離。但是我錯(cuò)了很多次。一個(gè)問題就是有可能peg不在多邊形內(nèi)，這要單獨(dú)判斷一下；還有一個(gè)問題找了很久才發(fā)現(xiàn)，原來題目中說滿足條件的多邊形不是純粹的凸多邊形，題目中的多邊形是“任意內(nèi)部兩點(diǎn)連線不會和多邊形的邊相交”，這樣如果多邊形的多個(gè)頂點(diǎn)存在共線的情況，其實(shí)也是可以的，但...
閱讀全文
posted @
2009-03-29 11:58 sdfond 閱讀(320) |
評論 (0) |
編輯 收藏最近又研究了線性篩素?cái)?shù)方法的擴(kuò)展，的確非常強(qiáng)大。
經(jīng)典的應(yīng)用是線性時(shí)間篩歐拉函數(shù)和求約數(shù)個(gè)數(shù)。我想了一下線性時(shí)間篩約數(shù)和(積性函數(shù))，也是可行的。
對于一個(gè)大于1的數(shù)n，可以寫成p1 ^ a1 * p2 ^ a2 * ... * pn ^ an，那么n的約數(shù)和就是(p1 ^ 0 + p1 ^ 1 + ... p1 ^ a1) * ... * (pn ^ 0 + ... + pn ^ an)，由此就可以有遞推關(guān)系了：
設(shè)f[i]表示i的約數(shù)和，e[i]表示i的最小素因子個(gè)數(shù)，t[i]表示(p1 ^ 0 + .. + p1 ^ a1)，p1是t的最小素因子，a1是p1的冪次，這樣對于i * p[j]，如果p[j]不是i的因子，那么根據(jù)積性條件，f[i*p[j]] = f[i] * (1 + p[j])，e[i] = 1，t[i] = 1 + p[j]；如果p[j]是i的因子，那么相當(dāng)于t[i]多了一項(xiàng)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)的時(shí)間計(jì)算出了f[i*p[j]]同時(shí)也計(jì)算出了附加信息。
這種方法還可以繼續(xù)推廣，例如可以記錄i的最小素因子，這樣就可以做到O(log n)時(shí)間的素因子分解。
posted @
2009-03-28 15:02 sdfond 閱讀(212) |
評論 (0) |
編輯 收藏