hdoj 4059 The Boss on Mars

类型：数学

题目：http://acm.hdu.edu.cn/showproblem.php?pid=4059

来源：2011 Asia Dalian Regional Contest

思路：参考：http://blog.sina.com.cn/s/blog_69c3f0410100vqfe.html

四次方和公式：

f(n) = 1^4 + 2^4 + … + n^4 = n * (n + 1) * (2 * n + 1) * (3 * n^2 + 3 * n - 1) / 30

答案ans为f(n)与所有与g(n)【g(n)为与n不互质且不大于n的正整数四次方和】的差

g(n)的求法：容斥原理。g(n)为所有质因数中任意奇数个质因数的积的倍数【不大于n的所有倍数】的四次方和与任意偶数个质因数的积的倍数的四次方和的差

！！！状态压缩枚举所有不同素因子的组合

！！！(a / b) % x = (a % (x * b)) / b

说明：(a / b) % x 
= a / b - floor((a / b) / x) * x 
= (a - floor((a / b) / x) * x * b) / b
= (a - floor(a / (b * x)) * x * b) / b
= (a % (x * b)) / b

// hdoj 4059 The Boss on Mars
// tle ac 62MS	344K
#include <iostream>
#include <string>
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
using namespace std;

#define FOR(i,a,b) for(i = (a); i < (b); ++i)
#define FORE(i,a,b) for(i = (a); i <= (b); ++i)
#define FORDE(i,a,b) for(i = (a); i >= (b); --i)
#define CLR(a,b) memset(a,b,sizeof(a))

const int MAXN = 10010;
const int INF = 0x7f7f7f7f;

__int64 n, ans, L, x, s, cas;
__int64 mod = 1000000007;
__int64 a[MAXN], prime[MAXN];

void get_prime() {
    int i, j;
    FOR(i, 2, MAXN)
        a[i] = 1;
    FOR(i, 2, MAXN)
        if(a[i]) {
            prime[L++]=i;
            for(j = i + i; j < MAXN; j += i)
                a[j] = 0;
        }
}

// a * b % x
__int64 mult(__int64 a, __int64 b, __int64 x) {
    __int64 s = 0;
    while(b) {
        if(b & 1)
            s = (s + a) % x;
        a <<= 1;
        b >>= 1;
        if(a > n)
            a = a % x;
    }
    return s;
}

//(a / b) % x = (a % (x * b)) / b
__int64 fact_4(__int64 x) {
    __int64 s, m;
    m = mod * 30;
    s = mult(3 * x, x, m);
    s = (s + 3 * x - 1) % m;
    s = mult(s, 2 * x + 1, m);
    s = mult(s, x + 1, m);
    s = mult(s, x, m);
    return s / 30;
}

int main() {
    __int64 i, t, k, j, p, q;

    get_prime();
    scanf("%I64d", &cas);
    while(cas--) {
        scanf("%I64d", &n);
        k = n;
        i = t = 0;
        while(k > 1 && i < L) {
            if(k % prime[i] == 0) {
                while(k % prime[i] == 0)
                    k /= prime[i];
                a[t++] = prime[i];
            }
            i++;
        }
        if(k > 1)
            a[t++] = k;
        s = 1 << t;
        ans = 0;
        FOR(i, 1, s) {
            k = i, p = 1, q = 0;
            FOR(j, 0, t - 1) {
                if(k & 1) {
                    p *= a[j];
                    q++;
                }
                k >>= 1;
            }
            if(q & 1)
                q = 1;
            else
                q = -1;
            x = mult(mult(p, p, mod), mult(p, p, mod), mod);
            ans = (ans + q * mult(x, fact_4(n / p), mod)) % mod;
        }
        ans=(fact_4(n) - ans + mod) % mod;
        if(n == 1)
            ans = 0;
        printf("%I64d\n", ans);
    }
    return 0;
}