def memoize(function): ''' stores results of old function calls ''' from functools import wraps memo = {} @wraps(function) def f(*args, **kwargs): key = args, tuple(sorted(kwargs.items())) if key not in memo: memo[key] = function(*args, **kwargs) return memo[key] f.memo = memo return f # constants mod = 10**9 + 7 K = 7 # functions @memoize def fac(n): ''' factorial ''' return 1 if n == 0 else n * fac(n-1) % mod @memoize def inv(n): ''' modular inverse ''' return 1 if n <= 1 else (mod - mod/n) * inv(mod%n) % mod def Ch(n,r): ''' binomial coefficient ''' return 1 if r == 0 else Ch(n,r-1) * (n-r+1) * inv(r) % mod @memoize def C(k,n): ''' number of subsets of size k with sum <= n ''' return 0 if n < k else 1 if k == 0 else (C(k,n-k) + C(k-1,n-k)) % mod @memoize def T(k,n): ''' sum of all subsets of size k with sum <= n ''' return 0 if n < k else 0 if k == 0 else (k*C(k,n) + T(k,n-k) + T(k-1,n-k)) % mod def matrix_inv(A): ''' matrix inverse of A ''' n = len(A) B = [[i==j for j in xrange(n)] for i in xrange(n)] for i in xrange(n): m = i while m < n and A[m][i] == 0: m += 1 assert m < n A[m], A[i] = A[i], A[m] B[m], B[i] = B[i], B[m] sc = pow(A[i][i], mod-2, mod) for j in xrange(n): A[i][j] = A[i][j] * sc % mod B[i][j] = B[i][j] * sc % mod for j in xrange(n): if i == j: continue sc = A[j][i] for k in xrange(n): A[j][k] = (A[j][k] - sc * A[i][k]) % mod B[j][k] = (B[j][k] - sc * B[i][k]) % mod return B @memoize def vand_inv(n): ''' inverse of vandermonde matrix of size n*n with a(i) = i. ''' return matrix_inv([[pow(i,j,mod) for j in xrange(n)] for i in xrange(n)]) def get_coeffs(a): ''' coefficients of polynomial P with P(i) = a[i] ''' n = len(a) B = vand_inv(n) return [sum(B[i][j] * a[j] for j in xrange(n)) % mod for i in xrange(n)] # compute polynomial coefficients coeffc = [None]*K coefft = [None]*K for k in xrange(1,K): fk = fac(k) coeffc[k] = [None]*fk coefft[k] = [None]*fk for b in xrange(fk): Vc = [C(k,a*fk+b) for a in xrange(k+1)] Vt = [T(k,a*fk+b) for a in xrange(k+2)] coeffc[k][b] = get_coeffs(Vc) coefft[k][b] = get_coeffs(Vt) def p_eval(P,x): ''' evaluate polynomial P at x ''' # horner's rule r = 0 for c in reversed(P): r = (x * r + c) % mod return r def CT(k,n): ''' Computes C_k(n) and T_k(n) ''' fk = fac(k) a, b = divmod(n, fk) Pc = coeffc[k][b] Pt = coefft[k][b] return p_eval(Pc, a), p_eval(Pt, a) def solve(k,n): c, t = CT(k-1,n) return (Ch(n,k) - (n+1)*c + t) % mod for cas in xrange(input()): n, k = map(int, raw_input().strip().split()) print solve(k, n)