1、HDUOJ 4675
题意：给定数n, m, k和数列{an}（1 <= ai <= m），求改变数列{an}中的k个数的值（不改变位置，且新值范围为1 <= ai' <= m，ai' != ai），使得数列的最大公约数为d的方案数为f[d]，输出f[1], f[2], f[3]....f[m]。
解法：最大公约数不好考虑，考虑将最大公约数转化为公约数。记满足题意且公约数为d的数列改变方案数为g[d]。
关键：g[d] = f[d] + f[2*d] + f[3*d] +... f[[m/d] * d]
下面考虑求g[d]，然后用消元的方法求f[d]。
对于一个固定的d，假设题目给定数列有 cnt 个 ai 为 d的倍数，且记t = m / d，则有：
1、cnt < n - k, 则g[d] = 0
2、cnt >= n - k，则将改变以后得到的新数列{bn}分为三种数，第一种，未改变的(bi = ai)，C(n-k, x)；第二种，ai 不是 d 的倍数，t ^ (n - cnt)；第三种，ai 是 d的倍数但bi != ai，(t - 1) ^ (cnt - n + k)。
tag：math, number theory, 计数, 欧拉公式, 乘法逆, good
1 /* 2 * Author: plum rain 3 * Created Time: 2013-08-16 11:10 4 * File Name: fuck.cpp 5 */ 6 #include<iostream> 7 #include<sstream> 8 #include<fstream> 9 #include<vector> 10 #include<list> 11 #include<deque> 12 #include<queue> 13 #include<stack> 14 #include<map> 15 #include<set> 16 #include<bitset> 17 #include<algorithm> 18 #include<cstdio> 19 #include<cstdlib> 20 #include<cstring> 21 #include<cctype> 22 #include<cmath> 23 #include<ctime> 24 #include<utility> 25 26 using namespace std; 27 28 #define CLR(x) memset(x, 0, sizeof(x)) 29 #define PB push_back 30 #define SZ(v) ((int)(v).size()) 31 #define INF 999999999999 32 #define zero(x) (((x)>0?(x):-(x))<eps) 33 #define out(x) cout<<#x<<":"<<(x)<<endl 34 #define tst(a) cout<<#a<<endl 35 #define CINBEQUICKER std::ios::sync_with_stdio(false) 36 37 typedef vector<int> VI; 38 typedef vector<string> VS; 39 typedef vector<double> VD; 40 typedef long long int64; 41 42 const double eps = 1e-8; 43 const double PI = atan(1.0)*4; 44 const int maxint = 2139062143; 45 const int N = 300005; 46 const int mod = 1000000007; 47 48 inline int Mymod( int a , int b ) { int x=a%b;if(x<0) x+=b;return x;} 49 50 int n, m, k; 51 int an[N]; 52 int64 f[N]; 53 int64 C[N]; 54 55 int64 Mypow(int64 a, int64 n) 56 { 57 int64 ret = 1, tmp = a % mod; 58 while (n){ 59 if (n&1) 60 ret = (ret * tmp) % mod; 61 tmp = (tmp * tmp) % mod; 62 n >>= 1; 63 } 64 return ret; 65 } 66 67 void Cinit() 68 { 69 CLR (C); 70 int64 x = n - k; 71 C[x] = 1; 72 for (int i = x; i < n; ++ i){ 73 C[i+1] = (C[i] * (i + 1)) % mod; 74 C[i+1] = (C[i+1] * Mypow(i+1-x, mod-2)) % mod; 75 } 76 } 77 78 79 void solve(int x) 80 { 81 int cnt = 0; 82 int64 sum = 0; 83 for (int i = 1; i*x <= m; ++ i){ 84 cnt += an[i*x]; 85 if (i > 1) sum = (sum + f[i*x]) % mod; 86 } 87 int t = m / x; 88 if (t == 1){ 89 if (cnt == n - k) f[x] = 1; 90 else f[x] = 0; 91 return ; 92 } 93 if (cnt < n - k){ 94 f[x] = 0; 95 return; 96 } 97 int64 tmp = (C[cnt] * Mypow(t, n-cnt)) % mod; 98 tmp = (tmp * Mypow(t-1, cnt - n + k)) % mod; 99 f[x] = (tmp - sum + mod) % mod; 100 } 101 102 int main() 103 { 104 // freopen("tst.in","r",stdin); 105 // freopen("fuck.out","w",stdout); 106 // std::ios::sync_with_stdio(false); 107 while (scanf ("%d%d%d", &n, &m, &k) != EOF){ 108 CLR (an); 109 int x; 110 for (int i = 0; i < n; ++ i){ 111 scanf ("%d", &x); 112 ++ an[x]; 113 } 114 115 CLR (f); 116 Cinit(); 117 for (int i = m; i >= 1; -- i) 118 solve(i); 119 120 printf ("%I64d", f[1]); 121 for (int i = 2; i <= m; ++ i) 122 printf (" %I64d", f[i]); 123 printf ("\n"); 124 } 125 return 0; 126 }
2、World Finals 2008, LA 4119
题意：给定一个关于正整数n的多项式P(n)和一个整数D，判断P(n) / D，是不是恒为整数。
解法：要运用差分数列的思想，数学归纳法思想。首先给出结论，设多项式的最高次项次数为k，则如果n = 1, 2, 3...n+1时P(n)都是D的倍数，那么P(n) / D恒为整数。
首先k = 0, 则只需要验证P(1)即可。
k = 1，设P(n) = a*n + b，P(n)为等差数列，已验证P(1), P(2)为D的倍数，则有数列中首项为D的倍数，公差d = P(2) - P(1)也为D的倍数，所以每一项都是D的倍数。
k = 2，设P(n) = a*n^2 + b*n + c，则P(m+1) - P（m) = 2*a*m + a + b，为一次多项式。已验证P(1)，P(2)，P(3)为D的倍数，则P(3) - P(2)，P(2) - P(1)，由上面证明的结论可以得到对任意m有P(m+1) - P(m)为D的倍数且P(1)为D的倍数，所以P(m)也为D的倍数，即P(n)每一项都为D的倍数。
以下可以循环论证了....
tag:math, number theory, 差分数列, good
1 /* 2 * Author: plum rain 3 * Created Time: 2013-08-23 21:00 4 * File Name: (good)-math-WF2008-LA-4119.cpp 5 */ 6 #include<iostream> 7 #include<sstream> 8 #include<fstream> 9 #include<vector> 10 #include<list> 11 #include<deque> 12 #include<queue> 13 #include<stack> 14 #include<map> 15 #include<set> 16 #include<bitset> 17 #include<algorithm> 18 #include<cstdio> 19 #include<cstdlib> 20 #include<cstring> 21 #include<cctype> 22 #include<cmath> 23 #include<ctime> 24 #include<utility> 25 26 using namespace std; 27 28 #define CLR(x) memset(x, 0, sizeof(x)) 29 #define PB push_back 30 #define SZ(v) ((int)(v).size()) 31 #define INF 999999999999 32 #define zero(x) (((x)>0?(x):-(x))<eps) 33 #define out(x) cout<<#x<<":"<<(x)<<endl 34 #define tst(a) cout<<#a<<endl 35 #define CINBEQUICKER std::ios::sync_with_stdio(false) 36 37 typedef vector<int> VI; 38 typedef vector<string> VS; 39 typedef vector<double> VD; 40 typedef long long int64; 41 42 const double eps = 1e-8; 43 const double PI = atan(1.0)*4; 44 const int maxint = 2139062143; 45 46 inline int Mymod( int a , int b ) { int x=a%b;if(x<0) x+=b;return x;} 47 48 struct P{ 49 int k, a, b; 50 //k * n ^ a 51 void clr(){ 52 k = 0; a = 0; 53 } 54 }; 55 56 P a[105]; 57 int tot; 58 string s; 59 int mod; 60 61 int64 Mypow(int64 p, int64 n) 62 { 63 int64 ret = 1; 64 while (n > 0){ 65 if (n & 1) 66 ret = (ret * p) % mod; 67 n >>= 1; 68 p = (p * p) % mod; 69 } 70 return ret; 71 } 72 73 void init() 74 { 75 int size = SZ (s); 76 int pos = 1; 77 for (int i = 0; i < 105; ++ i) 78 a[i].clr(); 79 bool add = true; 80 tot = -1; 81 while (pos < size){ 82 if (s[pos] == '-') 83 ++ pos, add = false; 84 if (s[pos] == '+') 85 ++ pos, add = true; 86 87 if (s[pos] == 'n'){ 88 a[++tot].k = 1; 89 a[tot].a = 1; 90 ++ pos; 91 } 92 else{ 93 a[++tot].k = 0; 94 while (s[pos] >= '0' && s[pos] <= '9'){ 95 a[tot].k = a[tot].k * 10 + s[pos] - '0'; 96 ++ pos; 97 } 98 } 99 if (!add) a[tot].k *= -1; 100 101 if (s[pos] == 'n'){ 102 a[tot].a = 1; 103 ++ pos; 104 } 105 106 if (s[pos] == '^'){ 107 ++ pos; 108 a[tot].a = 0; 109 while (s[pos] >= '0' && s[pos] <= '9'){ 110 a[tot].a = a[tot].a * 10 + s[pos] - '0'; 111 ++ pos; 112 } 113 } 114 115 if (s[pos] == ')'){ 116 pos += 2; 117 break; 118 } 119 } 120 ++ tot; 121 122 mod = 0; 123 while (pos < size){ 124 mod = mod * 10 + s[pos] - '0'; 125 ++ pos; 126 } 127 } 128 129 bool is_int(int x) 130 { 131 int64 tmp = 0; 132 for (int i = 0; i < tot; ++ i){ 133 tmp = (tmp + (((int64)a[i].k * Mypow((int64)x, (int64)a[i].a)) % mod)) % mod; 134 } 135 if (!(tmp % mod)) return true; 136 return false; 137 } 138 139 bool solve () 140 { 141 for (int i = 1; i <= a[0].a+1; ++ i) 142 if (!is_int(i)) return false; 143 return true; 144 } 145 146 int main() 147 { 148 // freopen("a.in","r",stdin); 149 // freopen("a.out","w",stdout); 150 // std::ios::sync_with_stdio(false); 151 s.clear(); 152 int test = 0; 153 while (cin >> s && s[0] != '.'){ 154 init (); 155 156 cout << "Case " << ++test << ": "; 157 if (solve()) cout << "Always an integer" << endl; 158 else cout << "Not always an integer" << endl; 159 } 160 return 0; 161 }
3、UVa 11426
题意：输入正整数n，求gcd(1, 2) + gcd(1, 3) + gcd(2, 3) + gcd(1, 4) +...gcd(n-1, n)，多组数据。 (2 <= n <= 4 * 10^6, case <= 100)
解法：首先，问题可以转化为求f(m) = gcd(1, m) + gcd(2, m) + gcd(3, m) +...+ gcd(m-1, m)。
用g(i, m)表示1-m中与m最大公约数为i的数的个数，则f(m) = sum(i * g(i, m))，i表示m的所有约数。
由于对gcd(x, y) = a，有gcd(x/a, y/a) = 1，所以g(i, m) = phi(m/i)，phi()为欧拉函数。
综上，问题得到解决。但是，如果依次计算f(m)，会超时，原因是因为对每个m都需要枚举它的所有约数i，速度较慢，但如果对每个约数i枚举它的倍数m(并更新f(m)的值)，时间复杂度则降到了允许的范围内(最后我的程序在OJ上跑了2.472s)。
tag: math, number theory, 欧拉函数, 转化最大公约数为互质, good
Ps：本篇题解第1题(hdu 4675)为转化最大公约数为公约数，可以参考对比一下
1 /* 2 * Author: plum rain 3 * Created Time: 2013-09-02 07:15 4 * File Name: math-Uva-11426.cpp 5 */ 6 #include<iostream> 7 #include<sstream> 8 #include<fstream> 9 #include<vector> 10 #include<list> 11 #include<deque> 12 #include<queue> 13 #include<stack> 14 #include<map> 15 #include<set> 16 #include<bitset> 17 #include<algorithm> 18 #include<cstdio> 19 #include<cstdlib> 20 #include<cstring> 21 #include<cctype> 22 #include<cmath> 23 #include<ctime> 24 #include<utility> 25 26 using namespace std; 27 28 #define CLR(x) memset(x, 0, sizeof(x)) 29 #define PB push_back 30 #define SZ(v) ((int)(v).size()) 31 #define INF 999999999999 32 #define zero(x) (((x)>0?(x):-(x))<eps) 33 #define out(x) cout<<#x<<":"<<(x)<<endl 34 #define tst(a) cout<<#a<<endl 35 #define CINBEQUICKER std::ios::sync_with_stdio(false) 36 37 typedef vector<int> VI; 38 typedef vector<string> VS; 39 typedef vector<double> VD; 40 typedef long long int64; 41 42 const double eps = 1e-8; 43 const double PI = atan(1.0)*4; 44 const int maxint = 2139062143; 45 const int N = 4000000; 46 47 inline int Mymod( int a , int b ) { int x=a%b;if(x<0) x+=b;return x;} 48 49 int64 s[N+1], f[N+1]; 50 int64 phi[N+1]; 51 52 void phi_table(int n) 53 { 54 CLR(phi); 55 phi[1] = 1; 56 for (int i = 2; i <= n; ++ i) 57 if (!phi[i]){ 58 for (int j = i; j <= n; j += i){ 59 if (!phi[j]) phi[j] = j; 60 phi[j] = phi[j] / i * (i-1); 61 } 62 } 63 } 64 65 int main() 66 { 67 phi_table(N); 68 69 CLR (f); 70 for (int i = 1; i <= N; ++ i) 71 for (int j = i*2; j <= N; j += i) 72 f[j] += i * phi[j / i]; 73 74 s[2] = f[2]; 75 for (int i = 3; i <= N; ++ i) 76 s[i] = s[i-1] + f[i]; 77 78 int n; 79 while (scanf ("%d", &n) == 1 && n) 80 printf ("%lld\n", s[n]); 81 82 return 0; 83 }
4、UVa 11754 - Code Feat
题意：有一个正整数N满足C个条件，每个条件都形如“它除以X的余数在集合{Y1,Y2...Yk}中”，所有条件中的X两两互素，求最小的S个解。(1 <= C <= 9, 1 <= S <= 10, 1 <= k <= 100, 0 <= Yi < X)
解法：首先，很容易想到两种解法。
一、中国剩余定理+搜索。对每个条件，枚举除X的余数为Yi(1 <= i <= sizeof{Yn})，然后用中国剩余定理解模方程，取最小的S个解。这种方法的时间复杂度为O(∏k)(∏k即为所有条件中集合Y的大小之积)
二、暴力枚举。选择所有条件中k / X最小的条件，先对该条件中的Yi排序，然后按照t = 0, 1, 2,...的顺序枚举所有tX + Yi(相同的t按照从小到大的顺序枚举Yi)，看看是否满足条件。具体时间复杂度不好估计，但由于∏k很大，所以很快就能找到解。(因为每次循环会找出所有t * x - (t+1) * x中的解，而k < x，这也是选择k / x小的条件的原因)
那么不难发现，上面两种思路刚好互补，结合一下就是正解了。
tag:maht, number theory, 中国剩余定理, good
1 /* 2 * Author: Plumrain 3 * Created Time: 2013-09-04 13:10 4 * File Name: math-UVa-117542.cpp 5 */ 6 #include<iostream> 7 #include<cstdio> 8 #include<algorithm> 9 #include<vector> 10 #include<set> 11 12 using namespace std; 13 14 #define PB push_back 15 #define SZ(v) ((int)(v).size()) 16 17 typedef long long int64; 18 19 const int Limit = 10000; 20 21 int n, m; 22 int64 a[10][105], x[10]; 23 int64 pro; 24 set<int64> b[10]; 25 vector<int64> ans; 26 27 int init(int64& tot) 28 { 29 int ret = 0; 30 pro = 1; tot = 1; 31 for (int i = 0; i < n; ++ i){ 32 scanf ("%lld%lld", &x[i], &a[i][0]); 33 pro *= x[i]; tot *= a[i][0]; 34 for (int j = 1; j <= a[i][0]; ++ j) 35 scanf ("%lld", &a[i][j]); 36 sort(a[i]+1, a[i]+a[i][0]+1); 37 38 if ((double)(a[i][0] / x[i]) < (double)(a[ret][0] / x[ret])) 39 ret = i; 40 } 41 return ret; 42 } 43 44 void gao1(int flag) 45 { 46 for (int i = 0; i < n; ++ i) if (i != flag){ 47 b[i].clear(); 48 for (int j = 1; j <= a[i][0]; ++ j) 49 b[i].insert(a[i][j]); 50 } 51 52 for (int t = 0; m != 0; ++ t) 53 for (int i = 1; i <= a[flag][0]; ++ i){ 54 int64 tmp = t * x[flag] + a[flag][i]; 55 56 bool ok = true; 57 for (int j = 0; j < n; ++ j) if (j != flag){ 58 if (!b[j].count(tmp % x[j])){ 59 ok = false; break; 60 } 61 } 62 63 if (ok){ 64 if (!tmp) continue; 65 printf ("%lld\n", tmp); 66 if (--m == 0) break; 67 } 68 } 69 } 70 71 void DFS(int d, int64 v) 72 { 73 if (d == n){ 74 v %= pro; 75 if (v < 0) v += pro; 76 ans.PB(v); 77 return ; 78 } 79 80 for (int i = 1; i <= a[d][0]; ++ i) 81 DFS(d+1, (v + a[d][i]) % pro); 82 } 83 84 void ext_gcd(int64 a, int64 b, int64& d, int64 &x, int64& y) 85 { 86 if (!b) d = a, x = 1, y = 0; 87 else{ 88 ext_gcd(b, a%b, d, y, x); 89 y -= x * (a / b); 90 } 91 } 92 93 void gao2() 94 { 95 for (int i = 0; i < n; ++ i){ 96 int64 w = pro / x[i]; 97 int64 d, xx, y; 98 ext_gcd (x[i], w, d, xx, y); 99 int64 tmp = (w * y) % pro; 100 for (int j = 1; j <= a[i][0]; ++ j) 101 a[i][j] = (a[i][j] * tmp) % pro; 102 } 103 104 ans.clear(); 105 DFS(0, 0); 106 107 sort(ans.begin(), ans.end()); 108 if (ans[0] != 0) 109 printf ("%lld\n", ans[0]), -- m; 110 int64 tmp = ans[0], add = 0; 111 int fla = 1, size = SZ (ans); 112 while (m > 0){ 113 if (fla == size) 114 fla = 0, add += pro; 115 if (ans[fla] + add == tmp || ans[fla] + add == 0){ 116 ++ fla; continue; 117 } 118 119 tmp = ans[fla] + add; 120 printf ("%lld\n", tmp); 121 ++ fla; -- m; 122 } 123 } 124 125 int main() 126 { 127 while (scanf ("%d%d", &n, &m) != EOF && n && m){ 128 int64 tot; 129 int flag = init (tot); 130 131 if (tot > Limit) gao1 (flag); 132 else gao2 (); 133 printf ("\n"); 134 } 135 return 0; 136 }
5、UVa 11916 - Emoogle Grid
题意：有这样一道题目，给一个M行N列的网格涂K种颜色(不一定所有颜色都涂)，其中B个格子不用上色，其他格子每格涂一种颜色，且同列相邻的格子不能同色，给出B个不上色格子的位置，求涂色方案数模100000007的结果R。现在，已知N，K，R，B和B个格子的位置(xi, yi)，求M。保证输入有解。
1 <= M, N <= 10^8，1 <= xi <= M，1 <= yi <= N，0 <= B <= 500，0 <= R < mod。
解法：其实，思路很简单，首先找到max{xi}，然后分类讨论，M可能为max{xi}，max{xi} + 1, 还有更大。前两种情况可以求出上色方案数取模看是否等于R，最后一种情况需要用到Shank的大步小步算法来解高阶同余方程。
tag: math, number theory, Shank's,
Ps：这道题不难，但是把很多数论部分的算法都考查到了，所以很不错！
PPs：这道题真是让我认识到了一点....刘汝佳的模板虽然好用，通用性强，效率高，但是有很多他个人的习惯。我用他的模板时改过一些地方，虽然是正确的改动，但是和他的习惯冲突了，所以当我多个刘汝佳模板一起用时就导致了很隐蔽的错误- -
Source Code就不贴了...下次写个比较好的代码再贴上来- -
6、POJ 3101 Astronomy
题意：有n个行星围绕着一个中心转动，称所有行星全部在一条直线上为“planet parade”。给出每个行星旋转一圈所需时间t[i]年，问两次planet parade间隔的时间为多长。分数形式输出。
2 <= n <= 1000，1 <= t[i] <= 10000。
解法：设半圈的长度为1，对于任意两颗行星i和j，他们每年旋转的长度差距为2 * abs(1/t[i] - 1/t[j])。则由n颗行星算出n-1个长度差x[i]/y[i]，题目即是要找到一个最小的分数T，使得T与任意长度差相乘之积都为整数。 所以T = lcm(y[i]) / gcd(x[i])。数据太大要用到高精度。
tag:math, 高精度, think
Ps：思路懂了，要用到高精度所以我就懒得写了.......
7、POJ 2429 GCD & LCM Inverse
题意：对于两个正整数a, b，给出gcd(a, b)和lcm(a, b)，求a和b。多组数据的话，输出a + b最小的一组。gcd(a, b), lcm(a, b) < 2^63。
解法：记mul = lcm(a, b)/ gcd(a, b)，题目相当于求两个互质的数a'和b'，使得a' * b' = mul，并且a' + b'最小。将mul分解质因数，然后枚举就好了，这样的时间复杂度是能接受的，250ms过。分解质因数的时候要用到miller_rabin和Pollard，因为数据太大。
tag:math, number theory, miller_rabin, Pollard
1 /* 2 * Author: Plumrain 3 * Created Time: 2013-10-24 09:16 4 * File Name: math-POJ-1811.cpp 5 */ 6 #include<iostream> 7 #include<cstdio> 8 #include<algorithm> 9 #include<ctime> 10 11 using namespace std; 12 13 const int MT = 5; 14 const int TIM = 240; 15 typedef long long int64; 16 const int64 SPE = 46856248255981; 17 18 int cnt; 19 int64 prm[5] = {2, 3, 7, 61, 24251}; 20 int64 fac[10000], all, an[10000], bn[10000]; 21 22 bool cmp(int64 a, int64 b) 23 { 24 return a < b; 25 } 26 27 int64 Mypow(int64 p, int64 n) 28 { 29 int64 ret = 1; 30 while (n){ 31 if (n & 1) ret *= p; 32 n >>= 1; 33 p *= p; 34 } 35 return ret; 36 } 37 38 int64 min(int64 a, int64 b) 39 { 40 return a < b ? a : b; 41 } 42 43 int64 gcd(int64 a, int64 b) 44 { 45 return b ? gcd(b, a%b) : a; 46 } 47 48 int64 mul_mod(int64 p, int64 q, int64 mod) 49 { 50 int64 ret = 0; 51 p %= mod; 52 while (q){ 53 if (q & 1){ 54 ret += p; 55 if (ret >= mod) ret -= mod; 56 } 57 p <<= 1; q >>= 1; 58 if (p >= mod) p -= mod; 59 } 60 return ret; 61 } 62 63 int64 pow_mod(int64 p, int64 n, int64 mod) 64 { 65 int64 ret = 1; 66 p %= mod; 67 while (n > 0){ 68 if (n & 1) ret = mul_mod(ret, p, mod); 69 n >>= 1; 70 p = mul_mod(p, p, mod); 71 } 72 return ret; 73 } 74 75 bool Wintess(int64 aa, int64 n, int64 mod) 76 { 77 int t = 0; 78 while ((n & 1) == 0){ 79 n >>= 1; 80 ++ t; 81 } 82 83 int64 x = pow_mod(aa, n, mod), y; 84 while (t --){ 85 y = mul_mod (x, x, mod); 86 if (y == 1 && x != 1 && x != mod-1) return 1; 87 x = y; 88 } 89 return (y != 1); 90 } 91 92 bool miller_rabin(int64 n, int tim) 93 { 94 if (n == 2) return true; 95 if (n == 1 || (n & 1) == 0 || n == SPE) return false; 96 97 for (int i = 0; i < tim; ++ i){ 98 int64 aa = prm[i]; 99 if (aa == n) return true; 100 if (Wintess(aa, n-1, n)) return false; 101 } 102 return true; 103 } 104 105 int64 pollard(int64 n, int c) 106 { 107 srand(time(NULL)); 108 int64 i = 1, k = 2, x = rand()%n; 109 int64 y = x; 110 while (1){ 111 ++ i; 112 x = (mul_mod(x, x, n) + c) % n; 113 int64 d = gcd(y-x, n); 114 if (d > 1 && d < n) return d; 115 if (y == x) return n; 116 if (i == k){ 117 y = x; 118 k <<= 1; 119 } 120 } 121 } 122 123 void get_prime(int64 n, int c) 124 { 125 if (n == 1) return; 126 if (miller_rabin(n, MT)){ 127 fac[cnt++] = n; 128 return; 129 } 130 131 int64 m = n; 132 while (m == n) m = pollard(m, c--); 133 get_prime(m, c); 134 get_prime(n/m, c); 135 } 136 137 void gao(int64 mul, int64 g) 138 { 139 int64 tmp = fac[0]; all = 0; 140 an[0] = fac[0]; bn[0] = 0; 141 for (int i = 0; i < cnt; ++ i){ 142 if (tmp != fac[i]){ 143 tmp = fac[i]; 144 an[++all] = tmp; 145 bn[all] = 1; 146 } 147 else 148 ++ bn[all]; 149 } 150 ++ all; 151 152 int64 ans = 1, sum = mul + 1; 153 for (int64 i = 0; i < (int64)1<<all; ++ i){ 154 int64 tmp = 1; 155 for (int64 j = 0; j < all; ++ j) 156 if (i & (int64)1<<j) 157 tmp *= Mypow(an[j], bn[j]); 158 159 if (sum > tmp + mul/tmp){ 160 sum = tmp + mul/tmp; 161 ans = tmp; 162 } 163 } 164 165 if (ans > mul/ans) ans = mul / ans; 166 printf ("%lld %lld\n", ans*g, mul*g/ans); 167 } 168 169 int main() 170 { 171 int64 g, l; 172 while (scanf ("%lld%lld", &g, &l) != EOF){ 173 int64 mul = l / g; 174 cnt = 0; 175 get_prime(mul, TIM); 176 sort(fac, fac+cnt, cmp); 177 gao(mul, g); 178 } 179 return 0; 180 }
8、POJ 2689 Prime Distance
题意：任意给定两个数l和r，求区间[l, r]上，相邻素数之差最小和最大的分别是哪两组素数。
解法：素数二次筛法。即首先将sqrt(r)内的所有素数全部筛出来，然后，用线性的方法筛掉[l, r]内的所有合数。注意一点，若r太大，而r - l不太大，在用bool v[]记录素数合数的时候可以将区间[l, r]平移到[0, r-l]。
tag:math, number theory, 素数二次筛法
1 /* 2 * Author: Plumrain 3 * Created Time: 2013-10-24 16:35 4 * File Name: math-POJ-2689.cpp 5 */ 6 #include<iostream> 7 #include<cstdio> 8 #include<cstring> 9 10 using namespace std; 11 12 #define CLR(x) memset(x, 0, sizeof(x)) 13 const int maxint = 2147483647; 14 const int N = 50000; 15 typedef long long int64; 16 17 bool vis[N+5], v[1000005]; 18 int64 prm[N+5]; 19 int all; 20 21 int64 max(int64 a, int64 b) 22 { 23 return a < b ? b : a; 24 } 25 26 void sieve(int n) 27 { 28 CLR (vis); 29 for (int64 i = 2; i <= n; ++ i) if (!vis[i]) 30 for (int64 j = i*i; j <= n; j += i) vis[j] = 1; 31 } 32 33 int primes(int n) 34 { 35 sieve(n); 36 int ret = 0; 37 for (int i = 2; i <= n; ++ i) 38 if (!vis[i]) prm[ret++] = i; 39 return ret; 40 } 41 42 int main() 43 { 44 all = primes(N); 45 int64 m, n; 46 while (scanf ("%lld%lld", &m ,&n) != EOF){ 47 m = max((int64)2, m); 48 CLR (v); 49 for (int i = 0; i < all && prm[i]*prm[i] <= n; ++ i){ 50 int64 t = m / prm[i] * prm[i]; 51 if (t != m) t += prm[i]; 52 while (t <= n){ 53 v[t - m] = 1; 54 t += prm[i]; 55 } 56 if (prm[i] >= m) v[prm[i]-m] = 0; 57 } 58 59 int64 a1 = 0, a2 = 0, b1 = 0, b2 = maxint, tmp = -1; 60 for (int64 i = 0; i <= n-m; ++ i) if (!v[i]){ 61 if (tmp == -1){tmp = i; continue;} 62 if (a2 - a1 < i - tmp){ 63 a1 = tmp; 64 a2 = i; 65 } 66 if (b2 - b1 > i - tmp){ 67 b2 = i; 68 b1 = tmp; 69 } 70 tmp = i; 71 } 72 if (!a1 && !a2) printf ("There are no adjacent primes.\n"); 73 else 74 printf ("%lld,%lld are closest, %lld,%lld are most distant.\n", b1+m, b2+m, a1+m, a2+m); 75 } 76 return 0; 77 }
9、ZOJ 2562 More Divisors
题意：每给一个n<=10^16，求n以内最大的反素数。反素数即是：对正整数x记其约数的个数记做g(x).例如g(1)=1,g(6)=4.如果某个正整数x满足:对于任意i(0<i<x),都有g(i)<g(x),则称x为反素数.
解法：反素数有一个性质，即是对任一个反素数p，p = (2^t0) * (3^t1) * (5^t2)...，且t0 >= t1 >= t2...这样的的话，就可以直接利用搜索了。
tag:math, 反素数，DFS
1 /* 2 * Author: Plumrain 3 * Created Time: 2013-10-25 00:18 4 * File Name: math-ZOJ-2562.cpp 5 */ 6 #include<iostream> 7 #include<cstdio> 8 #include<vector> 9 10 using namespace std; 11 12 #define PB push_back 13 typedef long long int64; 14 15 vector<int64> cur; 16 int64 cnt, ans, N; 17 int64 prm[15] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}; 18 19 int64 min(int64 a, int64 b) 20 { 21 return a < b ? a : b; 22 } 23 24 void DFS(int64 n, int64 pos, int64 num) 25 { 26 int64 tmp = 1; 27 for (int64 i = 0; i < (int64)cur.size(); ++ i) 28 tmp *= (cur[i] + 1); 29 if (tmp > cnt){ 30 ans = n; 31 cnt = tmp; 32 } 33 else if (tmp == cnt) 34 ans = min(ans, n); 35 36 if (pos >= 14) return; 37 38 int64 ttmp = prm[++pos] * n; 39 int64 ntmp = 1; 40 while (ttmp <= N && ntmp <= num){ 41 cur.PB(ntmp); 42 DFS (ttmp, pos, ntmp); 43 cur.pop_back(); 44 45 ttmp *= prm[pos]; 46 ++ ntmp; 47 } 48 } 49 50 int main() 51 { 52 while (scanf ("%lld", &N) != EOF){ 53 cur.clear(); 54 int64 tmp = 2, num = 1; 55 cnt = ans = 1; 56 while (tmp <= N){ 57 cur.PB(num); 58 DFS (tmp, 0, num); 59 cur.pop_back(); 60 tmp *= 2; 61 ++ num; 62 } 63 64 printf ("%lld\n", ans); 65 } 66 return 0; 67 }
扫一扫
1238
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
