数论专题训练

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 ("
");
124     }
125     return 0;
126 }

View Code

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 }

View Code

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
", s[n]);
81 
82     return 0;
83 }

View Code

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
", 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
", 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
", 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 ("
");
134     }
135     return 0;
136 }

View Code

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
", 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 }

View Code

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.
");
73         else
74             printf ("%lld,%lld are closest, %lld,%lld are most distant.
", b1+m, b2+m, a1+m, a2+m);
75     }
76     return 0;
77 }

View Code

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
", ans);
65     }
66     return 0;
67 }

View Code

------------------------------------------------------------------
现在的你，在干什么呢？
你是不是还记得，你说你想成为岩哥那样的人。