数论专题训练

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

/*
 * Author:  plum rain
 * Created Time:  2013-08-16 11:10
 * File Name: fuck.cpp
 */
#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<utility>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))
#define PB push_back
#define SZ(v) ((int)(v).size())
#define INF 999999999999
#define zero(x) (((x)>0?(x):-(x))<eps)
#define out(x) cout<<#x<<":"<<(x)<<endl
#define tst(a) cout<<#a<<endl
#define CINBEQUICKER std::ios::sync_with_stdio(false)

typedef vector<int> VI;
typedef vector<string> VS;
typedef vector<double> VD;
typedef long long int64;

const double eps = 1e-8;
const double PI = atan(1.0)*4;
const int maxint = 2139062143;
const int N = 300005;
const int mod = 1000000007;

inline int Mymod( int a , int b ) { int x=a%b;if(x<0) x+=b;return x;}

int n, m, k;
int an[N];
int64 f[N];
int64 C[N];

int64 Mypow(int64 a, int64 n)
{
    int64 ret = 1, tmp = a % mod;
    while (n){
        if (n&1)
            ret = (ret * tmp) % mod;
        tmp = (tmp * tmp) % mod;
        n >>= 1;
    }
    return ret;
}

void Cinit()
{
    CLR (C);
    int64 x = n - k;
    C[x] = 1;
    for (int i = x; i < n; ++ i){
        C[i+1] = (C[i] * (i + 1)) % mod;
        C[i+1] = (C[i+1] * Mypow(i+1-x, mod-2)) % mod;
    }
}


void solve(int x)
{
    int cnt = 0;
    int64 sum = 0;
    for (int i = 1; i*x <= m; ++ i){
        cnt += an[i*x];
        if (i > 1) sum = (sum + f[i*x]) % mod;
    }
    int t = m / x;
    if (t == 1){
        if (cnt == n - k) f[x] = 1;
        else f[x] = 0;
        return ;
    }
    if (cnt < n - k){
        f[x] = 0;
        return;
    }
    int64 tmp = (C[cnt] * Mypow(t, n-cnt)) % mod; 
    tmp = (tmp * Mypow(t-1, cnt - n + k)) % mod;
    f[x] = (tmp - sum + mod) % mod;
}

int main()
{
//    freopen("tst.in","r",stdin);
//    freopen("fuck.out","w",stdout);
//    std::ios::sync_with_stdio(false);
    while (scanf ("%d%d%d", &n, &m, &k) != EOF){ 
        CLR (an);
        int x;
        for (int i = 0; i < n; ++ i){
            scanf ("%d", &x);
            ++ an[x];
        }

        CLR (f);
        Cinit();
        for (int i = m; i >= 1; -- i)
            solve(i);

        printf ("%I64d", f[1]);
        for (int i = 2; i <= m; ++ i)
            printf (" %I64d", f[i]);
        printf ("
");
    }
    return 0;
}

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

/*
 * Author:  plum rain
 * Created Time:  2013-08-23 21:00
 * File Name: (good)-math-WF2008-LA-4119.cpp
 */
#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<utility>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))
#define PB push_back
#define SZ(v) ((int)(v).size())
#define INF 999999999999
#define zero(x) (((x)>0?(x):-(x))<eps)
#define out(x) cout<<#x<<":"<<(x)<<endl
#define tst(a) cout<<#a<<endl
#define CINBEQUICKER std::ios::sync_with_stdio(false)

typedef vector<int> VI;
typedef vector<string> VS;
typedef vector<double> VD;
typedef long long int64;

const double eps = 1e-8;
const double PI = atan(1.0)*4;
const int maxint = 2139062143;

inline int Mymod( int a , int b ) { int x=a%b;if(x<0) x+=b;return x;}

struct P{
    int k, a, b;
    //k * n ^ a
    void clr(){
        k = 0; a = 0;
    }
};

P a[105];
int tot;
string s;
int mod;

int64 Mypow(int64 p, int64 n)
{
    int64 ret = 1;
    while (n > 0){
        if (n & 1)
            ret = (ret * p) % mod;
        n >>= 1;
        p = (p * p) % mod; 
    }
    return ret;
}

void init()
{
    int size = SZ (s);
    int pos = 1;
    for (int i = 0; i < 105; ++ i)
        a[i].clr();
    bool add = true;
    tot = -1;
    while (pos < size){
        if (s[pos] == '-')
            ++ pos, add = false;
        if (s[pos] == '+')
            ++ pos, add = true;

        if (s[pos] == 'n'){
            a[++tot].k = 1;
            a[tot].a = 1;
            ++ pos;
        }
        else{ 
            a[++tot].k = 0;
            while (s[pos] >= '0' && s[pos] <= '9'){
                a[tot].k = a[tot].k * 10 + s[pos] - '0';
                ++ pos;
            }
        }
        if (!add) a[tot].k *= -1;

        if (s[pos] == 'n'){
            a[tot].a = 1;
            ++ pos;
        }

        if (s[pos] == '^'){
            ++ pos;
            a[tot].a = 0;
            while (s[pos] >= '0' && s[pos] <= '9'){
                a[tot].a = a[tot].a * 10 + s[pos] - '0';
                ++ pos;
            }
        }

        if (s[pos] == ')'){
            pos += 2;
            break;
        }
    }
    ++ tot;

    mod = 0;
    while (pos < size){ 
        mod = mod * 10 + s[pos] - '0';
        ++ pos;
    }
}

bool is_int(int x)
{
    int64 tmp = 0;
    for (int i = 0; i < tot; ++ i){
        tmp = (tmp + (((int64)a[i].k * Mypow((int64)x, (int64)a[i].a)) % mod)) % mod;
    }
    if (!(tmp % mod)) return true;
    return false;
}

bool solve ()
{
    for (int i = 1; i <= a[0].a+1; ++ i)
        if (!is_int(i)) return false;
    return true;
}

int main()
{
//    freopen("a.in","r",stdin);
//    freopen("a.out","w",stdout);
//    std::ios::sync_with_stdio(false);
    s.clear();
    int test = 0;
    while (cin >> s && s[0] != '.'){
        init ();

        cout << "Case " << ++test << ": ";
        if (solve()) cout << "Always an integer" << endl;
        else cout << "Not always an integer" << endl;
    }
    return 0;
}

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)为转化最大公约数为公约数，可以参考对比一下

/*
 * Author:  plum rain
 * Created Time:  2013-09-02 07:15
 * File Name: math-Uva-11426.cpp
 */
#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<utility>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))
#define PB push_back
#define SZ(v) ((int)(v).size())
#define INF 999999999999
#define zero(x) (((x)>0?(x):-(x))<eps)
#define out(x) cout<<#x<<":"<<(x)<<endl
#define tst(a) cout<<#a<<endl
#define CINBEQUICKER std::ios::sync_with_stdio(false)

typedef vector<int> VI;
typedef vector<string> VS;
typedef vector<double> VD;
typedef long long int64;

const double eps = 1e-8;
const double PI = atan(1.0)*4;
const int maxint = 2139062143;
const int N = 4000000;

inline int Mymod( int a , int b ) { int x=a%b;if(x<0) x+=b;return x;}

int64 s[N+1], f[N+1];
int64 phi[N+1];

void phi_table(int n)
{
    CLR(phi);
    phi[1] = 1;
    for (int i = 2; i <= n; ++ i)
        if (!phi[i]){
            for (int j = i; j <= n; j += i){
                if (!phi[j]) phi[j] = j;
                phi[j] = phi[j] / i * (i-1);
            }
        }
}

int main()
{
    phi_table(N);

    CLR (f);
    for (int i = 1; i <= N; ++ i)
        for (int j = i*2; j <= N; j += i)
            f[j] += i * phi[j / i];

    s[2] = f[2];
    for (int i = 3; i <= N; ++ i)
        s[i] = s[i-1] + f[i];
    
    int n;
    while (scanf ("%d", &n) == 1 && n)
        printf ("%lld
", s[n]);

    return 0;
}

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

/*
 * Author:  Plumrain
 * Created Time:  2013-09-04 13:10
 * File Name: math-UVa-117542.cpp
 */
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<vector>
#include<set>

using namespace std;

#define PB push_back
#define SZ(v) ((int)(v).size())

typedef long long int64;

const int Limit = 10000;

int n, m;
int64 a[10][105], x[10];
int64 pro;
set<int64> b[10];
vector<int64> ans;

int init(int64& tot)
{
    int ret = 0;
    pro = 1; tot = 1;
    for (int i = 0; i < n; ++ i){
        scanf ("%lld%lld", &x[i], &a[i][0]);
        pro *= x[i]; tot *= a[i][0];
        for (int j = 1; j <= a[i][0]; ++ j)
            scanf ("%lld", &a[i][j]);
        sort(a[i]+1, a[i]+a[i][0]+1);

        if ((double)(a[i][0] / x[i]) < (double)(a[ret][0] / x[ret]))
            ret = i;
    }
    return ret;
}

void gao1(int flag)
{
    for (int i = 0; i < n; ++ i) if (i != flag){
        b[i].clear();
        for (int j = 1; j <= a[i][0]; ++ j)
            b[i].insert(a[i][j]);
    }

    for (int t = 0; m != 0; ++ t)
        for (int i = 1; i <= a[flag][0]; ++ i){
            int64 tmp = t * x[flag] + a[flag][i];

            bool ok = true;
            for (int j = 0; j < n; ++ j) if (j != flag){
                if (!b[j].count(tmp % x[j])){
                    ok = false; break;
                }
            }

            if (ok){
                if (!tmp) continue;
                printf ("%lld
", tmp);
                if (--m == 0) break;
            }
        }
}

void DFS(int d, int64 v)
{
    if (d == n){ 
        v %= pro;
        if (v < 0) v += pro;
        ans.PB(v);
        return ;
    }

    for (int i = 1; i <= a[d][0]; ++ i)
        DFS(d+1, (v + a[d][i]) % pro);
}

void ext_gcd(int64 a, int64 b, int64& d, int64 &x, int64& y)
{
    if (!b) d = a, x = 1, y = 0;
    else{
        ext_gcd(b, a%b, d, y, x);
        y -= x * (a / b);
    }
}

void gao2()
{
    for (int i = 0; i < n; ++ i){ 
        int64 w = pro / x[i];
        int64 d, xx, y;
        ext_gcd (x[i], w, d, xx, y);
        int64 tmp = (w * y) % pro;
        for (int j = 1; j <= a[i][0]; ++ j)
            a[i][j] = (a[i][j] * tmp) % pro;
    }

    ans.clear();
    DFS(0, 0);

    sort(ans.begin(), ans.end());
    if (ans[0] != 0)
        printf ("%lld
", ans[0]), -- m;
    int64 tmp = ans[0], add = 0;
    int fla = 1, size = SZ (ans);
    while (m > 0){
        if (fla == size)
            fla = 0, add += pro;
        if (ans[fla] + add == tmp || ans[fla] + add == 0){
            ++ fla; continue;
        }

        tmp = ans[fla] + add;
        printf ("%lld
", tmp);
        ++ fla; -- m;
    }
}

int main()
{
    while (scanf ("%d%d", &n, &m) != EOF && n && m){    
        int64 tot;
        int flag = init (tot);

        if (tot > Limit) gao1 (flag); 
        else gao2 ();
        printf ("
");
    }
    return 0;
}

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

/*
 * Author:  Plumrain
 * Created Time:  2013-10-24 09:16
 * File Name: math-POJ-1811.cpp
 */
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<ctime>

using namespace std;

const int MT = 5;
const int TIM = 240;
typedef long long int64;
const int64 SPE = 46856248255981;

int cnt;
int64 prm[5] = {2, 3, 7, 61, 24251};
int64 fac[10000], all, an[10000], bn[10000];

bool cmp(int64 a, int64 b)
{
    return a < b;
}

int64 Mypow(int64 p, int64 n)
{
    int64 ret = 1;
    while (n){
        if (n & 1) ret *= p;
        n >>= 1;
        p *= p;
    }
    return ret;
}

int64 min(int64 a, int64 b)
{
    return a < b ? a : b;
}

int64 gcd(int64 a, int64 b)
{
    return b ? gcd(b, a%b) : a;
}

int64 mul_mod(int64 p, int64 q, int64 mod)
{
     int64 ret = 0;
     p %= mod;
     while (q){
         if (q & 1){
             ret += p;
             if (ret >= mod) ret -= mod;
         }
         p <<= 1; q >>= 1;
         if (p >= mod) p -= mod;
     }
     return ret;
}

int64 pow_mod(int64 p, int64 n, int64 mod)
{
    int64 ret = 1;
    p %= mod;
    while (n > 0){
        if (n & 1) ret = mul_mod(ret, p, mod);
        n >>= 1;
        p = mul_mod(p, p, mod);
    }
    return ret;
}

bool Wintess(int64 aa, int64 n, int64 mod)
{
    int t = 0;
    while ((n & 1) == 0){
        n >>= 1;
        ++ t;
    }

    int64 x = pow_mod(aa, n, mod), y;
    while (t --){
        y = mul_mod (x, x, mod);
        if (y == 1 && x != 1 && x != mod-1) return 1;
        x = y;
    }
    return (y != 1);
}

bool miller_rabin(int64 n, int tim)
{
    if (n == 2) return true;
    if (n == 1 || (n & 1) == 0 || n == SPE) return false;

    for (int i = 0; i < tim; ++ i){
        int64 aa = prm[i];
        if (aa == n) return true;
        if (Wintess(aa, n-1, n)) return false;
    }
    return true;
}

int64 pollard(int64 n, int c)
{
    srand(time(NULL));
    int64 i = 1, k = 2, x = rand()%n;
    int64 y = x;
    while (1){
        ++ i;
        x = (mul_mod(x, x, n) + c) % n;
        int64 d = gcd(y-x, n);
        if (d > 1 && d < n) return d;
        if (y == x) return n;
        if (i == k){
            y = x;
            k <<= 1;
        }
    }
}

void get_prime(int64 n, int c)
{
    if (n == 1) return;
    if (miller_rabin(n, MT)){
        fac[cnt++] = n;
        return;
    }

    int64 m = n;
    while (m == n) m = pollard(m, c--);
    get_prime(m, c);
    get_prime(n/m, c);
}

void gao(int64 mul, int64 g)
{ 
    int64 tmp = fac[0]; all = 0;
    an[0] = fac[0]; bn[0] = 0;
    for (int i = 0; i < cnt; ++ i){
        if (tmp != fac[i]){
            tmp = fac[i];
            an[++all] = tmp;
            bn[all] = 1;
        }
        else
            ++ bn[all];
    }
    ++ all;

    int64 ans = 1, sum = mul + 1;
    for (int64 i = 0; i < (int64)1<<all; ++ i){
        int64 tmp = 1;
        for (int64 j = 0; j < all; ++ j)
            if (i & (int64)1<<j)
                tmp *= Mypow(an[j], bn[j]);

        if (sum > tmp + mul/tmp){
            sum = tmp + mul/tmp;
            ans = tmp;
        }
    }

    if (ans > mul/ans) ans = mul / ans;
    printf ("%lld %lld
", ans*g, mul*g/ans);
}

int main()
{
    int64 g, l;
    while (scanf ("%lld%lld", &g, &l) != EOF){
        int64 mul = l / g;
        cnt = 0;
        get_prime(mul, TIM);            
        sort(fac, fac+cnt, cmp);
        gao(mul, g);
    }
    return 0;
}

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, 素数二次筛法

/*
 * Author:  Plumrain
 * Created Time:  2013-10-24 16:35
 * File Name: math-POJ-2689.cpp
 */
#include<iostream>
#include<cstdio>
#include<cstring>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))
const int maxint = 2147483647;
const int N = 50000;
typedef long long int64;

bool vis[N+5], v[1000005];
int64 prm[N+5];
int all;

int64 max(int64 a, int64 b)
{
    return a < b ? b : a;
}

void sieve(int n)
{
    CLR (vis);
    for (int64 i = 2; i <= n; ++ i) if (!vis[i])
        for (int64 j = i*i; j <= n; j += i) vis[j] = 1;
}

int primes(int n)
{
    sieve(n);
    int ret = 0;
    for (int i = 2; i <= n; ++ i)
        if (!vis[i]) prm[ret++] = i;
    return ret;
}

int main()
{
    all = primes(N);
    int64 m, n;
    while (scanf ("%lld%lld", &m ,&n) != EOF){
        m = max((int64)2, m);
        CLR (v);
        for (int i = 0; i < all && prm[i]*prm[i] <= n; ++ i){
            int64 t = m / prm[i] * prm[i];
            if (t != m) t += prm[i];
            while (t <= n){
                v[t - m] = 1;
                t += prm[i];
            }
            if (prm[i] >= m) v[prm[i]-m] = 0;
        }

        int64 a1 = 0, a2 = 0, b1 = 0, b2 = maxint, tmp = -1;
        for (int64 i = 0; i <= n-m; ++ i) if (!v[i]){
            if (tmp == -1){tmp = i; continue;}
            if (a2 - a1 < i - tmp){
                a1 = tmp;
                a2 = i;
            }
            if (b2 - b1 > i - tmp){
                b2 = i;
                b1 = tmp;
            }
            tmp = i;
        }
        if (!a1 && !a2) printf ("There are no adjacent primes.
");
        else
            printf ("%lld,%lld are closest, %lld,%lld are most distant.
", b1+m, b2+m, a1+m, a2+m);
    }
    return 0;
}

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

/*
 * Author:  Plumrain
 * Created Time:  2013-10-25 00:18
 * File Name: math-ZOJ-2562.cpp
 */
#include<iostream>
#include<cstdio>
#include<vector>

using namespace std;

#define PB push_back
typedef long long int64;

vector<int64> cur;
int64 cnt, ans, N;
int64 prm[15] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47};

int64 min(int64 a, int64 b)
{
    return a < b ? a : b;
}

void DFS(int64 n, int64 pos, int64 num)
{
    int64 tmp = 1;
    for (int64 i = 0; i < (int64)cur.size(); ++ i)
        tmp *= (cur[i] + 1);
    if (tmp > cnt){
        ans = n; 
        cnt = tmp;
    }
    else if (tmp == cnt) 
        ans = min(ans, n);

    if (pos >= 14) return; 

    int64 ttmp = prm[++pos] * n;
    int64 ntmp = 1;
    while (ttmp <= N && ntmp <= num){
        cur.PB(ntmp);
        DFS (ttmp, pos, ntmp);
        cur.pop_back();

        ttmp *= prm[pos];
        ++ ntmp;
    }
}

int main()
{
    while (scanf ("%lld", &N) != EOF){
        cur.clear();
        int64 tmp = 2, num = 1;
        cnt = ans = 1;
        while (tmp <= N){
            cur.PB(num);
            DFS (tmp, 0, num);
            cur.pop_back();
            tmp *= 2;
            ++ num;
        }

        printf ("%lld
", ans);
    }
    return 0;
}