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)