bzoj 4772 显而易见的数论——拆分数（五边形数定理）+线性筛

题目：https://www.lydsy.com/JudgeOnline/problem.php?id=4772

题解：https://blog.csdn.net/Dream_Lolita/article/details/82314788

关于 ( g[p^t] ) 的值是多少，提供自己的见解：

　　首先，和 ( p^t ) 互质的数有 ( p^{t-1} ) 段，每段有 ( p-1 ) 个，模 p 等于 1 ~ p-1 。

　　那么和 ( p^t ) 做 gcd 的就是 ( p^{t-1} ) 段模 p 等于 0 ~ p-2 的数。

　　把模 p 等于 1 ~ p-2 的数的贡献写在一起，就是 ( p^{t-1}*(p-2) ) ；

　　考虑剩下的那些 0 , p , 2p , 3p , ...... p^t-1*p ，那个 0 在这道题里可以写成 p^t 。

　　考虑 gcd 里贡献 p 的，有 ( p^{t-1} - p^{t-2} ) 个；贡献 p² 的，有 ( p^{t-2} - p^{t-3} ) 个，以此类推。

　　所以贡献就是 ( p*( p^{t-1} - p^{t-2} ) + p^2*( p^{t-2} - p^{t-3} ) + ... + p^{t-1}*( p^1 - p^0 ) + p^t ) 

　　乘开就是 ( ( p^t - p^{t-1} ) + ( p^t - p^{t-1} ) + ... + ( p^t - p^{t-1} ) + p^t = (t-1)*( p^t - p^{t-1} ) + p^t )

　　所以 ( g[p^t] = p^{t-1}*(p-2) + (t-1)*(p^t - p^{t-1}) + p^t )

　　如果想把 ( g[p^{t-1}] ) 代入式子里，就会变成：

　　　　( g[p^{t-1}]=p^{t-2} + (t-2)*(p^{t-1} - p^{t-2}) + p^t-1 )

　　　　( g[p^t] = p*g[p^{t-1}]+(p^t - p^{t-1}) )

筛的时候可以记录一下 mindiv 的 p 是 p 的几次方，就可以方便地知道 i 是不是 ( p^t ) 或者 i 是由哪个互质的数乘起来的了。

注意异或运算要加括号。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
int rdn()
{
  int ret=0;bool fx=1;char ch=getchar();
  while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return fx?ret:-ret;
}
int Mx(int a,int b){return a>b?a:b;}
int Mn(int a,int b){return a<b?a:b;}
const int N=2005,M=1e5+5,M2=1e7+5,mod=1e9+7;
int upt(int x,int md=mod){while(x>=md)x-=md;while(x<0)x+=md;return x;}
int pw(int x,int k,int md=mod)
{int ret=1;while(k){if(k&1)ret=(ll)ret*x%md;x=(ll)x*x%md;k>>=1;}return ret;}

int type,n,k,a[M],f[N][N],g[M2],pri[M2],mdv[M2];
int p[N],nm[N][N],ct[M],jc[N],jcn[N]; bool vis[M2];
int gcd(int a,int b){return b?gcd(b,a%b):a;}
bool flag=0;
int F(int x,int y)
{
  if(type==1)return 1;
  else if(type==2)return gcd(x,y);
  else return upt(pw(x,y,k)+pw(y,x,k)+(x^y),k);//(x^y) token!!!
}
int C(int n,int m){return (ll)jc[n]*jcn[m]%mod*jcn[n-m]%mod;}
void init(int mx)
{
  p[0]=1;
  for(int i=1;i<=n;i++)
    for(int j=1;;j++)
      {
    int k0=j*(3*j-1)>>1, k1=j*(3*j+1)>>1;
    int fx=(j&1)?1:-1;
    if(k0>i&&k1>i)break;
    if(k0<=i)p[i]=upt(p[i]+fx*p[i-k0]);
    if(k1<=i)p[i]=upt(p[i]+fx*p[i-k1]);
      }

  jc[0]=1;for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod;
  jcn[n]=pw(jc[n],mod-2,mod);
  for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;
  for(int i=1;i<=n;i++)
    for(int j=n-i;j>=i;j--)
      {
    int ret=0;
    if(j==i)
      {
        for(int k=2;k*i<=n;k++)
          {
        int tmp=upt(p[n-k*i]-((k+1)*i<=n?p[n-(k+1)*i]:0));
        ret=(ret+(ll)C(k,2)*tmp)%mod;
          }
      }
    else
      {
        for(int k0=1,d0=i;d0+j<=n;k0++,d0+=i)
          for(int k1=1,d1=j;d0+d1<=n;k1++,d1+=j)
        ret=upt(ret+p[n-d0-d1]);
      }
    int d=F(i,j)%k; ct[d]=upt(ct[d]+ret);
      }

  int cnt=0; g[1]=1;
  for(int i=2;i<=mx;i++)
    {
      if(!vis[i])g[i]=upt(2*i-2),pri[++cnt]=i,mdv[i]=i;
      for(int j=1,d;j<=cnt&&(d=i*pri[j])<=mx;j++)
    {
      vis[d]=1; int p=pri[j];
      if(i%pri[j]==0)
        {
          mdv[d]=mdv[i]*p;
          if(mdv[d]==d)g[d]=((ll)g[i]*p+d-i)%mod;
          else g[d]=(ll)g[d/mdv[d]]*g[mdv[d]]%mod;
          break;
        }
      g[d]=(ll)g[i]*g[p]%mod; mdv[d]=p;
    }
    }
}
int main()
{
  type=rdn();n=rdn();k=rdn(); int mx=0;
  for(int i=0;i<k;i++)a[i]=rdn(),mx=Mx(mx,a[i]);
  init(mx); int ans=0;
  for(int i=0;i<k;i++)
    ans=(ans+(ll)g[a[i]]*ct[i])%mod;
  printf("%d
",ans); return 0;
}