Read problems statements in Mandarin chinese, Russian and Vietnamese as well. 

Chef is tired of solving boring math problems by hand. One of these problems is summing up the products of elements from each k-subset of the set [n]. Here, a k-subset is a subset containing exactly k elements and [n] refers to the set which contains all integers between 1 and n (inclusive). More formally:

Let’s denote this number by f(n, k). Note that f(n, 0) = 1.

Since calculating f(n, k) is too boring, Chef wants to know how the divisibility of f(n, k) by a given prime p depends on k. Specifically, for a fixed n, you should compute the number of ways to select k (0 ≤ k ≤ n) so that f(n, k) isn't divisible by p. After a while of thinking, Chef realized there might be too many ways to do that. Therefore, you should compute this number modulo 10^{9}+7.

------ Input ------ 

The first line of the input contains a single integer T denoting the number of test cases. The description of T test cases follows.
The first and only line of each test case contains two space-separated integers n and p.

------ Output ------ 

For each test case, print a single line containing one number — the number of ways to select k, modulo 10^{9}+7.

------ Constraints ------ 

$1 ≤ T ≤ 4$
$1 ≤ n < 10^{501}$
$n does not contain leading zeroes$
$2 ≤ p ≤ 100,000$
$p is prime$

------ Subtasks ------ 

Subtask #1 (10 points): n ≤ 5,000

Subtask #2 (40 points): n ≤ 100,000

Subtask #3 (50 points): original constraints

----- Sample Input 1 ------ 
1
4 2
----- Sample Output 1 ------ 
2
----- explanation 1 ------ 
Example case 1: The values of f(4, k) are as follows:

f(4, 0) = 1
f(4, 1) = 1+2+3+4 = 10
f(4, 2) = 1·2+2·3+3·4+4·1+1·3+2·4 = 35
f(4, 3) = 1·2·3+2·3·4+1·3·4+1·2·4 = 50
f(4, 4) = 1·2·3·4 = 24

Only k = 0 and k = 2 give numbers indivisible by p = 2, so the answer is 2.

Difficulty: VERY_HARD
Tags: Combinatorics

Complete search

Divide and conquer

Mathematics

Polynomials and generating functions
Skill types: Complete search

Starter code:

Solutions:
import numpy as np
mod = int(1000000000.0 + 7)

def stir(p, st, en):
	if st >= en:
		return [1, en]
	m = (st + en) // 2
	ret = np.polymul(stir(p, st, m), stir(p, m + 1, en))
	return np.mod(ret, p)
for _ in range(int(input())):
	(n, p) = map(int, input().split())
	nst = 1 if n % p == p - 1 else np.sum(stir(p, 1, n % p) != 0)
	n = (n + 1) // p
	while n > 0:
		(nst, n) = (nst * (n % p + 1) % mod, n // p)
	print(nst)


import numpy as np

def poly(n, p):
	a = np.array([[1, i] for i in range(1, n + 1)])
	while n != 1:
		b = []
		for i in range(0, n // 2):
			b.append(np.polymul(a[2 * i], a[2 * i + 1]))
			for j in range(len(b[i])):
				b[i][j] = b[i][j] % p
		if len(a) % 2 == 1:
			b.append(a[n - 1])
		n = (n + 1) // 2
		a = np.array(b)
	return a

def hola(n, p):
	ho = poly(n, p)
	hol = 0
	for i in range(0, len(ho[0])):
		if ho[0][i] % p != 0:
			hol = hol + 1
	return hol
t = int(input())
while t > 0:
	(n, p) = [int(x) for x in input().split()]
	t = t - 1
	n = n + 1
	b = n % p
	if b > 1:
		ans = hola(b - 1, p)
	else:
		ans = 1
	a = 1
	mod = 1000000007
	while n > 0:
		n = n // p
		a = a * (n % p + 1) % mod
	print(a * ans % mod)


import numpy as np
mod = 10 ** 9 + 7

def fn(alist, p):
	if len(alist) > 1:
		mid = len(alist) // 2
		lefthalf = alist[:mid]
		righthalf = alist[mid:]
		p1 = fn(lefthalf, p)
		p2 = fn(righthalf, p)
		p3 = p1 * p2
		for x in range(p3.order + 1):
			p3[x] = p3[x] % p
		return p3
	elif len(alist) == 1:
		pp = np.poly1d(alist[0])
		return pp
	else:
		alist = [1]
		return np.poly1d(alist)
t = int(input())
while t > 0:
	(n, p) = map(int, input().split())
	ans = 1
	N = (n + 1) % p
	n = (n + 1) // p
	while n > 0:
		ans = ans * ((n % p + 1) % mod) % mod
		n = n // p
	l = []
	for i in range(N):
		temp = []
		temp.append(1)
		temp.append(i)
		l.append(temp)
	mul = fn(l, p)
	mul1 = 0
	for x in range(mul.order + 1):
		if mul[x] % p != 0:
			mul1 = mul1 + 1
	print(ans * mul1 % mod)
	t = t - 1


import math
mod = 10 ** 9 + 7
from numpy import polymul

def ro(x):
	while x & x - 1:
		x = (x | x >> 1) + 1
	return max(x, 1)

def bas(a, p):
	(r1, r2) = divmod(a, p)
	k = r2
	ad = 1
	while r1 > 0:
		(r1, r2) = divmod(r1, p)
		ad = ad * (r2 + 1) % mod
	return ([k], ad)

def mul(a, b):
	fin = polymul(a, b)
	fin = [i % p for i in fin]
	return fin

def repmul(k):
	while len(k) != 1:
		temp = []
		i = 0
		while i < len(k):
			if i + 1 < len(k):
				temp.append(mul(k[i], k[i + 1]))
				i += 2
			else:
				temp.append(k[i])
				i += 1
		k = temp[:]
	return k[0]
for _ in range(int(input())):
	(nn, p) = [int(i) for i in input().split()]
	mil = []
	ck = 1
	ans = 1
	(lk, ans) = bas(nn + 1, p)
	if lk[0] == 0 or lk[0] == 1:
		print(ans)
		continue
	rg = lk[0] - 1
	nn = rg
	k = [[1, (i + 1) % p] for i in range(nn)]
	k = repmul(k)
	c = 0
	for i in k:
		if i % p != 0:
			c += 1
	ans = ans * c % mod
	print(ans)


import numpy
import sys
osszes = []
try:
	while True:
		currs_ = input().split(' ')
		currs = []
		for i in currs_:
			if len(i) == 0 or not i[0].isdigit():
				continue
			else:
				currs.append(i)
		osszes += currs
except EOFError as e:
	a = 0
szamok = [int(i) for i in osszes]
mod = 1000000000.0 + 7
t = szamok[0]
T = 0
it = 1
while T < t:
	n = szamok[it]
	p = szamok[it + 1]
	it += 2
	n += 1
	pinary = []
	tmp = n
	while tmp > 0:
		pinary.append(tmp % p)
		tmp //= p

	def szorzo(L, R):
		if L > R:
			return [0]
		if L == R:
			return [1, L]
		mid = (L + R) // 2
		(res1, res2) = (szorzo(L, mid), szorzo(mid + 1, R))
		res = numpy.polymul(res1, res2)
		for i in range(0, len(res)):
			res[i] = int(res[i]) % p
		return res
	ans = 1
	if pinary[0] > 1:
		cnt = 0
		res = szorzo(1, pinary[0] - 1)
		for i in range(0, len(res)):
			if res[i] > 0:
				cnt += 1
		ans = cnt
	for i in range(1, len(pinary)):
		ans = ans * (pinary[i] + 1) % mod
	print(int(ans))
	T += 1


from itertools import combinations

def mult_polynom(polynom, number, prime):
	for elem_idx in range(len(polynom)):
		polynom[elem_idx] *= number
	polynom.append(0)
	for elem_idx in range(len(polynom) - 1, 0, -1):
		polynom[elem_idx] += polynom[elem_idx - 1] // number
		polynom[elem_idx] %= prime
	polynom[0] %= prime

def solve_tail(tail, prime):
	ans = 0
	curr_product = [1]
	for i in range(1, tail + 1):
		mult_polynom(curr_product, i, prime)
	for el in curr_product:
		if el != 0:
			ans += 1
	return ans

def solve(number=1, prime=2):
	mod = 10 ** 9 + 7
	(number, prime) = (int(x) for x in input().split())
	p_numeric = []
	num = number
	if number % prime == prime - 1:
		number //= prime
		number += 1
		p_numeric.append(0)
	while number > 0:
		p_numeric.append(number % prime)
		number //= prime
	ans = 1
	for digit in p_numeric[1:]:
		ans *= digit + 1
		ans %= mod
	if p_numeric[0] < 10003:
		ans *= solve_tail(p_numeric[0], prime)
	elif prime - p_numeric[0] > prime // 2:
		ans *= p_numeric[0]
	else:
		ans *= p_numeric[0] + 1
	ans %= mod
	print(ans)
	return ans

def naive_solve(number, prime):
	sums = []
	for amount in range(0, number + 1):
		curr_sum = 0
		for comb in combinations(range(1, number + 1), amount):
			curr_prod = 1
			for elem in comb:
				curr_prod *= elem
			curr_sum += curr_prod
		curr_sum %= prime
		sums.append(curr_sum)
	ans = 0
	for el in sums:
		if el != 0:
			ans += 1
	return ans

def test():
	primes = [2, 3, 5, 7, 11, 13, 17, 19]
	was = []
	for number in range(1, 21):
		for prime in primes:
			naiv_ans = naive_solve(number, prime)
			ans = solve(number, prime)
			if naiv_ans != ans:
				print('WA', number, prime)
				print(ans, ' != ', naiv_ans)
				was.append((number, prime, ans, naiv_ans))
			else:
				print('OK')
	print(*was, sep='\n')
reqs_nr = int(input())
for req in range(reqs_nr):
	solve()


import sys
import numpy as np

def getints():
	return map(int, sys.stdin.readline().split())

def fftconv(l0, l1, thr=512):
	n0 = len(l0)
	n1 = len(l1)
	if n0 < thr or n1 < thr:
		return np.mod(np.convolve(l0, l1), p)
	n = max(n0, n1)
	k = 0
	while 1 << k < 2 * n - 1:
		k += 1
	N = 1 << k
	h = p // 2
	l0 = np.append((l0 + h) % p - h, np.zeros(N - n0, dtype=np.float64))
	l1 = np.append((l1 + h) % p - h, np.zeros(N - n1, dtype=np.float64))
	t0 = np.fft.rfft(l0)
	t1 = np.fft.rfft(l1)
	aa = np.fft.irfft(t0 * t1)
	return np.floor(aa[:2 * n - 1] + 0.5).astype(int) % p

def nnz(N):

	def cr(a, b):
		d = b - a
		if d == 0:
			return np.array([1])
		if d == 1:
			return np.array([a, 1])
		if d == 2:
			return np.array([a * (a + 1) % p, 2 * a + 1, 1])
		if d == 3:
			return np.array([a * (a + 1) * (a + 2) % p, (3 * a * (a + 2) + 2) % p, 3 * a + 3, 1])
		m = (a + b) // 2
		return fftconv(cr(a, m), cr(m, b))
	l0 = cr(1, N + 1)
	if N < 10:
		l0 = np.mod(l0, p)
	s = 0
	for k in range(N + 1):
		s += l0[k] != 0
	return s

def basep(n, p):
	l = []
	while n > 0:
		l.append(n % p)
		n //= p
	return l

def prod(l, M=None):
	t = 1
	for x in l:
		t *= x
		if M != None:
			t = t % M
	return t
(T,) = getints()
M = 10 ** 9 + 7
for _ in range(T):
	(n, p) = getints()
	bb = basep(n, p)
	r = 0
	while r < len(bb) and bb[r] == p - 1:
		r += 1
	if r == 0:
		b = nnz(bb[0]) * prod([x + 1 for x in bb[1:]], M)
	else:
		bb.append(0)
		bb[r] += 1
		b = prod([x + 1 for x in bb[r:]], M)
	print(b % M)


import numpy as np

def get_poly(le, ri, p):
	if le == ri:
		return np.poly1d(np.array([le, 1]).astype(np.int64))
	mi = (le + ri) // 2
	p0 = get_poly(le, mi, p)
	p1 = get_poly(mi + 1, ri, p)
	pol = np.polymul(p0, p1).c
	pol = [x % p for x in pol]
	return np.poly1d(pol)

def get_nnz(n, p):
	if n == 0:
		return 1
	pol = get_poly(1, n, p)
	return sum([x != 0 for x in pol])

def main():
	mod = 1000000000 + 7
	n_case = int(input())
	for _ in range(n_case):
		(n, p) = map(int, input().split())
		r = n % p
		n //= p
		if r == p - 1:
			r = 0
			n = n + 1
		ans = int(get_nnz(r, p))
		while n != 0:
			r = n % p
			ans = ans * (r + 1) % mod
			n //= p
		print(ans)
main()