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

Description

给出一正整数$A$的质因子分解形式$A=p_1^{r_1}p_2^{r_2}...p_n^{r_n}$，其中$|r_i-r_j|\le1$，问$\frac{1}{A-1}$在$\prod\limits_{i=1}^np_i$进制下小数点后第$k$位的值

Input

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

$(1\le n\le 10^6,1\le k<max\{r_1^2,r_2^2,...,r_n^2\},1\le \prod\limits_{i=1}^np_i\le 10^6,1\le r_i\le 5\cdot10^4)$

Output

输出$\frac{1}{A-1}$在$\prod\limits_{i=1}^np_i$进制下小数点后第$k$位的值

Sample Input

1

2 1

2 5

2 2

Sample Output

0

Solution

$\frac{1}{A-1}=-1+\frac{1}{1-\frac{1}{A}}=-1+\sum\limits_{i=0}^\infty(\frac{1}{A})^i=\sum\limits_{i=1}^\infty(\frac{1}{A})^i$

记$L=max\{r_1,...,r_n\}$，则$\frac{1}{A}$在$base=\prod\limits_{i=1}^np_i$进制下可以表示为一个$L$位小数，第$L$位为$x=p_1^{L-r_1}...p_n^{L-r_n}$

由于$1\le k<L^2$，$(\frac{1}{A})^{L+1}=x^{L+1}base^{-L^2-L}<base^{-L^2+1}\le base^{-k}$，故在算第$k$位数时只需要算$x^m$即可，其中$m=\lfloor \frac{k}{L}\rfloor+1\le L$，而不用担心$\sum\limits_{i=L+1}^\infty(\frac{1}{A})^i$产生的进位影响，计算$x^m$用快速幂加$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;
}