JZPKIL-CSDN博客

本文链接：https://blog.csdn.net/AcrossTheSky/article/details/18901741

终于A掉JZPKIL了。。感觉以后可以当数论题模板了(Bernoulli数、Pollard's rho、Miller Rabin、*快速乘、快速幂、逆元、组合数)

最神奇的是gyz的O(1)快速乘

证明见:http://tieba.baidu.com/p/2175344676

刷到R1了好开心

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <set>
#include <cmath>
#include <cstdlib>
#define mod 1000000007
#define abs(x) ((x)>0?(x):(-(x)))
using namespace std;
typedef long long LL;
LL p[20000],prime[20000];
int num[20000],tot;
LL P[20000],Q[20000],B[20000];
LL b[3010];
int t[3010][3010],c[3010][3010];
int cnt;
LL mul(LL x,LL y,LL z){return (x*y-(LL)(x/(long double)z*y+1e-3)*z+z)%z;}
inline LL power(LL p,LL n,LL mo)
{
	LL ans=1;
	for(;n;n>>=1,p=mul(p,p,mo))
		if(n&1)
			ans=mul(ans,p,mo);
	return ans;
}
LL power1(LL a,LL b,LL mo)
{
    if (b==0) return 1;
    if (b==1) return a%mo;
    LL m=power1(a,b/2,mo);
    m=(m*m)%mo;
    if (b&1) m=(m*a)%mo;
    return m;
}
inline bool Miller_Rabin(int p,LL n)
{
	int s=0;
	LL d=n-1;
	while(!(d&1))
		++s,d>>=1;
	LL w=power(p,d,n);
	if(w==1)
		return true;
	for(int i=0;i<=s-1;++i)
	{
		if(w==n-1)
			return true;
		w=mul(w,w,n);
	}
	return false;
}
inline bool isPrime(LL n)
{
	static int nd[]={2,3,5,7,11,13,17,19,23};
	if(n==1)
		return false;
	if(find(nd,nd+9,n)!=nd+9)
		return true;
	if(!(n&1))
		return false;
	for(int i=0;i<9;++i)
		if(!Miller_Rabin(nd[i],n))
			return false;
	return true;
}
LL gcd(LL a,LL b)
{
    if (b==0) return a;
    else return gcd(b,a%b);
}
void frac(LL n)
{
	if(isPrime(n))
	{
		p[++cnt]=n;
		return ;
	}
	int c=16831;
	while(1)
	{
		LL x1,x2;
		int i=1,k=2;
		x1=x2=1;
		while(1)
		{
			x1=mul(x1,x1,n)+c;
			x1%=n;
			LL d=gcd(abs(x1-x2),n);
			if(d!=1&&d!=n)
			{
				frac(d),frac(n/d);
				return ;
			}
			if(x1==x2)
				break;
			if(++i==k)
				k<<=1,x2=x1;
		}
		c--;
	}
}
LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if (b==0)
    {
        x=1,y=0;
        return a;
    }
    LL r=exgcd(b,a%b,x,y);
    LL t=x;
    x=y;
    y=t-a/b*y;
    return r;
}
LL anti(LL n)
{
    LL x,y;
    exgcd(n,mod,x,y);
    x=x%mod;
    if (x<0) x+=mod;
    return x;
}
void init()
{
    int n=3000;
    c[0][0]=1;
    for (int i=1;i<=n+1;++i)
    {
        c[i][0]=1;
        for (int j=1;j<=i;++j)
        {
            c[i][j]=c[i-1][j]+c[i-1][j-1];
            if (c[i][j]>=mod) c[i][j]-=mod;
        }
    }
    for (int i=0;i<=n;++i)
    {
        b[i]=i+1;
        for (int j=0;j<i;++j)
        {
            b[i]=(b[i]-b[j]*c[i+1][j])%mod;
            if (b[i]<0) b[i]+=mod;
        }
        b[i]=(b[i]*anti(c[i+1][i]))%mod;
    }
    for (int i=0;i<=n;++i)
    {
        LL a=anti(i+1);
        for (int j=0;j<=i;++j)
            t[i][j]=(a*c[i+1][j]%mod*b[j])%mod;
    }
}
void add(LL n,LL m)
{
    prime[++tot]=m;
    num[tot]=0;
    while (n%m==0)
    {
        n/=m;
        num[tot]++;
    }
}
LL kil[10000],jzp[10000],m[10000],super[100][100];
LL calc(LL n,LL k,int x,int y)
{
    LL res=1LL;
    for (int i=1;i<=tot;++i)
    {
        LL d=1LL,tmp=0;
        LL kk=1;
        for (int j=0;j<=num[i];++j)
        {
            LL a=super[i][num[i]-j];
            if (num[i]-j>=1) a-=(kil[i]*super[i][num[i]-j-1])%mod;
            a%=mod;
            if (a<0) a+=mod;
            tmp+=kk*a;
            kk*=jzp[i];
            kk%=mod;
            tmp%=mod;
            d*=prime[i];
        }
        res=res*tmp;
        res%=mod;
    }
    return res;
}
LL work(LL n,int x,int y)
{
    cnt=0,tot=0;
    frac(n);
    sort(p+1,p+cnt+1);
    add(n,p[1]);
    for (int i=2;i<=cnt;++i)
        if (p[i]!=p[i-1]) add(n,p[i]);
    LL ans=0;
	for (int i=1;i<=tot;++i)
	{
        kil[i]=power1(prime[i]%mod,y,mod);
        jzp[i]=power1(prime[i]%mod,x,mod);
        m[i]=power1(prime[i],num[i],n+1);
		LL d=1;
		for (int j=0;j<=num[i];++j,d*=prime[i]) super[i][j]=d%mod;
	}
    for (int i=y;i>=0;--i)
    {
        ans+=(t[y][i]*calc(n,i,x,y))%mod;
        ans%=mod;
        if (ans<0) ans+=mod;
		for (int i=1;i<=tot;++i)
		{
			LL d=1;
			for (int j=0;j<=num[i];++j,d*=prime[i]) super[i][j]=(super[i][j]*(d%mod))%mod;
		}
    }
	n%=mod;
    ans=(ans*power1(n,y,mod))%mod;
    return ans;
}
LL kill(LL n,int x,int y)
{
    LL ans=0,p=n%mod,d=p;
    for(int i=y;i>=0;--i)
    {
        ans=ans+t[y][i]*d;
        if(ans>=mod)
            ans%=mod;
		d=d*p;
		if(d>=mod)
			d%=mod;
    }
    return ans*power1(p,y,mod)%mod;
}
int main()
{
    int T;LL n;
    int x,y;
    init();
    cin>>T;
    while (T--)
    {
        cin>>n>>x>>y;
		//frac(n);
        if (x==y) cout<<kill(n,x,y)<<endl;
        else cout<<work(n,x,y)<<endl;
    }
    return 0;
}