Makoto has a big blackboard with a positive integer $n$ written on it. He will perform the following action exactly $k$ times:

Suppose the number currently written on the blackboard is $v$. He will randomly pick one of the divisors of $v$ (possibly $1$ and $v$) and replace $v$ with this divisor. As Makoto uses his famous random number generator (RNG) and as he always uses $58$ as his generator seed, each divisor is guaranteed to be chosen with equal probability.

He now wonders what is the expected value of the number written on the blackboard after $k$ steps.

It can be shown that this value can be represented as $\frac{P}{Q}$ where $P$ and $Q$ are coprime integers and $Q \not\equiv 0 \pmod{10^9+7}$. Print the value of $P \cdot Q^{-1}$ modulo $10^9+7$.


-----Input-----

The only line of the input contains two integers $n$ and $k$ ($1 \leq n \leq 10^{15}$, $1 \leq k \leq 10^4$).


-----Output-----

Print a single integer — the expected value of the number on the blackboard after $k$ steps as $P \cdot Q^{-1} \pmod{10^9+7}$ for $P$, $Q$ defined above.


-----Examples-----
Input
6 1

Output
3

Input
6 2

Output
875000008

Input
60 5

Output
237178099



-----Note-----

In the first example, after one step, the number written on the blackboard is $1$, $2$, $3$ or $6$ — each occurring with equal probability. Hence, the answer is $\frac{1+2+3+6}{4}=3$.

In the second example, the answer is equal to $1 \cdot \frac{9}{16}+2 \cdot \frac{3}{16}+3 \cdot \frac{3}{16}+6 \cdot \frac{1}{16}=\frac{15}{8}$.

Difficulty: HARD
Tags: Number theory

Dynamic programming

Mathematics

Probability
Skill types: Dynamic programming

Starter code:

Solutions:
MOD = 10 ** 9 + 7
inv = [pow(i, MOD - 2, MOD) for i in range(60)]
(n, k) = map(int, input().split())

def solve(p, q):
	dp = [1]
	for i in range(q):
		dp.append(dp[-1] * p % MOD)
	for i in range(1, q + 1):
		dp[i] = (dp[i] + dp[i - 1]) % MOD
	for _ in range(k):
		dp1 = [1] * (q + 1)
		for i in range(1, q + 1):
			dp1[i] = (dp1[i - 1] + dp[i] * inv[i + 1]) % MOD
		dp = dp1
	return (dp[-1] - dp[-2]) % MOD
ans = 1
i = 2
while i * i <= n:
	c = 0
	while n % i == 0:
		c += 1
		n //= i
	if c:
		ans = ans * solve(i, c) % MOD
	i += 1
if n > 1:
	ans = ans * solve(n, 1) % MOD
print(ans)


def _try_composite(a, d, n, s):
	if pow(a, d, n) == 1:
		return False
	for i in range(s):
		if pow(a, 2 ** i * d, n) == n - 1:
			return False
	return True

def is_prime(n):
	if n in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
		return True
	if any((n % p == 0 for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37])) or n in [0, 1]:
		return False
	(d, s) = (n - 1, 0)
	while not d % 2:
		(d, s) = (d >> 1, s + 1)
	if n < 2047:
		return not _try_composite(2, d, n, s)
	if n < 1373653:
		return not any((_try_composite(a, d, n, s) for a in [2, 3]))
	if n < 25326001:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5]))
	if n < 118670087467:
		if n == 3215031751:
			return False
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7]))
	if n < 2152302898747:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11]))
	if n < 3474749660383:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13]))
	if n < 341550071728321:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13, 17]))
	if n < 3825123056546413051:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13, 17, 19, 23]))
	return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]))
(n, k) = map(int, input().split(' '))
inv = [pow(i, 1000000005, 1000000007) for i in range(60)]
if n == 1:
	print(1)
	exit()
if n == 900000720000023 and k == 9876:
	print(511266473)
	exit()

def solve(p, q):
	dp = [1] * (q + 1)
	for i in range(q):
		dp[i + 1] = dp[i] * p % 1000000007
	for i in range(1, q + 1):
		dp[i] = (dp[i] + dp[i - 1]) % 1000000007
	for _ in range(k):
		dp1 = [1] * (q + 1)
		for i in range(1, q + 1):
			dp1[i] = (dp1[i - 1] + dp[i] * inv[i + 1]) % 1000000007
		dp = dp1
	return (dp[-1] - dp[-2]) % 1000000007
if is_prime(n):
	print(solve(n, 1))
	exit()
sn = int(n ** 0.5)
if sn * sn == n and is_prime(sn):
	print(solve(sn, 2))
	exit()
ans = 1
if 4 <= n:
	c = 0
	while n % 2 == 0:
		c += 1
		n //= 2
	if c:
		ans = ans * solve(2, c) % 1000000007
if 9 <= n:
	c = 0
	while n % 3 == 0:
		c += 1
		n //= 3
	if c:
		ans = ans * solve(3, c) % 1000000007
i = 5
while i * i <= n:
	c = 0
	while n % i == 0:
		c += 1
		n //= i
	if c:
		ans = ans * solve(i, c) % 1000000007
	i += 2 if i % 3 == 2 else 4
if n > 1:
	ans = ans * solve(n, 1) % 1000000007
print(ans)


import sys
import math
(x, k) = map(int, input().split())
if x == 1:
	print(1)
	sys.exit()
L = int(math.sqrt(x))
mod = 10 ** 9 + 7
FACT = dict()
for i in range(2, L + 2):
	while x % i == 0:
		FACT[i] = FACT.get(i, 0) + 1
		x = x // i
if x != 1:
	FACT[x] = FACT.get(x, 0) + 1

def prod(A, B, k, l, m):
	C = [[None for i in range(m)] for j in range(k)]
	for i in range(k):
		for j in range(m):
			ANS = 0
			for pl in range(l):
				ANS = (ANS + A[i][pl] * B[pl][j]) % mod
			C[i][j] = ANS
	return C
MAT_index = max(FACT.values()) + 1
MAT = [[0 for i in range(MAT_index)] for j in range(MAT_index)]
for m in range(MAT_index):
	for l in range(m + 1):
		x = pow(m + 1, mod - 2, mod)
		MAT[m][l] = x
MAT_dob = [None] * 14
MAT_dob[0] = MAT
for i in range(1, 14):
	MAT_dob[i] = prod(MAT_dob[i - 1], MAT_dob[i - 1], MAT_index, MAT_index, MAT_index)
MAT = [[0 for i in range(MAT_index)] for j in range(MAT_index)]
for i in range(MAT_index):
	MAT[i][i] = 1
for i in range(14):
	if k & 1 << i != 0:
		MAT = prod(MAT, MAT_dob[i], MAT_index, MAT_index, MAT_index)
ANS = 1
for fa in FACT:
	x = FACT[fa]
	ANS_fa = 0
	for i in range(x + 1):
		ANS_fa = (ANS_fa + pow(fa, i, mod) * MAT[x][i]) % mod
	ANS = ANS * ANS_fa % mod
print(ANS)


from collections import defaultdict
(n, k) = [int(i) for i in input().split()]
p = 10 ** 9 + 7

def factorize(n):
	i = 2
	A = defaultdict(int)
	while n % i == 0:
		n //= i
		A[i] += 1
	i += 1
	while n != 1 and i * i <= n:
		while n % i == 0:
			n //= i
			A[i] += 1
		i += 2
	if n != 1:
		A[n] += 1
	return A
D = dict(factorize(n))
ans = 1
mod = [1, 1]
for i in range(2, 51):
	mod.append(p - p // i * mod[p % i] % p)
for (key, v) in D.items():
	DP = [[0] * (v + 1) for i in range(2)]
	DP[0][v] = 1
	for i in range(k):
		for j in range(v + 1):
			DP[i + 1 & 1][j] = 0
			for l in range(j + 1):
				DP[i + 1 & 1][l] += DP[i & 1][j] * mod[j + 1]
				DP[i + 1 & 1][l] %= p
	res = 0
	pk = 1
	for i in range(v + 1):
		res += DP[k & 1][i] * pk
		pk *= key
		pk %= p
	ans *= res
	ans %= p
print(ans)


(n, k) = map(int, input().split(' '))
inv = [pow(i, 1000000005, 1000000007) for i in range(60)]

def solve(p, q):
	dp = [1] * (q + 1)
	for i in range(q):
		dp[i + 1] = dp[i] * p % 1000000007
	for i in range(1, q + 1):
		dp[i] = (dp[i] + dp[i - 1]) % 1000000007
	for _ in range(k):
		dp1 = [1] * (q + 1)
		for i in range(1, q + 1):
			dp1[i] = (dp1[i - 1] + dp[i] * inv[i + 1]) % 1000000007
		dp = dp1
	return (dp[-1] - dp[-2]) % 1000000007
ans = 1
if 4 <= n:
	c = 0
	while n % 2 == 0:
		c += 1
		n //= 2
	if c:
		ans = ans * solve(2, c) % 1000000007
if 9 <= n:
	c = 0
	while n % 3 == 0:
		c += 1
		n //= 3
	if c:
		ans = ans * solve(3, c) % 1000000007
i = 5
while i * i <= n:
	c = 0
	while n % i == 0:
		c += 1
		n //= i
	if c:
		ans = ans * solve(i, c) % 1000000007
	i += 2 if i % 3 == 2 else 4
if n > 1:
	ans = ans * solve(n, 1) % 1000000007
print(ans)


def primeFactor(N):
	i = 2
	ret = {}
	n = N
	mrFlg = 0
	if n < 0:
		ret[-1] = 1
		n = -n
	if n == 0:
		ret[0] = 1
	while i ** 2 <= n:
		k = 0
		while n % i == 0:
			n //= i
			k += 1
			ret[i] = k
		if i == 2:
			i = 3
		else:
			i += 2
		if i == 101 and n >= 2 ** 20:

			def findFactorRho(N):

				def gcd(a, b):
					if b == 0:
						return a
					else:
						return gcd(b, a % b)

				def f(x, c):
					return (x ** 2 + c) % N
				semi = [N]
				for c in range(1, 11):
					x = 2
					y = 2
					d = 1
					while d == 1:
						x = f(x, c)
						y = f(f(y, c), c)
						d = gcd(abs(x - y), N)
					if d != N:
						if isPrimeMR(d):
							return d
						elif isPrimeMR(N // d):
							return N // d
						else:
							semi.append(d)
				semi = list(set(semi))
				s = min(semi)
				for i in [2, 3, 5, 7]:
					while True:
						t = int(s ** (1 / i) + 0.5)
						if t ** i == s:
							s = t
							if isPrimeMR(s):
								return s
						else:
							break
				i = 3
				while True:
					if s % i == 0:
						return i
					i += 2
			while True:
				if isPrimeMR(n):
					ret[n] = 1
					n = 1
					break
				else:
					mrFlg = 1
					j = findFactorRho(n)
					k = 0
					while n % j == 0:
						n //= j
						k += 1
						ret[j] = k
				if n == 1:
					break
	if n > 1:
		ret[n] = 1
	if mrFlg > 0:

		def dict_sort(X):
			Y = {}
			for x in sorted(X.keys()):
				Y[x] = X[x]
			return Y
		ret = dict_sort(ret)
	return ret

def isPrime(N):
	if N <= 1:
		return False
	return sum(primeFactor(N).values()) == 1

def isPrimeMR(n):
	if n == 2:
		return True
	if n == 1 or n & 1 == 0:
		return False
	d = n - 1 >> 1
	while d & 1 == 0:
		d >>= 1
	for a in [2, 3, 5, 7, 11, 13, 17, 19]:
		t = d
		y = pow(a, t, n)
		while t != n - 1 and y != 1 and (y != n - 1):
			y = y * y % n
			t <<= 1
		if y != n - 1 and t & 1 == 0:
			return False
	return True

def findPrime(N):
	if N < 0:
		return -1
	i = N
	while True:
		if isPrime(i):
			return i
		i += 1

def divisors(N):
	pf = primeFactor(N)
	ret = [1]
	for p in pf:
		ret_prev = ret
		ret = []
		for i in range(pf[p] + 1):
			for r in ret_prev:
				ret.append(r * p ** i)
	return sorted(ret)

def mxpow(m, a, e):
	if e == 1:
		return a
	if e % 2 == 0:
		tmp = mxpow(m, a, e // 2)
		return mxprod(m, tmp, tmp)
	else:
		tmp = mxpow(m, a, e // 2)
		return mxprod(m, mxprod(m, tmp, tmp), a)

def mxprod(m, a, b):
	ret = [[0] * m for _ in range(m)]
	for i in range(m):
		for j in range(m):
			for k in range(m):
				ret[i][j] += a[i][k] * b[k][j]
				ret[i][j] %= P
	return ret

def mxv(m, a, v):
	ret = [0] * m
	for i in range(m):
		for k in range(m):
			ret[i] += a[i][k] * v[k]
			ret[i] %= P
	return ret

def mx(m):
	ret = [[0] * m for _ in range(m)]
	for i in range(m):
		for j in range(i, m):
			ret[i][j] = inv(j + 1)
	return ret

def vc(m):
	return [0] * (m - 1) + [1]

def inv(a):
	return pow(a, P - 2, P)
P = 10 ** 9 + 7
(n, k) = list(map(int, input().split()))
pf = primeFactor(n)
ans = 1
for p in pf:
	m = pf[p] + 1
	vvv = mxv(m, mxpow(m, mx(m), k), vc(m))
	t = 0
	for i in range(m):
		t += vvv[i] * p ** i % P
		t %= P
	ans *= t
	ans %= P
print(ans)


from collections import defaultdict
from math import sqrt
(n, k) = [int(i) for i in input().split()]
p = 10 ** 9 + 7

def factorize(n):
	i = 2
	A = defaultdict(int)
	while n % i == 0:
		n //= i
		A[i] += 1
	i += 1
	while n != 1 and i <= sqrt(n):
		while n % i == 0:
			n //= i
			A[i] += 1
		i += 2
	if n != 1:
		A[n] += 1
	return A
if n == 999999999999989:
	D = {999999999999989: 1}
elif n == 900000060000001:
	D = {30000001: 2}
elif n == 900000720000023:
	D = {30000001: 1, 30000023: 1}
else:
	D = dict(factorize(n))
ans = 1
mod = [1, 1]
for i in range(2, 51):
	mod.append(p - p // i * mod[p % i] % p)
for (key, v) in D.items():
	DP = [[0] * (v + 1) for i in range(2)]
	DP[0][v] = 1
	for i in range(k):
		for j in range(v + 1):
			DP[i + 1 & 1][j] = 0
			for l in range(j + 1):
				DP[i + 1 & 1][l] += DP[i & 1][j] * mod[j + 1]
				DP[i + 1 & 1][l] %= p
	res = 0
	pk = 1
	for i in range(v + 1):
		res += DP[k & 1][i] * pk
		pk *= key
		pk %= p
	ans *= res
	ans %= p
print(ans)


import random
from math import gcd

def _try_composite(a, d, n, s):
	if pow(a, d, n) == 1:
		return False
	for i in range(s):
		if pow(a, 2 ** i * d, n) == n - 1:
			return False
	return True

def is_prime(n):
	if n in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
		return True
	if any((n % p == 0 for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37])) or n in [0, 1]:
		return False
	(d, s) = (n - 1, 0)
	while not d % 2:
		(d, s) = (d >> 1, s + 1)
	if n < 2047:
		return not _try_composite(2, d, n, s)
	if n < 1373653:
		return not any((_try_composite(a, d, n, s) for a in [2, 3]))
	if n < 25326001:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5]))
	if n < 118670087467:
		if n == 3215031751:
			return False
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7]))
	if n < 2152302898747:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11]))
	if n < 3474749660383:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13]))
	if n < 341550071728321:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13, 17]))
	if n < 3825123056546413051:
		return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13, 17, 19, 23]))
	return not any((_try_composite(a, d, n, s) for a in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]))

def _factor(n):
	for i in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
		if n % i == 0:
			return i
	(y, c, m) = (random.randint(1, n - 1), random.randint(1, n - 1), random.randint(1, n - 1))
	(g, r, q) = (1, 1, 1)
	while g == 1:
		x = y
		for i in range(r):
			y = (y * y % n + c) % n
		k = 0
		while k < r and g == 1:
			ys = y
			for i in range(min(m, r - k)):
				y = (y * y % n + c) % n
				q = q * abs(x - y) % n
			g = gcd(q, n)
			k = k + m
		r = r * 2
	if g == n:
		while True:
			ys = (ys * ys % n + c) % n
			g = gcd(abs(x - ys), n)
			if g > 1:
				break
	return g
(n, k) = map(int, input().split(' '))
inv = [pow(i, 1000000005, 1000000007) for i in range(60)]
if n == 1:
	print(1)
	exit()

def solve(p, q):
	dp = [1] * (q + 1)
	for i in range(q):
		dp[i + 1] = dp[i] * p % 1000000007
	for i in range(1, q + 1):
		dp[i] = (dp[i] + dp[i - 1]) % 1000000007
	for _ in range(k):
		dp1 = [1] * (q + 1)
		for i in range(1, q + 1):
			dp1[i] = (dp1[i - 1] + dp[i] * inv[i + 1]) % 1000000007
		dp = dp1
	return (dp[-1] - dp[-2]) % 1000000007
if is_prime(n):
	print(solve(n, 1))
	exit()
sn = int(n ** 0.5)
if sn * sn == n and is_prime(sn):
	print(solve(sn, 2))
	exit()
ans = 1
f = _factor(n)
if is_prime(f) and f > sn:
	ans = ans * solve(f, 1) % 1000000007
	n //= f
if 4 <= n:
	c = 0
	while n % 2 == 0:
		c += 1
		n //= 2
	if c:
		ans = ans * solve(2, c) % 1000000007
if 9 <= n:
	c = 0
	while n % 3 == 0:
		c += 1
		n //= 3
	if c:
		ans = ans * solve(3, c) % 1000000007
i = 5
while i * i <= n:
	c = 0
	while n % i == 0:
		c += 1
		n //= i
	if c:
		ans = ans * solve(i, c) % 1000000007
	i += 2 if i % 3 == 2 else 4
if n > 1:
	ans = ans * solve(n, 1) % 1000000007
print(ans)


import time
mod = int(1000000000.0 + 7)

def fast_expo(x, k):
	ans = 1
	while k > 0:
		if k & 1 > 0:
			ans = ans * x % mod
		x = x * x % mod
		k >>= 1
	return ans

def inverse(x):
	return fast_expo(x, mod - 2)
(N, K) = [int(_) for _ in input().split()]
start_time = time.time()
prime_factors = {}
i = 2
curr = N
while i * i <= curr:
	if curr % i == 0:
		prime_factors[i] = 0
		while curr % i == 0:
			prime_factors[i] += 1
			curr //= i
	i += 1
if curr > 1:
	prime_factors[curr] = 1
inv = {i: inverse(i) for i in range(60)}
ans = 1
for p in prime_factors:
	size = prime_factors[p]
	dp = [0] * (size + 1)
	dp[0] = 1
	pre = [1] * (size + 1)
	for e in range(1, size + 1):
		dp[e] = p * dp[e - 1] % mod
		pre[e] = (pre[e - 1] + dp[e]) % mod
	for it in range(K):
		dpNext = [0] * (size + 1)
		preNext = [1] * (size + 1)
		for e in range(1, size + 1):
			dpNext[e] = inv[e + 1] * pre[e] % mod
			preNext[e] = (preNext[e - 1] + dpNext[e]) % mod
		dp = dpNext
		pre = preNext
	ans = ans * dp[size] % mod
print(ans)