zoj 3547 The Boss on Mars 第36届ACM大连预选赛I题

以下转载于：

http://blog.csdn.net/xieshimao/article/details/6840731

数论题！

求与N不互质的数的K次方（K=4），反过来想若知道与N互质的K次方和，那所求就容易多了哦。

观察到与n互质的数的性质

比如12=2*2*3

那么与12不互质的数就有2,3,4,6,8,9,10,12

其实就是2的所有倍数，以及3的所有倍数

所以可以先求一个1到12的所有数的四次方和。这个有公式：

n*(n+1)*(2*n+1)*(3*n*n+3*n-1)/30

注意对与除以30可以看成是乘以30的逆元（对于1000000007的逆元）。

求的所有的四次方和之后当然要减去那些不互质的数的四次方

也就是说分别剪去了2和3的倍数的四次方，注意这里2和3的公倍数被多减去了，所以要加回来

所以容斥原理也要用

/*
 * zoj3547.c
 *
 *  Created on: 2011-10-3
 *      Author: bjfuwangzhu
 */

#include<stdio.h>
#include<string.h>
#include<math.h>
#define LL long long
#define nmod 1000000007
#define nmax 10010
#define nnum 1500
int prime[nnum], flag[nmax], pfactor[nnum];
int plen, pflen;
LL inverse, x, y;
void extend_gcd(LL a, LL b) {
    if (b == 0) {
        x = 1, y = 0;
        return;
    }
    extend_gcd(b, a % b);
    LL tx = x;
    x = y, y = tx - a / b * y;
}
void init() {
    int i, j;
    memset(flag, -1, sizeof(flag));
    for (i = 2, plen = 0; i < nmax; i++) {
        if (flag[i]) {
            prime[plen++] = i;
        }
        for (j = 0; (j < plen) && (i * prime[j] < nmax); j++) {
            flag[i * prime[j]] = 0;
            if (i % prime[j] == 0) {
                break;
            }
        }
    }
    extend_gcd(30, nmod);
    inverse = (x % nmod + nmod) % nmod;
}
LL getSumPow(int n) {
    LL res, temp;
    res = 1;
    res = res * n % nmod;
    res = res * (n + 1) % nmod;
    res = res * (n * 2 + 1) % nmod;
    temp = (LL) n;
    temp = ((temp * temp % nmod * 3 + temp * 3 % nmod - 1) % nmod + nmod) % nmod;
    res = res * temp % nmod;
    res = res * inverse % nmod;
    return res;
}
LL getPow(int n) {
    LL res;
    res = (LL) n;
    res = res * res % nmod * res % nmod * res % nmod;
    return res;
}
LL dfs(int start, int n) {
    LL res;
    int i, temp;
    for (i = start, res = 0; i < pflen; i++) {
        temp = pfactor[i];
        res = (res + getSumPow(n / temp) * getPow(temp) % nmod) % nmod;
        res = ((res - dfs(i + 1, n / temp) * getPow(temp) % nmod) % nmod + nmod)
                % nmod;
    }
    return res;
}
int main() {
#ifndef ONLINE_JUDGE
    freopen("data.in", "r", stdin);
#endif
    int t, n, N, i, te;
    LL res, temp;
    init();
    while (~scanf("%d", &t)) {
        while (t--) {
            scanf("%d", &n);
            te = (int) sqrt(n * 1.0);
            N = n;
            for (i = 0, pflen = 0; (i < plen) && (prime[i] <= te); i++) {
                if (n % prime[i] == 0) {
                    pfactor[pflen++] = prime[i];
                    while (n % prime[i] == 0) {
                        n /= prime[i];
                    }
                }
            }
            if (n > 1) {
                pfactor[pflen++] = n;
            }
            if (n == N) {
                printf("%lld\n", getSumPow(n - 1));
                continue;
            }
            temp = dfs(0, N);
            res = ((getSumPow(N) - temp) % nmod + nmod) % nmod;
            printf("%lld\n", res);
        }
    }
    return 0;
}

 以下转载于：

http://hi.baidu.com/tomspirit/blog/item/6b41263b4ef927ef3b87cecf.html

/*题意：找出比n小且与n互质的数的4次方和4: sum=1+3*3*3*3=82解法：
 * 容斥定理+状态压缩枚举
 * 比如12:sum = Sum(1^4+...12^4)%MOD - (所有非互质数的4次方和)%MOD
 * 求解Sum(1^4+...+12^4)有公式：1/30*n*(n+1)(2n+1)(3n^2+3n-1)，
 * 因为式子中有除，取余扩展欧几里得来做。现在的问题是所有非互质数的4次方和。
 * 可以先求12的素因子：2,3。在非互质数中，
 * ①含2素因子的数为2,4,6,8,10；含3素因子的数为3,6,9,12
 * ②含2,3的数为6,12可以看出非互质数 = ① - ②；
 * 满足容斥定理。对于2(含2因子)分别是2,4,6,8,10,12
 * ，转化为题意,2^4+4^4+6^4+8^4+10^4+12^4 = 2^4*(1^4+2^4+3^4+4^4+5^4+6^4)
 * 这样就可以利用上面的公式求解，但是这只是一种情况。我们可以用状态压缩来枚举所有的情况。
 * 对于12，有4中情况，全不选，只选2，只选3和2,3都选,
 * 那么相对于的和为0*(0)，2^4(1^4+...+6^4)，3^4(1^4+...+4^4)，6^4(1^4+...+2^4)
 * 这样才把所有的非互质数的4次方和找到，并利用容斥定理（奇数个素因子为+,偶数个素因子为-），
 * 答案就出来了。*/
#include<stdio.h>
#include<cstring>
#include<vector>
#include <map>
#include<iostream>
using namespace std;
const long long MOD = 1000000007;
const int N = 10010;
vector<long long> a;
bool f[N];
long long prim[N];
void Extend_euclid(long long a, long long b, long long &x, long long &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return;
    }
    long long tx, ty;
    Extend_euclid(b, a % b, tx, ty);
    x = ty;
    y = tx - a / b * x;
}
void getPrime() {
    long long i, j;
    memset(f, 1, sizeof(f));
    for (i = 2; i * i <= N; i++)
        if (f[i])
            for (j = i * i; j <= N; j += i)
                f[j] = 0;
    f[1] = 0;
    long long pt = 0;
    for (i = 2; i <= N; i++)
        if (f[i])
            prim[pt++] = i;
}
long long getCnt(long long nn) {
    long long res = 0, i;
    for (i = 0; prim[i] * prim[i] <= nn; ++i) {
        if (nn % prim[i] == 0) {
            a.push_back(prim[i]);
            nn /= prim[i];
            while (nn % prim[i] == 0) {
                nn /= prim[i];
            }
            if (nn == 1)
                break;
        }
    }
    if (nn > 1)
        a.push_back(nn);
    return a.size();
}
long long difNum, n, sum;
bool bit[10];
long long pow4(long long aa) {
    return aa * aa % MOD * aa % MOD * aa % MOD;
}
int main() {
    long long t, x, y, len, nn;
    int i, j;
    getPrime();
    Extend_euclid(30, MOD, x, y);
    x = (x % MOD + MOD) % MOD;
    scanf("%lld", &t);
    while (t--) {
        a.clear();
        scanf("%lld", &n);
        if (n == 1) {
            printf("0\n");
            continue;
        }
        nn = n;
        sum = nn * (nn + 1) % MOD * (2 * nn % MOD + 1) % MOD
                * (nn * nn % MOD * 3 % MOD + 3 * nn % MOD - 1) % MOD * x % MOD;
        difNum = getCnt(n);
        long long Limit = 1 << difNum;
        for (i = 0; i < Limit; ++i) {
            long long sum1 = 1, cnt = 0;
            for (j = 0; j < difNum; ++j) {
                if ((1 << j) & i) {
                    sum1 *= a[j];
                    ++cnt;
                }
            }
            if (cnt)
                nn = n / sum1;
            else
                nn = 0;
            if (cnt & 1)
                sum = (sum
                        - pow4(sum1)
                                * (nn * (nn + 1) % MOD * (2 * nn % MOD + 1)
                                        % MOD
                                        * (nn * nn % MOD * 3 % MOD
                                                + 3 * nn % MOD - 1) % MOD * x
                                        % MOD) % MOD + MOD) % MOD;
            else
                sum = (sum
                        + pow4(sum1)
                                * (nn * (nn + 1) % MOD * (2 * nn % MOD + 1)
                                        % MOD
                                        * (nn * nn % MOD * 3 % MOD
                                                + 3 * nn % MOD - 1) % MOD * x
                                        % MOD) % MOD + MOD) % MOD;
        }
        printf("%lld\n", sum);
    }
    return 0;
}