You are given ${D}$ datasets where each dataset is in the form of two integers, $m$ and ${a}$, such that:
$n=\prod\limits_{i=1}^m p_i^{a+i},\text{where}~p_i~\text{is the}~i^{th}~\text{prime}.$

For each dataset, find and print the following on a new line: 

$\sum_{d|n}\sigma_0(d)$  

where $\sigma_0(x)$ is the count of divisors of ${x}$. As the answer can be quite large, print the result of this value modulo $(10^9+7)$.  

Input Format

The first line contains an integer, ${D}$, denoting the number of datasets. 

Each line ${i}$ of the ${D}$ subsequent lines contains two space-separated integers describing the respective values of $m_i$ and $a_i$ for dataset ${i}$.  

Constraints

$1\leq D\leq10^5$  
$1\leq m\leq10^5$  
$0\leq a\leq10^5$

Output Format

For each dataset, print a single integer denoting the result of the summation above modulo $(10^9+7)$ on a new line.

Sample Input
3
2 0
3 0
2 4

Sample Output
18
180
588

Explanation

For the first dataset where $m=2$ and ${a}=0$,
$\begin{aligned}&n=2^1\times3^2\\\Rightarrow&2\times9\\\Rightarrow&18$
$18}$ has the following divisors: $\{1,2,3,6,9,18\}$. Therefore, the result is:

$\begin{aligned}&\sigma_0(1)+\sigma_0(2)+\sigma_0(3)+\sigma_0(6)+\sigma_0(9)+\sigma_0(18)\\ \Rightarrow&1+2+2+4+3+6\\ \Rightarrow&18\end{aligned}$

Thus we print the value of $18\%(10^9+7)=18$ on a new line.

Difficulty: MEDIUM_HARD
Tags: Number theory

Starter code:

Solutions:
def gcd_ext(a, b):
	s = 0
	t = 1
	while b:
		(q, r) = divmod(a, b)
		(s, t) = (t, s - q * t)
		(a, b) = (b, r)
	return s

def cache_factorials(n, p):
	factorials = [1]
	for i in range(1, n + 1):
		factorials.append(factorials[-1] * i % p)
	inv_factorials = [gcd_ext(p, factorials[n])]
	for i in range(n, 0, -1):
		inv_factorials.append(inv_factorials[-1] * i % p)
	inv_factorials.reverse()
	return (factorials, inv_factorials)

def main():
	p = 10 ** 9 + 7
	(f, inv_f) = cache_factorials(200002, p)
	for i in range(int(input())):
		(m, a) = map(int, input().split())
		result = f[a + m + 1] * f[a + m + 2] * inv_f[a + 1] * inv_f[a + 2]
		result *= gcd_ext(p, pow(2, m, p))
		print(result % p)
main()


d = int(input())
mod = 10 ** 9 + 7
i2 = pow(2, mod - 2, mod)
_modfac = [1]

def modfac(x):
	while len(_modfac) <= x:
		_modfac.append(_modfac[-1] * len(_modfac) % mod)
	return _modfac[x]
for _ in range(d):
	(m, a) = map(int, input().split())
	a += 2
	x = m + a
	res = modfac(x) * modfac(x - 1) % mod * pow(modfac(a) * modfac(a - 1), mod - 2, mod) % mod * pow(i2, m, mod) % mod
	print(res)


D = int(input())
MOD = 1000000007

def modinv(a, m):
	(u1, u2, u3) = (1, 0, a)
	(v1, v2, v3) = (0, 1, m)
	while v3 != 0:
		q = u3 // v3
		(v1, v2, v3, u1, u2, u3) = (u1 - q * v1, u2 - q * v2, u3 - q * v3, v1, v2, v3)
	return u1 % m
factors = [1] * 200008
for I in range(2, 200008):
	factors[I] = I * factors[I - 1] % MOD
for _ in range(D):
	(M, A) = [int(X) for X in input().split(' ')]
	ans = factors[M + A + 2] * factors[M + A + 1] * modinv(factors[A + 2], MOD) * modinv(factors[A + 1], MOD) * modinv(pow(2, M, MOD), MOD)
	print(ans % MOD)


def solve(datasets):
	M = max(map(sum, datasets))
	negPowersOf2 = [1]
	facts = [1]
	factsInv = [1]
	for i in range(1, M + 3):
		negPowersOf2.append(negPowersOf2[-1] * 500000004 % 1000000007)
		facts.append(facts[-1] * i % 1000000007)
		factsInv.append(pow(facts[-1], 1000000005, 1000000007))
	A = []
	for (m, a) in datasets:
		A.append(facts[a + m + 1] * facts[a + m + 2] * factsInv[a + 1] * factsInv[a + 2] * negPowersOf2[m] % 1000000007)
	return A
d = int(input())
D = [tuple((int(x) for x in input().split())) for i in range(d)]
A = solve(D)
print(*A, sep='\n')


from functools import reduce
from itertools import accumulate
cache = list(accumulate((pow(x, 2, 1000000007) for x in range(1, 200002)), lambda x, y: x * y % 1000000007))
print(*((lambda m, a: reduce(lambda x, y: x * y % 1000000007, (cache[m + a], pow(cache[a + 1], 1000000005, 1000000007), a + 2, a + m + 2, pow(2, 1000000006 - m, 1000000007))))(*map(int, input().split())) for _ in range(int(input()))), sep='\n')


MODULUS = 10 ** 9 + 7
T = int(input())

def triangle(n):
	return n * (n + 1) // 2
dp = [1]

def divide(x, y, mod):
	if x % y == 0:
		return x // y
	x += mod * divide(-x % y, mod % y, y)
	return x // y % mod
for test in range(T):
	(m, a) = [int(i) for i in input().split()]
	prod = 1
	while len(dp) <= m + a:
		dp.append(dp[-1] * triangle(len(dp) + 1) % MODULUS)
	print(divide(dp[m + a], dp[a], MODULUS))


N = 200003
p = 1000000007
t = [n * (n + 1) // 2 % p for n in range(N)]
cumprod = [1] * N
for i in range(2, N):
	cumprod[i] = cumprod[i - 1] * t[i] % p
cumprod_inv = [0] * N
for i in range(1, N):
	cumprod_inv[i] = pow(cumprod[i], p - 2, p)
for _ in range(int(input())):
	(m, a) = map(int, input().split())
	print(cumprod[a + m + 1] * cumprod_inv[a + 1] % p)


mod = int(1000000000.0 + 7)
c = [0] * 200005
for i in range(1, 200005):
	c[i] = c[i - 1] + i
for i in range(2, 200005):
	c[i] = c[i - 1] * c[i] % mod
for _ in range(int(input())):
	(m, a) = map(int, input().split())
	print(c[m + a + 1] * pow(c[a + 1], mod - 2, mod) % mod)


from functools import reduce
from itertools import accumulate
Dmax = 100000
Mmax = 100000
MOD = 1000000007
mul_mod = lambda x, y: x * y % MOD
s = list(accumulate(range(1, Dmax + Mmax + 1)))
p = list(accumulate(s, mul_mod))

def inverser(a):
	(t, newt) = (0, 1)
	(r, newr) = (MOD, a)
	while newr:
		quotient = r // newr
		(t, newt) = (newt, t - quotient * newt)
		(r, newr) = (newr, r - quotient * newr)
	if t < 0:
		t += MOD
	return t
for _ in range(int(input())):
	(m, a) = map(int, input().split())
	print(p[a + m] * inverser(p[a]) % MOD)


import sys
modulo = 1000000007
fact_mod_limit = 200002
fact_mod_inv_limit = 100002
pow_mod_inv_limit = 100000

def mod_inv(a):
	b = modulo
	bb = b
	(s, ss) = (1, 0)
	while b != 0:
		q = a // b
		(a, b) = (b, a % b)
		(s, ss) = (ss, s - q * ss)
	return s % bb
fact_mod = [1]
for i in range(fact_mod_limit):
	x = fact_mod[i] * (i + 1) % modulo
	fact_mod.append(x)
fact_mod_inv = [mod_inv(1)]
for i in range(fact_mod_inv_limit):
	x = fact_mod_inv[i] * mod_inv(i + 1) % modulo
	fact_mod_inv.append(x)
pow_mod_inv = [mod_inv(1)]
mod_inv_2 = mod_inv(2)
for i in range(pow_mod_inv_limit):
	x = pow_mod_inv[i] * mod_inv_2 % modulo
	pow_mod_inv.append(x)

def solve(m, a):
	solution = 1
	solution = solution * fact_mod[a + m + 2] % modulo
	solution = solution * fact_mod_inv[a + 1] % modulo
	solution = solution * fact_mod[a + m + 1] % modulo
	solution = solution * fact_mod_inv[a + 2] % modulo
	solution = solution * pow_mod_inv[m] % modulo
	return solution
tc_len = int(sys.stdin.readline())
for tc in range(tc_len):
	(m, a) = tuple((int(x) for x in sys.stdin.readline().split()))
	print(solve(m, a))


import math

def extended_gcd(aa, bb):
	(lastremainder, remainder) = (abs(aa), abs(bb))
	(x, lastx, y, lasty) = (0, 1, 1, 0)
	while remainder:
		(lastremainder, (quotient, remainder)) = (remainder, divmod(lastremainder, remainder))
		(x, lastx) = (lastx - quotient * x, x)
		(y, lasty) = (lasty - quotient * y, y)
	return (lastremainder, lastx * (-1 if aa < 0 else 1), lasty * (-1 if bb < 0 else 1))

def modinv(a, m):
	(g, x, y) = extended_gcd(a, m)
	if g != 1:
		raise ValueError
	return x % m
M = 10 ** 9 + 7
d = int(input())
results = [1]
for i in range(1, 2 * 10 ** 5 + 1):
	mul = i * (i + 1) // 2
	results.append(results[-1] * mul % M)
for tc in range(d):
	(m, a) = map(int, input().split())
	total = 1
	b = results[a + 1]
	k = results[a + m + 1]
	result = k % M * (modinv(b, M) % M) % M
	print(result)


mod = 1000000007
fact = [1 for i in range(200006)]
for i in range(1, 200006):
	fact[i] = fact[i - 1] * i % mod
d = int(input())
for _ in range(d):
	(a, offset) = [int(i) for i in input().split()]
	temp1 = fact[offset + 2] * fact[offset + 2] % mod
	temp2 = fact[offset + 1 + a] * fact[offset + 1 + a] % mod
	temp1 = pow(temp1, mod - 2, mod)
	temp = temp1 * temp2 % mod
	num = temp * (offset + 2) * (offset + 2 + a) % mod
	den = pow(2, a, mod)
	den = pow(den, mod - 2, mod)
	ans = num * den % mod
	print(ans)


MOD = 1000000007

def Extended_Euclid(a, b):
	if b == 0:
		return (a, 1, 0)
	else:
		(r, x, y) = Extended_Euclid(b, a % b)
		return (r, y, x - a // b * y)
d = int(input())
n = 200001
answers = [1] * n
x = 3
for i in range(1, n):
	a = i + 1
	multipier = a * (a + 1) // 2 % MOD
	answers[i] = answers[i - 1] * multipier % MOD
for _ in range(d):
	(m, a) = [int(el) for el in input().split()]
	numerator = answers[m + a]
	(r, x, denominator) = Extended_Euclid(MOD, answers[a])
	denominator %= MOD
	ans = numerator * denominator % MOD
	print(ans)


mod = 10 ** 9 + 7
prod = [1] * (2 * 10 ** 5 + 1)
for i in range(1, len(prod)):
	prod[i] = (prod[i - 1] * (i + 1) * (i + 2) >> 1) % mod
D = int(input())
for d in range(D):
	(m, a) = map(int, input().split())
	s = prod[m + a] * pow(prod[a], mod - 2, mod) % mod
	print(s)


mod = 10 ** 9 + 7
inv2 = pow(2, mod - 2, mod)
facs = [1 for _ in range(2 * 10 ** 5 + 1)]
for i in range(1, 2 * 10 ** 5 + 1):
	facs[i] = facs[i - 1] * i * i % mod
for _ in range(int(input())):
	(m, a) = [int(x) for x in input().split()]
	res = facs[a + m + 1] * pow(facs[a + 2], mod - 2, mod) * (a + 2) * (a + m + 2) * pow(inv2, m, mod) % mod
	print(res)


def bitize(x):
	L = []
	while x > 0:
		L.append(x % 2)
		x //= 2
	if not L:
		L = [0]
	return L

def pow_mod(a, b, m):
	b_bits = bitize(b)
	a_powers = [0 for x in b_bits]
	a_powers[0] = a % m
	if b_bits[0]:
		total = a % m
	else:
		total = 1
	for i in range(1, len(b_bits)):
		a_powers[i] = a_powers[i - 1] * a_powers[i - 1] % m
		if b_bits[i]:
			total *= a_powers[i]
			total %= m
	return total % m

def precompute_factorials(L, modulus):
	factorials = [1, 2]
	factorial_inverses = [1, 500000004]
	for i in range(2, L + 1):
		factorials.append(factorials[i - 1] * i * (i - 1) % modulus)
	return factorials

def compute_divisor_count(m, a, modulus, factorials, factorial_inverses, two_m_inverses):
	if factorial_inverses[a + 2] != -1:
		x = factorials[m + a + 2] * factorial_inverses[a + 2] * two_m_inverses[m] % modulus
	else:
		factorial_inverses[a + 2] = pow_mod(factorials[a + 2], modulus - 2, modulus)
		x = factorials[m + a + 2] * factorial_inverses[a + 2] * two_m_inverses[m] % modulus
	return (x, factorial_inverses)
modulus = pow(10, 9) + 7
factorials = precompute_factorials(2 * pow(10, 5) + 1, modulus)
factorial_inverses = [-1 for _ in factorials]
two_m_inverses = [1 for _ in factorials]
two_inverse = pow_mod(2, modulus - 2, modulus)
two_m_inverses[0] = 1
two_m_inverses[1] = two_inverse
for i in range(2, 2 * pow(10, 5) + 1):
	two_m_inverses[i] = two_inverse * two_m_inverses[i - 1] % modulus
Q = int(input())
for _ in range(0, Q):
	(m, a) = [int(x) for x in input().split()]
	(x, factorial_inverses) = compute_divisor_count(m, a, modulus, factorials, factorial_inverses, two_m_inverses)
	print(x)


import sys
from time import time

def modinv(n, p):
	return pow(n, p - 2, p)

def solve(n, a, m, ff):
	u = pow(ff[a + n + 2], 2, m)
	v = pow(2, n, m) * (a + 2)
	v %= m
	v *= a + n + 2
	v %= m
	v *= pow(ff[a + 1], 2, m)
	v %= m
	rval = u * modinv(v, m) % m
	return rval
tm = time()
m = pow(10, 9) + 7
(ff, n) = ([1], pow(10, 5))
for i in range(1, 2 * n + 3):
	ff.append(ff[-1] * i % m)
t = int(input())
for _ in range(t):
	(n, a) = map(int, input().split())
	print(solve(n, a, m, ff))
print('time elapsed: %.5f' % (time() - tm,), file=sys.stderr)


import sys
M = 1000000007

def modinv(a, m):
	(lastremainder, remainder) = (a, m)
	(x, lastx) = (0, 1)
	while remainder > 1:
		(lastremainder, (quotient, remainder)) = (remainder, divmod(lastremainder, remainder))
		(x, lastx) = (lastx - quotient * x, x)
	return x % m
rem_fact = [1]

def fact(n):
	try:
		return rem_fact[n]
	except IndexError:
		z = len(rem_fact)
		while z <= n:
			rem_fact.append(rem_fact[-1] * z % M)
			z += 1
		return rem_fact[n]

def xv(m, a):
	return pow(fact(m + a + 1), 2, M) * modinv(pow(fact(a + 2), 2, M), M) * (m + a + 2) * (a + 2) * modinv(1 << m, M) % M

def solve(line):
	r = []
	for i in range(int(next(line))):
		(m, a) = map(int, next(line).split())
		r.append(xv(m, a))
	print('\n'.join(map(str, r)))
solve(sys.stdin)


mod = 1000000007
fact = [1 for i in range(200006)]
for i in range(1, 200006):
	fact[i] = fact[i - 1] * i % mod
d = int(input())
for _ in range(d):
	(a, offset) = [int(i) for i in input().split()]
	temp1 = fact[offset + 2] * fact[offset + 2] % mod
	temp2 = fact[offset + 1 + a] * fact[offset + 1 + a] % mod
	temp1 = pow(temp1, mod - 2, mod)
	temp = temp1 * temp2 % mod
	num = temp * (offset + 2) * (offset + 2 + a) % mod
	den = pow(2, a, mod)
	den = pow(den, mod - 2, mod)
	ans = num * den % mod
	print(ans)


import os, sys, bisect

def sorted_list_add(a, x):
	bisect.insort(a, x)

def sorted_list_max_less_equal(a, x):
	i = bisect.bisect_left(a, x)
	if i == len(a):
		return a[-1] if a else None
	else:
		return a[i]

def modular_power(a, n, m):
	a = a % m
	if a < 0:
		a += m
	prod = 1
	pwr = a
	while n != 0:
		if n % 2 != 0:
			prod = prod * pwr % m
		n >>= 1
		if n == 0:
			break
		pwr = pwr * pwr % m
	return prod

def modular_inverse(a, m):
	m0 = m
	y = 0
	x = 1
	if m == 1:
		return 0
	while a > 1:
		q = a // m
		t = m
		m = a % m
		a = t
		t = y
		y = x - q * y
		x = t
	if x < 0:
		x = x + m0
	return x
factorial_cache = {}
factorial_cache_sorted = {}

def hr_factorial(n, modulus):
	dbg = False
	cache = factorial_cache.get(modulus)
	if cache is None:
		factorial_cache[modulus] = {}
		cache = factorial_cache[modulus]
		factorial_cache_sorted[modulus] = []
	res = cache.get(n)
	if res is not None:
		if dbg:
			print('Using cached factorial {}'.format(n), file=sys.stderr)
		return res
	kstart = None
	kstart = sorted_list_max_less_equal(factorial_cache_sorted[modulus], n)
	if kstart is None:
		kstart = 1
		p = 1
	else:
		p = cache[kstart]
	if dbg:
		print('Computing factorial {} from {}'.format(n, kstart), file=sys.stderr)
	for k in range(kstart + 1, n + 1):
		p = p * k % modulus
	cache[n] = p
	sorted_list_add(factorial_cache_sorted[modulus], n)
	return p

def hr_subfactorial(k1, k2, m):
	return hr_factorial(k2, m) * modular_inverse(hr_factorial(k1 - 1, m), m) % m

def solve(m, a):
	MOD = 10 ** 9 + 7
	ONE_HALF = 5 * 10 ** 8 + 4
	p = hr_subfactorial(a + 3, a + m + 1, MOD)
	p = p * p * (a + 2) * (a + m + 2) * modular_power(ONE_HALF, m, MOD) % MOD
	return p

def solveInputs(inputs):
	MOD = 10 ** 9 + 7
	factorialArgs = {a + 3 - 1 for (m, a) in inputs} | {a + m + 1 for (m, a) in inputs}
	for x in sorted(factorialArgs):
		hr_factorial(x, MOD)
	return [solve(m, a) for (m, a) in inputs]

def main():
	fptr = open(os.environ['OUTPUT_PATH'], 'w') if 'OUTPUT_PATH' in os.environ else sys.stdout
	d = int(input())
	inputs = []
	for d_itr in range(d):
		ma = input().split()
		m = int(ma[0])
		a = int(ma[1])
		inputs.append((m, a))
	results = solveInputs(inputs)
	for res in results:
		fptr.write(str(res) + '\n')
	if fptr != sys.stdout:
		fptr.close()
main()


import os
from functools import reduce
import sys
import math as mt
from itertools import combinations
from itertools import compress

def primes(n):
	sieve = bytearray([True]) * (n // 2)
	for i in range(3, int(n ** 0.5) + 1, 2):
		if sieve[i // 2]:
			sieve[i * i // 2::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
	return [2, *compress(range(3, n, 2), sieve[1:])]
primelist = primes(100000)

def factorization(n):
	pf = []
	for p in primeslist:
		if p * p > n:
			break
		count = 0
		while not n % p:
			n //= p
			count += 1
		if count > 0:
			pf.append((p, count))
	if n > 1:
		pf.append((n, 1))
	return pf

def divisors(n):
	divs = [1]
	for (p, e) in factorization(n):
		divs += [x * p ** k for k in range(1, e + 1) for x in divs]
	return divs

def modInverse(a, m):
	return power(a, m - 2, m)

def power(x, y, p):
	res = 1
	x = x % p
	while y > 0:
		if y & 1 == 1:
			res = res * x % p
		y = y >> 1
		x = x * x % p
	return res

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

def modDivide(a, b, m):
	a = a % m
	inv = modInverse(b, m)
	if inv == -1:
		return
	else:
		return inv * a % m
D = int(input())
L = [0] * 200005
L[0] = 1
M = 10 ** 9 + 7
for i in range(1, 200005):
	L[i] = L[i - 1] % M * ((i + 1) * (i + 2) // 2) % M % M
for _ in range(D):
	(m, a) = map(int, input().strip().split(' '))
	print(modDivide(L[m + a], L[a], M))


mod = int(1000000000.0 + 7)
c = [0] * 200005
for i in range(1, 200005):
	c[i] = c[i - 1] + i
for i in range(2, 200005):
	c[i] = c[i - 1] * c[i] % mod
for _ in range(int(input())):
	(m, a) = map(int, input().split())
	print(c[m + a + 1] * pow(c[a + 1], mod - 2, mod) % mod)


def gcd_ext(a, b):
	s = 0
	t = 1
	while b:
		(q, r) = divmod(a, b)
		(s, t) = (t, s - q * t)
		(a, b) = (b, r)
	return s

def cache_factorials(n, p):
	factorials = [1]
	for i in range(1, n + 1):
		factorials.append(factorials[-1] * i % p)
	inv_factorials = [gcd_ext(p, factorials[n])]
	for i in range(n, 0, -1):
		inv_factorials.append(inv_factorials[-1] * i % p)
	inv_factorials.reverse()
	return (factorials, inv_factorials)

def main():
	p = 10 ** 9 + 7
	(f, inv_f) = cache_factorials(200002, p)
	for i in range(int(input())):
		(m, a) = map(int, input().split())
		result = f[a + m + 1] * f[a + m + 2] * inv_f[a + 1] * inv_f[a + 2]
		result *= gcd_ext(p, pow(2, m, p))
		print(result % p)
main()


def gcd_ext(a, b):
	s = 0
	t = 1
	while b:
		(q, r) = divmod(a, b)
		(s, t) = (t, s - q * t)
		(a, b) = (b, r)
	return s

def cache_factorials(n, p):
	factorials = [1]
	for i in range(1, n + 1):
		factorials.append(factorials[-1] * i % p)
	inv_factorials = [gcd_ext(p, factorials[n])]
	for i in range(n, 0, -1):
		inv_factorials.append(inv_factorials[-1] * i % p)
	inv_factorials.reverse()
	return (factorials, inv_factorials)

def main():
	p = 10 ** 9 + 7
	(f, inv_f) = cache_factorials(200002, p)
	for i in range(int(input())):
		(m, a) = map(int, input().split())
		result = f[a + m + 1] * f[a + m + 2] * inv_f[a + 1] * inv_f[a + 2]
		result *= gcd_ext(p, pow(2, m, p))
		print(result % p)
main()