HDU 6137 Engineering of the Clones（快速幂+NTT）

Description

给出一正整数A的质因子分解形式A=pr11pr22...prnn，其中|ri−rj|≤1，问1A−1在∏i=1npi进制下小数点后第k位的值

Input

首先输入一整数T表示用例组数，每组用例首先输入两整数n和k，然后输入n个互不相同的素数pi，之后输入n个整数ri

(1≤n≤106,1≤k<max{r21,r22,...,r2n},1≤∏i=1npi≤106,1≤ri≤5⋅104)

Output

输出1A−1在∏i=1npi进制下小数点后第k位的值

Sample Input

1

2 1

2 5

2 2

Sample Output

0

Solution

1A−1=−1+11−1A=−1+∑i=0∞(1A)i=∑i=1∞(1A)i

记L=max{r1,...,rn}，则1A在base=∏i=1npi进制下可以表示为一个L位小数，第L位为x=pL−r11...pL−rnn

由于1≤k<L2，(1A)L+1=xL+1base−L2−L<base−L2+1≤base−k，故在算第k位数时只需要算xm即可，其中m=⌊kL⌋+1≤L，而不用担心∑i=L+1∞(1A)i产生的进位影响，计算xm用快速幂加NTT即可

Code

#include<cstdio>#include<iostream>#include<cstring>#include<algorithm>#include<cmath>#include<vector>#include<queue>#include<map>#include<set>#include<ctime>using namespace std;typedef long long ll;typedef long double ld;const int maxbit=17,maxlen=1<<maxbit,maxn=100005;const ll mod=4179340454199820289,g=3;ll wn[maxlen],inv2[maxbit+1];ll mul(ll a,ll b){    return (a*b-(ll)(a/(ld)mod*b+1e-3)*mod+mod)%mod;}ll mod_pow(ll a,ll b){    ll ans=1;    while(b)    {        if(b&1)ans=mul(ans,a);        a=mul(a,a);        b>>=1;    }    return ans;}void init(){    wn[0]=1,wn[1]=mod_pow(g,(mod-1)>>maxbit);    for(int i=2;i<maxlen;i++)wn[i]=mul(wn[i-1],wn[1]);    inv2[0]=1,inv2[1]=(mod+1)/2;    for(int i=2;i<=maxbit;i++)inv2[i]=mul(inv2[i-1],inv2[1]);//预处理2^i的逆元 }void ntt(ll *x,int len,int sta) {    for(int i=0,j=0;i<len;i++)    {        if(i>j)swap(x[i],x[j]);        for(int l=len>>1;(j^=l)<l;l>>=1);    }    for(int i=1,d=1;d<len;i++,d<<=1)        for(int j=0;j<len;j+=d<<1)            for(int k=0;k<d;k++)            {                ll t=mul(wn[(maxlen>>i)*k],x[j+k+d]);                x[j+d+k]=x[j+k]-t<0?x[j+k]-t+mod:x[j+k]-t;                x[j+k]=x[j+k]+t>=mod?x[j+k]+t-mod:x[j+k]+t;            }    if(sta==-1)    {        reverse(x+1,x+len);        int bitlen=0;        while((1<<bitlen)<len)bitlen++;        ll val=inv2[bitlen];        for(int i=0;i<len;i++)x[i]=mul(x[i],val);    }}void Deal(ll *x,int len,int base){    for(int i=0;i<len;i++)        if(x[i]>=base)            x[i+1]+=x[i]/base,x[i]%=base;}int T,n,p[maxn],r[maxn];ll k,A[maxlen],B[maxlen],Ans[maxlen];ll Solve(int L,int m,int pos,int x,int base){    init();    memset(A,0,sizeof(A));    memset(Ans,0,sizeof(Ans));    A[0]=x,Ans[0]=1;    int len=2;    while(m)    {        ntt(A,len,1);        if(m&1)        {            ntt(Ans,len,1);            for(int i=0;i<len;i++)Ans[i]=mul(Ans[i],A[i]);            ntt(Ans,len,-1);            Deal(Ans,len,base);        }        for(int i=0;i<len;i++)A[i]=mul(A[i],A[i]);        ntt(A,len,-1);        Deal(A,len,base);        len<<=1;        m>>=1;    }    return Ans[pos];}int main(){    scanf("%d",&T);    while(T--)    {        scanf("%d%I64d",&n,&k);        int base=1;        for(int i=1;i<=n;i++)        {            scanf("%d",&p[i]);            base*=p[i];        }        for(int i=1;i<=n;i++)scanf("%d",&r[i]);        int L=r[1];        for(int i=1;i<=n;i++)L=max(L,r[i]);        int x=1;        for(int i=1;i<=n;i++)            if(r[i]<L)x*=p[i];        k--;        printf("%I64d\n",Solve(L,k/L+1,L-1-k%L,x,base));    }    return 0;}