BZOJ 4772 显而易见的数论（生成函数+组合数学+数论）

Description

Input

Output

Sample Input
1
3 3
0 1 2
Sample Output
4
Solution


Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int maxn=2001,maxm=1e7+1,mod=1e9+7;
int p[maxn][maxn];//p[i][j]表示i^j%k
int gcd[maxn][maxn];//gcd[i][j]表示i和j的最大公倍数
int sum[maxn];//sum[i]表示i的拆分数
int num[maxn][maxn];//num[i][j]表示i和j在多少种拆分方案中存在
int res=0,prime[666666],temp[maxm];
bool mark[maxm];
//res统计素数个数,prime记录素数,mark标记素数,temp[i]记录i的最小素因子及其幂指数 
int cnt[maxm];//cnt[i]表示F(pi,pj)%k=i的方案数 
int type,n,k,g[maxm],a[maxm];
void inc(int &x,int y)
{
    x=x+y>=mod?x+y-mod:x+y;
}
void dec(int &x,int y)
{
    x=x-y<0?x-y+mod:x-y;
}
void init()
{
    for(int i=1;i<=n;i++)
    {
        p[i][0]=1;
        for(int j=1;j<=n;j++)p[i][j]=(ll)p[i][j-1]*i%k;
    }
    for(int i=1;i<=n;i++)gcd[0][i]=gcd[i][0]=gcd[i][i]=i,gcd[1][i]=gcd[i][1]=1;
    for(int i=2;i<=n;i++)
        for(int j=2;j<i;j++)
        {
            if(!gcd[i][j])gcd[i][j]=gcd[j][i-j];
            gcd[j][i]=gcd[i][j];
        }
    sum[0]=1;
    for(int i=1;i<=n;i++)
    {
        for(int j=1,w=1;w<=i;w+=3*j+1,j++)
            if(j&1)inc(sum[i],sum[i-w]);
            else dec(sum[i],sum[i-w]);
        for(int j=1,w=2;w<=i;w+=3*j+2,j++)
            if(j&1)inc(sum[i],sum[i-w]);
            else dec(sum[i],sum[i-w]);
    }
    g[0]=0,g[1]=1;
    for(int i=2;i<=1e7;i++)
    {
        if(!mark[i])prime[res++]=i,temp[i]=i,g[i]=2*i-2;
        for(int j=0;j<res&&i*prime[j]<=1e7;j++)
        {
            mark[i*prime[j]]=1;
            if(i%prime[j])
            {
                temp[i*prime[j]]=prime[j];
                g[i*prime[j]]=(ll)g[i]*g[prime[j]]%mod;
            }
            else
            {
                temp[i*prime[j]]=temp[i]*prime[j];
                if(temp[i]!=i)g[i*prime[j]]=(ll)g[i/temp[i]]*g[temp[i]*prime[j]]%mod;
                else g[i*prime[j]]=((ll)prime[j]*g[i]+i*prime[j]-i)%mod;
                break;
            }//g[p^k]=(k+1)(p^k-p^(k-1))
        } 
    }
} 
int f(int x,int y)
{
    if(type==1)return 1%k;
    if(type==2)return gcd[x][y]%k;
    if(type==3)return (p[x][y]+p[y][x]+(x^y))%k;
}
int main()
{
    scanf("%d%d%d",&type,&n,&k);
    for(int i=0;i<k;i++)scanf("%d",&a[i]);
    init();
    for(int i=1;i<=n;i++)
        for(int j=i+1;i+j<=n;j++)
        {
            int t=f(i,j);
            for(int ni=1;ni*i+j<=n;ni++)
                for(int nj=1;ni*i+nj*j<=n;nj++)
                    inc(cnt[t],sum[n-ni*i-nj*j]); 
        }
    for(int i=1;i<=n;i++)
    {
        int t=f(i,i);
        for(int ni=1;ni*i<=n;ni++)
        {
            int s=sum[n-ni*i];
            if((ni+1)*i<=n)dec(s,sum[n-(ni+1)*i]);
            inc(cnt[t],(ll)ni*(ni-1)/2*s%mod);
        }
    }
    int ans=0;
    for(int i=0;i<k;i++)inc(ans,(ll)cnt[i]*g[a[i]]%mod);
    printf("%d\n",ans);
    return 0;
}