hdu_6088 Rikka with Rock-paper-scissors

题目:

  As we know, Rikka is poor at math. Yuta is worrying about this situation, so he gives Rikka some math tasks to practice. There is one of them:

  Alice and Bob are going to play a famous game: Rock-paper-scissors. Both of them don’t like to think a lot, so both of them will use the random strategy: choose rock/paper/scissors in equal probability.

  They want to play this game times, then they will calculate the score s in the following way: if Alice wins times, Bob wins times, and in the remaining games they make a tie, the score will be the greatest common divisor of and .

  Know Yuta wants to know the expected value of . It is easy to find that the answer must be an integer.

  It is too difficult for Rikka. Can you help her?

  Note: If one of is , we define the greatest common divisor of and as .

思路:

  每一次胜负平的概率相等,所以每次两人赢了 把的方案数就是 。那么题目所求的就是:
  根据 ,把 作为 带入,得到,这样就把 去掉了。
  再交换这几个 ,把 提到最外面,得到 ,再把 都变成连续的数字:。这样后面那个分式就能看成两个次数分别为 的多项式相乘。这样复杂度就能做到 ,大约是 的。

代码:

#include <bits/stdc++.h>
#define mset(a,b) memset(a,b,sizeof a)
#define rep(i,a,b) for(int i(a),i##_END_(b);i<=i##_END_;i++)
#define drep(i,a,b) for(int i(a),i##_END_(b);i>=i##_END_;i--)
template<class T> inline bool tomax(T &a,T b){return a<b?a=b,1:0;}
template<class T> inline bool tomin(T &a,T b){return b<a?a=b,1:0;}
typedef long long ll;
const int N=262144+5,M=1e5+5;
const double PI=acos(-1);

int qkpow(int x,int k,int m){
    int ans=1;
    for(;k;k>>=1,x=(ll)x*x%m)
        if(k&1)ans=(ll)ans*x%m;
    return ans;
}
int inv(int a,int m){
    return qkpow(a,m-2,m);
}

namespace MTT{
    struct Comp{
        double r,i;
        Comp operator+ (const Comp _) const {
            return (Comp){r+_.r,i+_.i};
        }
        Comp operator- (const Comp _) const {
            return (Comp){r-_.r,i-_.i};
        }
        Comp operator* (const Comp _) const {
            return (Comp){r*_.r-i*_.i,r*_.i+i*_.r};
        }
    };
    int rev[N];
    void DFT(Comp a[],Comp w[],int n){
        rep(i,0,n-1)
            if(i<rev[i])
                std::swap(a[i],a[rev[i]]);
        for(int l=2,m=1;l<=n;l<<=1,m<<=1)
            for(Comp *p=a;p!=a+n;p+=l)
                rep(i,0,m-1){
                    Comp t=w[n/l*i]*p[i+m];
                    p[i+m]=p[i]-t;
                    p[i]=p[i]+t;
                }
    }
    Comp a[N],b[N],c[N],d[N],omg[N];
    int log2(int x){
        int r=0;
        while((1<<r)<x)r++;
        return r;
    }
    void calc(int _a[],int _b[],int res[],int n,int m,int mod){
        int lim=log2(n+m),r=1<<lim;
        rep(i,0,r-1)a[i]=b[i]=(Comp){0,0};
        rep(i,0,n-1)a[i]=(Comp){_a[i]&32767,_a[i]>>15};
        rep(i,0,m-1)b[i]=(Comp){_b[i]&32767,_b[i]>>15};
        rep(i,0,r-1){
            omg[i]=(Comp){cos(PI*2/r*i),sin(PI*2/r*i)};
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(lim-1));
        } 
        DFT(a,omg,r);DFT(b,omg,r);
        rep(i,0,r-1){
            int j=(r-i)&(r-1);
            c[j]=(Comp){(a[i].r+a[j].r)/2,(a[i].i-a[j].i)/2}*b[i];
            d[j]=(Comp){(a[i].i+a[j].i)/2,(a[j].r-a[i].r)/2}*b[i];
        }
        DFT(c,omg,r);DFT(d,omg,r);
        rep(i,0,n+m-2){
            ll u=c[i].r/r+0.5,v=c[i].i/r+0.5;
            ll x=d[i].r/r+0.5,y=d[i].i/r+0.5;
            u%=mod;v%=mod;x%=mod;y%=mod;
            res[i]=(u+(v+x<<15)+(y<<30))%mod;
        }
    }
}

int phi[M],pri[M];
bool np[M];
void getPhi(int n){
    int tot=0;
    phi[1]=1;
    rep(i,2,n){
        if(!np[i]){
            phi[i]=i-1;
            pri[++tot]=i;
        }
        rep(j,1,tot){
            if(i*pri[j]>n)break;
            np[i*pri[j]]=true;
            if(i%pri[j])
                phi[i*pri[j]]=phi[i]*phi[pri[j]];
            else {
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
        }
    }
}

int a[N],b[N],res[N],fac[M],ifac[M];
int main(){
    getPhi(1e5);
    int cas,n,m;
    scanf("%d",&cas);
    while(cas--){
        scanf("%d%d",&n,&m);

        fac[0]=1;
        rep(i,1,n)fac[i]=(ll)fac[i-1]*i%m;
        ifac[n]=inv(fac[n],m);
        drep(i,n-1,0)ifac[i]=(ll)ifac[i+1]*(i+1)%m;

        int ans=0;
        rep(d,1,n){
            int x=n/d;
            a[0]=b[0]=0;
            rep(i,0,x)
                a[i]=b[i]=ifac[i*d];
            MTT::calc(a,b,res,x+1,x+1,m);
            rep(i,1,x){
                int k=(ll)res[i]*fac[n]%m*ifac[n-i*d]%m;
                k=(ll)k*phi[d]%m;
                (ans+=k)%=m;
            }
        }
        ans=(ll)ans*qkpow(3,n,m)%m;
        printf("%d\n",ans);
    }
    return 0;
}