BZOJ 2627 JZPKIL

题面：

2627: JZPKIL

Time Limit: 50 Sec  Memory Limit: 128 MB
Submit: 347  Solved: 104
[Submit][Status][Discuss]

Description

Input

第一行，询问个数T。
　　下面T行，每行三个整数，n, x, y。

Output

　　T行，每行一个整数，表示相应的询问的答案

Sample Input

5
6 0 0
6 0 1
6 1 0
6 1 1
1000000000 50 50

Sample Output

6
66
15
126
393442025

HINT

数据规模和约定
30%的数据，x=y
另30%的数据，n<=10^9, x, y<=100
100%的数据，T<=100, 1<=n<=10^18, 0<=x, y<=3000

这题是个50sec卡评测反演好题。

$quad sum_{i=1}^{n}(i,n)^{x}[i,n]^{y}$

$=sum_{i=1}^{n}(i,n)^{x-y}(icdot n)^{y}$

$=n^{y}sum_{i=1}^{n}(i,n)^{x-y}i^{y}$

$=n^{y}sum_{d|n}d^{x}sum_{i=1}^{frac{n}{d}}sum_{k|i,k|frac{n}{d}}mu(k)i^{y}$

$=n^{y}sum_{d|n}d^{x}sum_{k|frac{n}{d}}mu(k)k^{y}sum_{i=1}^{frac{n}{kcdot d}i^{y}}$

$=n^{y}sum_{d|n}d^{x}sum_{k|frac{n}{d}}mu(k)k^{y}sum_{i=0}^{y}frac{1}{y+1}dbinom{y+1}{i}B_{i}(frac{n}{kcdot d})^{y+1-i}$

$=sum_{i=0}^{y}t_{i}n^{y}sum_{d|n}d^{x}sum_{k|frac{n}{d}}mu(k)k^{y}(frac{n}{kcdot d})^{y+1-i}$

$t_{i}=frac{1}{y+1}dbinom{y+1}{i}B_{i}$ ($B_{i}$为伯努利数)

将$n$分解质因数$n=p_{1}^{k_1}cdot p_{2}^{k_2}cdot p_{3}^{k_3}cdots $

令$f(p_{j}^{k_j},i)=sum_{d|p_{j}^{k_j}}d^{p_{j}^{k_j}}sum_{k|frac{p_{j}^{k_j}}{d}}mu(k)k^{y}(frac{p_{j}^{k_j}}{kcdot d})^{y+1-i}$

$f$为积性函数,可以用每个质因子的答案相乘得到$n$的答案。

但是由于$n$很大,我们只能用$Pollard_rho$分解大数因子。

时间复杂度为$O(y^{2}+Tcdot n^{frac{1}{4}}+Tcdot ycdot log;n)$

#include <iostream>
#include <stdio.h>
#include <algorithm>
using namespace std;
#define LL long long
template<typename __>
inline void read(__ &s)
{
    char ch=getchar();
    for(s=0;ch<'0'||ch>'9';ch=getchar());
    for(;ch>='0'&&ch<='9';s=s*10+ch-'0',ch=getchar());
}
namespace Pollard_rho
{
    inline LL mul(LL x,LL y,LL mod)
    {
        return (x*y-(LL)(x/(long double)mod*y+0.5)*mod+mod)%mod;
    }
    inline LL qpow(LL x,LL y,LL mod)
    {
        LL ans=1;
        for(;y;y>>=1,x=mul(x,x,mod))
            if(y&1)
                ans=mul(ans,x,mod);
        return ans;
    }
    inline bool judge(LL n,LL p)
    {
        if(n<=p)
            return 1;
        if(n%p==0)
            return 0;
        LL x=n-1;
        int s=-1;
        while(!(x&1))
            x>>=1,++s;
        x=qpow(p,x,n);
        if(x==1||x==n-1)
            return 1;
        while(s--)
        {
            x=mul(x,x,n);
            if(x==1)
                return 0;
            if(x==n-1)
                return 1;
        }
        return 0;
    }
    inline bool miller_rabin(LL n)
    {
        return judge(n,2)&&judge(n,3)&&judge(n,5)&&judge(n,7)&&judge(n,11)&&judge(n,13)&&judge(n,17)&&judge(n,19)&&judge(n,23);
    }
    LL *P;
    int Pcnt;
    inline LL pollard_rho(LL n)
    {
        LL x=rand()%(n-1)+1,xx=mul(x,x,n)+1,t=1;
        while(t==1)
        {
            x=mul(x,x,n)+1;
            xx=mul(xx,xx,n)+1;
            xx=mul(xx,xx,n)+1;
            t=__gcd(x-xx+n,n);
        }
        return t;
    }
    inline void divide(LL n)
    {
        if(n==1)
            return ;
        else
            if(n%2==0)
                P[++Pcnt]=2,divide(n/2);
            else
                if(miller_rabin(n))
                    P[++Pcnt]=n;
                else
                {
                    LL x=pollard_rho(n);
                    divide(x);
                    divide(n/x);
                }
    }
    inline void init(LL n,LL *p,int &cnt)
    {
        P=p;
        Pcnt=0;
        divide(n);
        cnt=Pcnt;
    }
}
#define mod 1000000007
inline int mul(int x,int y)
{
    return ((LL)x*y)%mod;
}
inline int add(int x,int y)
{
    x+=y;
    if(x>=mod)
        x-=mod;
    return x;
}
inline int dec(int x,int y)
{
    x-=y;
    if(x<0)
        x+=mod;
    return x;
}
inline int qpow(int x,int y)
{
    x%=mod;
    int ans=1;
    for(;y;y>>=1,x=mul(x,x))
        if(y&1)
            ans=mul(ans,x);
    return ans;
}
#define maxn 3002
int B[maxn],C[maxn][maxn],a[maxn],pw[20][maxn],iv[20][maxn];
int p[20],c[20],x,y,N;
LL n;
int tot;
inline void init()
{
    int i,j;
    C[0][0]=1;
    for(i=1;i<=3001;i++)
    {
        C[i][0]=1;
        for(j=1;j<=i;j++)
            C[i][j]=add(C[i-1][j],C[i-1][j-1]);  
    }
    B[0]=1;
    for(i=1;i<=3001;i++)
    {
        for(j=0;j<i;++j)
            B[i]=dec(B[i],mul(C[i+1][j],B[j]));
        B[i]=mul(B[i],qpow(i+1,mod-2));
    }
}
inline void pre()
{
    for(int i=0;i<=y+1;i++)
        a[i]=0;
    for(int i=0;i<=y;i++)
        a[y+1-i]=mul(C[y+1][i],B[i]);
    int tmp=qpow(y+1,mod-2);
    for(int i=0;i<=y+1;i++)
        a[i]=mul(a[i],tmp);
    a[y]=add(a[y],1);
    if(!y)
        a[0]=dec(a[0],1);
}
inline void get_frac()
{
    static LL frac[97];
    int cnt=0;
    tot=0;
    Pollard_rho::init(n,frac,cnt);
    sort(frac+1,frac+cnt+1);
    for(int i=1;i<=cnt;i++)
        if(frac[i]!=frac[i-1])
            p[++tot]=frac[i]%mod,c[tot]=1;
        else
            ++c[tot];
    for(int i=1;i<=tot;i++)
    {
        pw[i][0]=iv[i][0]=1;
        int t=max(max(x,y),c[i])+1;
        for(int j=1;j<=t;j++)
            pw[i][j]=mul(pw[i][j-1],p[i]);
        int v=qpow(p[i],mod-2);
        for(int j=1;j<=t;j++)
            iv[i][j]=mul(iv[i][j-1],v);
    }
}
inline int geo(int i,int x)
{
    int c=::c[i];
    int p;
    if(x<0)
        p=iv[i][-x];
    else
        p=pw[i][x];
    if(p==1)
        return c;
    int ans=0,t=p;
    for(int j=1;j<=c;j++)
    {
        ans=add(ans,t);
        t=mul(t,p);
    }
    return ans;
}
inline void input()
{
    read(n);
    read(x);
    read(y);
    unsigned int rand_seed=18230742;
    rand_seed+=(n+(x<<16)+y)+(rand_seed<<2);
    srand(rand_seed);
}
inline int calc(int x,int y)
{
    int ans=x<0?qpow(qpow(N,mod-2),-x):qpow(N,x);
    for(int i=1;i<=tot;i++)
    {
        int t1=geo(i,-x);
        int t2=dec(1,y>=0?pw[i][y]:qpow(p[i],mod-2));
        ans=mul(ans,add(1,mul(t1,t2)));
    }
    return ans;
}
void solve()
{
    pre();
    get_frac();
    N=n%mod;
    if(!N)
        return (void)(puts("0"));
    int ans=0;
    int t=1,tmp;
    for(int i=0;i<=y+1;i++)
    {
        tmp=calc(x-i,y-i);
        ans=add(ans,mul(mul(a[i],t),tmp));
        t=mul(t,N);
    }
    ans=mul(ans,qpow(N,y));
    printf("%d
",ans);
}
int main()
{
    init();
    int T;
    read(T);
    while(T--)
    {
        input();
        solve();
    }
}

PS：$Pollard_rho$写丑了会TLE,需要大力卡一波常数才能过。