bzoj 4961: 除除除

Description

我们定义一种操作是将一个正整数n(n>1)用某个整数d替换，这个数必须是n的约数(1≤d≤n)。给你一个正整数n，

你需要确定操作进行的期望次数，如果我们希望不断地使用这种操作来将n变成1，假设每次操作选择每个可能的d

的概率均等。为了便于计算，输入将给出n和它的所有不同质因数p_1,p_2,?p_m，保证n恰好有m个不同的质因数。

为了便于输出，设答案是有理数a/b，并且有bc≡1(mod10^9+7)，你只需要输出ac对10^9+7取模的值。例如当n=351

384000时，期望运算的次数为

1384855049944986283970053414177036273994739277918823/282971529872677632598150446595770345000925504317000≈4.893973081

但你只需要输出321468106即可。

Input

输入包含多组测试数据，以EOF结束。

对于每组测试数据：

第一行包含两个正整数n和m，其中m表示n的不同质因数个数，满足2≤n≤10^24。

第二行包含m个质数p_1,p_2,...,p_m，对于i=1,2,...,m满足2≤p_i≤10^6。

约200000组数据。

Output

对于每组测试数据，输出一行一个整数，表示题目要求输出的值。

对一个询问，答案只和每个质因子的幂次构成的可重集有关，在数据范围内只有约170000个本质不同的询问，因此可以预处理递推出答案

递推式为（f为答案，d为约数个数）：

$f(n)=\frac{d(n)+\sum_{i|n and i<n}f(d)}{d(n)-1}$

为了优化这个递推式，将n的质因子的幂次降序排列为p(n,1),p(n,2),p(n,3),...，求和部分考虑记录一个前缀和g(n,k)辅助计算，表示满足 当x<=k时p(a,x)=p(n,x)，否则p(a,x)<=p(n,x) 的f(a)之和，按p(n)的字典序升序处理，适当用hash存储g可以使转移复杂度与g的状态数（约1300000）同阶

#include<cstdio>
typedef unsigned int u32;
typedef unsigned long long u64;
u64 pp[107];
const u64 _eq=1844677;
const u32 P=1e9+7;
struct num{
    u64 v0,v1,v2;
    void push(u32 x){
        v0=v0*10+x,v1*=10,v2*=10;
        v1+=v0>>32,v0=u32(v0);
        v2+=v1>>32,v1=u32(v1);
    }
    bool div(u32 x){
        u64 a0=v0,a1=v1,a2=v2;
        a1+=a2%x<<32,a2/=x;
        a0+=a1%x<<32,a1/=x;
        if(a0%x)return 0;
        v0=a0/x,v1=a1,v2=a2;
        return 1;
    }
    bool read(){
        v0=v1=v2=0;
        int c=getchar();
        while(c<48){
            if(c<0)return 0;
            c=getchar();
        }
        while(c>47)push(c-48),c=getchar();
        return 1;
    }
};
int _(){
    int x=0,c=getchar();
    while(c<48)c=getchar();
    while(c>47)x=x*10+c-48,c=getchar();
    return x;
}
const u32 M=1<<22;
struct hmp{
    u64 hx[M];
    int hy[M],y0;
    int&operator[](u64 x){
        if(!x)return y0;
        int w=x&(M-1);
        while(hx[w]){
            if(hx[w]==x)return hy[w];
            w=(w+1237)&(M-1);
        }
        hx[w]=x;
        return hy[w];
    }
}H;
int ps[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73};
int ss[33];
u32 iv[100007],va=0;
u32 inv(u32 a){
    if(a<=100000)return iv[a];
    u32 v=1;
    for(u32 n=P-2;n;n>>=1,a=u64(a)*a%P)if(n&1)v=u64(v)*a%P;
    return v;
}
void dfs(int w,int d,double s,u64 h){
    if(w){
        u64 h1=h,s0=0;
        u32 c=1;
        for(int i=0;i<w;++i){
            int x=ss[i];
            c*=(x+1);
            s0+=H[h1+pp[x-1]-pp[x]];
            h1+=(_eq-1)*pp[x];
        }
        va+=H[h1]=(s0+c)%P*inv(c-1)%P;
        for(int i=w-1;i>=0;--i){
            int x=ss[i];
            u64 h2=h1+(1-_eq)*pp[x];
            H[h2]=(H[h1]+H[h2+pp[x-1]-pp[x]])%P;
            h1=h2;
        }
    }
    for(int i=1;i<=d;++i){
        s*=ps[w];
        if(s>1.01e24)return;
        ss[w]=i;
        dfs(w+1,i,s,h+pp[i]);
    }
}
int main(){
    num x;
    u32 m;
    pp[1]=1;
    iv[1]=1;
    for(int i=2;i<=100000;++i)iv[i]=u64(P-P/i)*iv[P%i]%P;
    for(u32 i=2;i<=100;++i)pp[i]=pp[i-1]*149;
    dfs(0,100,1,0);
    while(x.read()){
        u64 h=0;
        for(m=_();m;--m){
            u32 y=_(),t=0;
            while(x.div(y))++t;
            h+=pp[t];
        }
        printf("%d\n",H[h*_eq]);
    }
    return 0;
}