zju3547

题意：给出n(1<=n<=10^8)，求小于n的，求所有与n互质的数字的四次幂的加和是多少。

分析：容斥原理

首先要知道四次幂求和公式，1^4+2^4+...+n^4=n*(n+1)*(2n+1)*(3n^2+3n-1)/30

先求所有小于等于n的数字的四次幂和，然后减去那些不互质的即可。

这个减去的过程用到了容斥原理。

先对n分解质因子，每个不同的质因子只保留一个。

然后分别枚举这些质因子的组合情况，由奇数个因子组成的数要减去，由偶数个因子组成的数要加上。

对于一个因子组合的乘积a，我们需要一次性计算a^4+(2a)^4 + (3a)^4+...

将其转化为a^4 * (1^4+2^4+...)即可。

这道题还有一个难点，就是公式中有除法（除以30），却还要进行模运算。

除法是不支持模运算的，因此我们要将除法转化为乘法，除以30变为乘以30的逆元。

逆元的意思是，如果a、b互为mod c下的逆元，则a * b = 1 (mod c)。

求逆元可以用扩展欧几里德gcd(30,MOD,x,y)，把x/gcd(30,MOD)整理到0～MOD-1范围内即为30的逆元。

具体原因查阅扩展欧几里德算法。

#include <cstdio>
using namespace std;

#define D(x) 

const int MOD = (int)(1e9) + 7;
const int MAX_FACTOR = 40;

int n;
int factor_num;
long long factor[MAX_FACTOR];
long long inverse;

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

long long to_forth(long long value)
{
    long long ret = value;
    ret = ret * ret % MOD;
    ret = ret * ret % MOD;
    return ret;
}

long long cal(long long value)
{
    long long num = n / value;
    long long ret = 1;
    ret = ret * num % MOD * (num + 1) % MOD;
    ret = ret * (2 * num + 1) % MOD;
    ret = ret * ((num * num % MOD * 3 % MOD + 3 * num % MOD - 1) % MOD) % MOD;
    if (ret / 30 != ret * inverse % MOD)
    {
        D(printf("#%lld %lld
", ret / 30, ret * inverse % MOD));
    }else
    {
        D(printf("**
"));
    }
    ret = ret * inverse % MOD;

    ret = ret * to_forth(value) % MOD;

    return ret;
}

void get_factors()
{
    factor_num = 0;
    int m = n;
    for (int i = 2; i * i <= m; i++)
    {
        if (m % i == 0)
            factor[factor_num++] = i;
        while (m % i == 0)
        {
            m /= i;
        }
    }
    if (m != 1)
    {
        factor[factor_num++] = m;
    }
}

long long work()
{
    long long ans = 0;
    for (int i = 1; i < (1 << factor_num); i++)
    {
        int num = 0;
        long long temp = 1;
        int index = 0;
        for (int mask = 1; mask <= i; mask <<= 1, index++)
        {
            if ((mask & i) == 0)
            {
                continue;
            }
            num++;
            temp *= factor[index];
        }
        D(printf("temp=%lld
", temp));
        if (num & 1)
            ans += cal(temp);
        else
            ans -= cal(temp);
        ans = (ans % MOD + MOD) % MOD;
    }
    ans = ((cal(1) - ans) % MOD + MOD) % MOD;
    return ans;
}

void gcd_extend(long long a,long long b,long long &g,long long &x,long long &y)
{
    if (!b)
    {
        g = a;
        x = 1;
        y = 0;
        return;
    }
    gcd_extend(b, a % b, g, y, x);
    y -= a / b * x;
}

int main()
{
    long long x, y, g;
    gcd_extend(30, MOD, g, x, y);
    D(printf("%lld %lld %lld
", x, y, g));
    x = (x % MOD + MOD) % MOD;
    inverse = x / g;
    D(printf("%lld
", inverse));
    D(printf("%lld
", inverse * 30 % MOD));
    int t;
    scanf("%d", &t);
    while (t--)
    {
        scanf("%d", &n);
        if (n == 1)
        {
            puts("0");
            continue;
        }
        get_factors();
        int ans_int = work();
        printf("%d
", ans_int);
    }
    return 0;
}