hdu4059 容斥原理的好题

1.首先要知道4次方求和的公式

2.注意求和公式中有除以30应该乘以30的逆

3.容斥原理

即 原来的所有和 -  1个质因数的倍数 + 2个质因数的倍数 - 3个质因数的倍数 + 4个质因数的倍数……

每个质因数的贡献就是getsum(n/link[i]) * four(link[i])，

比如质因数是2时，就是 (1^4+2^4+3^4+4^4......+(n/2)^4) * (2^4)

还要减掉 dfs(i+1,n/link[i]) * four(link[i]) ，乘以的数是four(link[i])

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <set>
#include <algorithm>
#include <ctime>
#include <vector>
#include <string>
#include <stack>
#include <queue>
using namespace std;
/*ll ext_gcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1;y=0;
        return a;
    }
    ll x1,y1,d;
    d=ext_gcd(b,a%b,x1,y1);
    x=y1;
    y=x1-a/b*x;
    return d;
}
ll linear_equations(ll a,ll b,ll n) // a=b(mod n)
{
    ll x,y,d;
    d=ext_gcd(a,n,x,y);
    return (x%mod+mod)%mod;
}*/
#define mod 1000000007
#define ll long long
const int inv=233333335;
#define N 20000
int prime[N],is_prime[N],prime_num;
int link[100],cnt;
void init()
{
    for(int i=2;i<N;++i)
        if(!is_prime[i])
    {
        prime[prime_num++]=i;
        for(int j=i;j<N;j+=i)
            is_prime[j]=1;
    }
}
void get_divisor(ll n)
{
    cnt=0;
    for(int i=0;prime[i]*prime[i]<=n;++i)
        if(n%prime[i]==0)
        {
            link[cnt++]=prime[i];
            while(n%prime[i]==0)
                n/=prime[i];
        }
    if(n>1)
        link[cnt++]=n;
}
ll getsum(ll n)
{
    ll res=1;
    res*=n*(n+1)%mod;
    res*=(6*n*n%mod*n%mod+9*n*n%mod+n-1)%mod;
    res%=mod;
    res*=inv;
    res%=mod;
    return res;
}
ll four(ll n)
{
    ll a=n*n%mod;
    return a*a%mod;
}

ll dfs(int idx,ll n)
{
    ll res=0,temp;
    for(int i=idx;i<cnt;++i)
    {
        res=((res+(temp=getsum(n/link[i])*four(link[i])%mod))%mod+mod)%mod;
        res=((res-(temp=dfs(i+1,n/link[i])*four(link[i])%mod))%mod+mod)%mod;
    }
    return (res%mod+mod)%mod;
}
ll n;
int main ()
{
    int ncase;
    scanf("%d",&ncase);
    init();
    while(ncase--)
    {
        scanf("%I64d",&n);
        get_divisor(n);
        ll ans=getsum(n);
        ll t=dfs(0,n);
        printf("%I64d\n",((ans-t)%mod+mod)%mod);
    }
    return 0;
}