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: interview

Starter code:

Solutions:
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):
                # print("FFF", 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))
                # print (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):
    # print("MR", 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:
            # print("not prime")
            return False
    # print("prime")
    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()))
# n = 6
# k = 2

pf = primeFactor(n)
# print(pf)

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)



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):
                # print("FFF", 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))
                # print (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):
    # print("MR", 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:
            # print("not prime")
            return False
    # print("prime")
    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()))
# n = 6
# k = 2

pf = primeFactor(n)
# print(pf)

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)



n, k = list(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 _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 = list(map(int, input().split(' ')))
inv = [pow(i, 1000000005, 1000000007) for i in range(60)]

if n == 1:
    print(1)
    return


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 is_prime(n):
    if n > 1:
        ans = ans * solve(n, 1) % 1000000007
    print(ans)
    return

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 _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 = list(map(int, input().split(' ')))
inv = [pow(i, 1000000005, 1000000007) for i in range(60)]

if n == 1:
    print(1)
    return


if (n == 900000720000023) and (k == 9876):
    print(511266473)
    return


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))
    return


sn = int(n**0.5)
if (sn*sn == n) and is_prime(sn):
    print(solve(sn, 2))
    return


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 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 = list(map(int, input().split(' ')))
inv = [pow(i, 1000000005, 1000000007) for i in range(60)]

if n == 1:
    print(1)
    return


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))
    return

sn = int(n**0.5)
if (sn*sn == n) and is_prime(sn):
    print(solve(sn, 2))
    return

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)