[莫比乌斯反演 伯努利数] BZOJ 2627 JZPKIL

ydc的题解 orzz
因为这道题才去学习了伯努利数 可以看这里

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

namespace RHO{
  int prime[9]={2,3,5,7,11,13,17,19,23};
  unsigned long long RR;
  int Pcnt;
  long long p[70];
  inline long long R(long long Mo){return (RR+=4179340454199820289ll)%Mo;}
  inline long long mul(long long x,long long y,long long Mo){
    long long tmp=(x*y-(long long)((long double)x/Mo*y+0.5)*Mo);
    return tmp<0?tmp+Mo:tmp;
  }
  inline long long Pow(long long x,long long k,long long Mo){
    long long ans=1;
    for(;k;k>>=1,x=mul(x,x,Mo))if(k&1)ans=mul(ans,x,Mo);
    return ans;
  }
  inline bool MR(long long n){
    if(n<=1)return false;
    for(int i=0;i<9;i++)if(n==prime[i])return true;
    long long d=n-1;int tmp=0;
    while((d&1)==0)d>>=1,tmp++;
    for(int i=0;i<9;i++){
      long long x=Pow(prime[i],d,n),p=x;
      for(int j=1;j<=tmp;j++){
    x=mul(x,x,n);
    if((x==1)&&(p!=1)&&(p!=n-1))return false;
    p=x;
      }
      if(x!=1)return false;
    }
    return true;
  }
  inline long long f(long long x,long long c,long long Mo){return (mul(x,x,Mo)+c)%Mo;}
  inline long long gcd(long long x,long long y){return x==0?y:gcd(y%x,x);}
  inline long long check(long long c,long long n){
    long long x=R(n),y=f(x,c,n),p=n;
    while((x!=y)&&((p==n)||(p==1))){
      if(x>y)p=gcd(n,x-y);
      else p=gcd(n,y-x);
      x=f(x,c,n);y=f(f(y,c,n),c,n);
    }
    return p;
  }
  inline void rho(long long n){
    if(n<=1)return;
    if(MR(n)){p[++Pcnt]=n;return;}
    while(true){
      long long tmp=check(R(n-1)+1,n);
      if(tmp!=n && tmp!=1){rho(tmp),rho(n/tmp);return;}
    }
  }
}

const int M=3005;
const int P=1e9+7;
int C[M][M],inv[M],B[M];
int T[M];

inline void Pre(int m){
  inv[1]=1; for (int i=2;i<=m+1;i++) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
  C[0][0]=1;
  for (int i=1;i<=m+1;i++){
    C[i][0]=1; for (int j=1;j<=i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
  }
  B[0]=1;
  for (int i=1;i<=m;i++){
    ll ret=0;
    for (int j=0;j<i;j++) ret+=(ll)C[i+1][j]*B[j]%P;
    ret%=P;
    B[i]=(P-ret*inv[i+1]%P)%P;
  }
}

ll n; int x,y;
int Pcnt; ll p[70]; int q[70];

inline ll Pow(ll a,int b){
  a%=P; ll ret=1;
  for (;b;b>>=1,a=a*a%P) if (b&1) ret=ret*a%P;
  return ret;
}
inline ll Inv(ll a){ return Pow(a,P-2); }

inline ll calc(ll p,int t,int x,int k){
  if (t==1) return (Pow(p,x)+Pow(p,k))%P;
  if (t==0) return 1;
  if (x==k) return Pow(p,t*k)*(t+1)%P;
  ll tt=x>k?Pow(p,x-k):Inv(Pow(p,k-x));
  return Pow(p,t*k)*(Pow(tt,t+1)-1)%P*Inv(tt-1)%P;
  /*ll ret=0;
  for (int a=0;a<=t;a++)
    ret+=Pow(p,a*x+(t-a)*k);
    return ret%P;*/
}

inline ll Solve(int K){
  ll ret=1;
  for (int i=1;i<=Pcnt;i++){
    ll tmp=(calc(p[i],q[i],x,K)+P-Pow(p[i],y)*calc(p[i],q[i]-1,x,K)%P)%P;
    ret=ret*tmp%P;
  }
  return ret;
}

int main(){
  int T;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  Pre(3000); scanf("%d",&T);
  while (T--){
    scanf("%lld%d%d",&n,&x,&y);
    RHO::Pcnt=0; RHO::rho(n);
    Pcnt=RHO::Pcnt; for (int i=1;i<=Pcnt;i++) p[i]=RHO::p[i];
    ll tmp=n;
    sort(p+1,p+Pcnt+1); Pcnt=unique(p+1,p+Pcnt+1)-p-1;
    for (int i=1;i<=Pcnt;i++){
      q[i]=0;
      while (tmp%p[i]==0) tmp/=p[i],q[i]++;
    }
    ll ans=0;
    for (int k=1;k<=y+1;k++){
      ll tk=(ll)inv[y+1]*C[y+1][k]%P*B[y+1-k]%P;
      if (k==y) tk=(P-tk)%P;
      ans+=tk*Solve(k)%P;
    }
    ans%=P; ans=ans*Pow(n,y)%P;
    printf("%lld\n",ans);
  }
  return 0;
}