luogu4714 「数学」约数个数和（组合数学+miller_rabin判素数）

luogu6月月赛E。
当时绝对是脑抽了
其实$O(wK)$的暴力递推还是可以写的呀…
令$n=\prod p_i^{q_i}$
我们发现可以分开讨论每一个$p_i^{q_i}$，最后乘起来就好了，且只跟次数q有关
设$f[k][q]$表示$p^q$的k次因数个数。
那么有递推式$f[k][q]=\sum\limits_{i=0}^qf[k-1][i]$
$f[0][q]=q+1$

考虑这个东西的组合意义，可以看做有q+1个空档，插k+1个板，板可重叠
于是有$f[k][q]=C_{q+1+k}^{k+1}=C_{q+1+k}^{q}$
于是最后答案就是$\prod C_{q_i+1+k}^{q_i}$
虽然K很大，但是$q_i<60$，我们可以直接乘出来，下面阶乘逆元。

所以现在唯一的问题就是大数分解质因数了。
我们可以筛出1e6以内的质数，然后用这些数试除了n之后，n只有四种情况：1,p,pq,p^2.
1和p^2可以判掉，因此我们只需要判断一下剩下的n是一个大质数还是两个不同质数的乘积即可。
Miller-Rabin即可。

我一定要吐槽一下我的快速乘！我没取mod又wa了！毁我青春，费我时间！
因为会有精度误差，所以必须取mod！

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define inf 0x3f3f3f3f
#define N 1000010
#define eps 1e-8
inline char gc(){
    static char buf[1<<16],*S,*T;
    if(S==T){T=(S=buf)+fread(buf,1,1<<16,stdin);if(T==S) return EOF;}
    return *S++;
}
inline ll read(){
    ll x=0,f=1;char ch=gc();
    while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=gc();}
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=gc();
    return x*f;
}
int prime[N>>3],tot=0,cnt[110],inv[110],ans=1;
bool notprime[N];
inline void init(){
    notprime[1]=1;
    for(int i=2;i<=N-10;++i){
        if(!notprime[i]) prime[++tot]=i;
        for(int j=1;prime[j]*i<=N-10;++j){
            notprime[prime[j]*i]=1;
            if(i%prime[j]==0) break;
        }
    }
}
inline ll mul(ll x,ll y,ll mod){
    ll res=x*y-(ll)((long double)x/mod*y+eps)*mod;res%=mod;return res<0?res+mod:res;
}
inline ll ksm(ll x,ll k,ll mod){
    ll res=1;for(;k;k>>=1,x=mul(x,x,mod)) if(k&1) res=mul(res,x,mod);return res;
}
inline bool check_p(ll a,ll x,int y,ll n){
    x=ksm(a,x,n);if(x==1||x==n-1) return 1;
    for(int i=1;i<y;++i){
        x=mul(x,x,n);
        if(x==1) return 0;
        if(x==n-1) return 1;
    }return 0;
}
inline bool miller_rabin(ll n){
    if(n==2) return 1;
    if(n%2==0) return 0;
    ll x=n-1;int y=0;
    while(x%2==0) x/=2,++y;
    for(int i=1;i<=10000;++i)
        if(!check_p(rand()%(n-1)+1,x,y,n)) return 0;return 1;
}
inline int C(ll n,int m,int mod){
    int res=1;
    for(int i=0;i<m;++i)
        res=mul(res,n-i,mod);
    res=(ll)res*inv[m]%mod;
    return res;
}
int main(){
//  freopen("a.in","r",stdin);
    ll n=read(),K=read();init();int m=0;
    for(int i=1;i<=tot;++i){
        if(n%prime[i]) continue;++m;
        while(n%prime[i]==0) n/=prime[i],++cnt[m];
    }if(n!=1){
        int nn=sqrt(n);if((ll)nn*nn==n) cnt[++m]=2;
        else{
            if(miller_rabin(n)) cnt[++m]=1;
            else cnt[++m]=1,cnt[++m]=1;
        }
    }int mod=998244353;inv[0]=inv[1]=1;
    for(int i=2;i<=100;++i) inv[i]=(ll)inv[mod%i]*(mod-mod/i)%mod;
    for(int i=1;i<=100;++i) inv[i]=(ll)inv[i-1]*inv[i]%mod;
    for(int i=1;i<=m;++i) ans=(ll)ans*C(K+cnt[i]+1,cnt[i],mod)%mod;
    printf("%d\n",ans);
    return 0;
}