積性函式 $f(p\ast q)=f(p$

#include <bits/stdc++.h>
#define
 
 mem(ar,num) memset(ar,num,sizeof(ar))
#define me(ar) memset(ar,0,sizeof(ar))
#define lowbit(x) (x&(-x))
#define Pb push_back
#define  FI first
#define  SE second
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define IOS ios::sync_with_stdio(false)
#define
 
 DEBUG cout<<endl<<"DEBUG"<<endl; 
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
const int    prime = 999983;
const int    INF = 0x7FFFFFFF;
const LL     INFF =0x7FFFFFFFFFFFFFFF;
const double pi = acos(-1.0);
const double inf = 1e18;
const double eps = 1e-6;
const LL     mod = 1e9 + 7;
LL qpow(LL a,LL b){LL s=1;while(b>0){if(b&1)s=s*a%mod;a=a*a%mod;b>>=1;}return s;}
LL gcd(LL a,LL b) {return b?gcd(b,a%b):a;}
int dr[2][4] = {1,-1,0,0,0,0,-1,1};
typedef pair<int,int> P;
const int maxn = 64;
LL dp[maxn];
LL sum[maxn];
LL n,k;
LL inv[maxn];
typedef long long ll;
typedef unsigned long long u64;
u64 pmod(u64 a,u64 b,u64 p) {
   u64 d=(u64)floor(a*(long double)b/p+0.5);
   ll ret=a*b-d*p;
   if (ret<0) ret+=p;
   return ret;
}
 
LL F(LL p,LL r){
	  dp[0] = 1;
      for(int i = 0;i <= r; ++i)
            dp[i] = qpow(p,i);
        sum[0] = 1;
        for(int i = 1;i <= r; ++i)
            sum[i] = (dp[i]+sum[i-1])%mod;
        for(int j = 1;j <= k; ++j){
         for(int i = 1;i <= r; ++i){
            dp[i] = pmod(inv[i+1],sum[i],mod);
            if(sum[i-1]+dp[i] >= mod)
            	sum[i] = sum[i-1]+dp[i]-mod;
            else
            	sum[i] = sum[i-1]+dp[i];
            // sum[i] = (sum[i-1]+dp[i])%mod;//
         }
        }
        return dp[r];
}
int main(void)
{
    for(int i = 1;i <= 64; ++i)
      inv[i] = qpow(i,mod-2);
    cin>>n>>k;
    if(n == 1)
        return 0*puts("1");
    //LL m = sqrt(n);
    LL r;
    LL ans = 1;
    for(LL i = 2;i*i <= n; ++i){
        if(n % i == 0){
            r = 0;
            while(n%i == 0)
                n /= i,r++;
           ans = ans*F(i,r)%mod;
        }
    }
    if(n != 1)
       ans = ans*F(n,1)%mod;
    
    cout<<ans<<endl;


   return 0;
}