Vlad and the Nice Paths (hard version)
This is hard version of the problem, it differs from the easy one only by constraints on $$$n$$$ and $$$k$$$.Vlad found a row of $$$n$$$ tiles and the integer $$$k$$$. The tiles are indexed from left to right and the $$$i$$$-th tile has the color $$$c_i$$$. After a little thought, he decided what to do with it.You can start from any tile and jump to any number of tiles right, forming the path $$$p$$$. Let's call the path $$$p$$$ of length $$$m$$$ nice if:  $$$p$$$ can be divided into blocks of length exactly $$$k$$$, that is, $$$m$$$ is divisible by $$$k$$$;  $$$c_{p_1} = c_{p_2} = \ldots = c_{p_k}$$$;  $$$c_{p_{k+1}} = c_{p_{k+2}} = \ldots = c_{p_{2k}}$$$;  $$$\ldots$$$  $$$c_{p_{m-k+1}} = c_{p_{m-k+2}} = \ldots = c_{p_{m}}$$$; Your task is to find the number of nice paths of maximum length. Since this number may be too large, print it modulo $$$10^9 + 7$$$.
In the first sample, it is impossible to make a nice path with a length greater than $$$0$$$.In the second sample, we are interested in the following paths:  $$$1 \rightarrow 3 \rightarrow 4 \rightarrow 5$$$  $$$2 \rightarrow 4 \rightarrow 5 \rightarrow 7$$$  $$$1 \rightarrow 3 \rightarrow 5 \rightarrow 7$$$  $$$1 \rightarrow 3 \rightarrow 4 \rightarrow 7$$$ In the third example, any path of length $$$8$$$ is nice.
The first line of each test contains the integer $$$t$$$ ($$$1 \le t \le 10^4$$$) — the number of test cases in the test.The first line of each test case contains two integers $$$n$$$ and $$$k$$$ ($$$1 \le k \le n \le 5000$$$) — the number of tiles in a row and the length of the block.The second line of each test case contains $$$n$$$ integers $$$c_1, c_2, c_3, \dots, c_n$$$ ($$$1 \le c_i \le n$$$) — tile colors.It is guaranteed that the sum of $$$n^2$$$ over all test cases does not exceed $$$25 \cdot 10^6$$$.
Print $$$t$$$ numbers, each of which is the answer to the corresponding test case — the number of nice paths of maximum length modulo $$$10^9 + 7$$$.
def make_nCr_mod(max_n=2 * 10 ** 5, mod=10 ** 9 + 7):
    max_n = min(max_n, mod - 1)
    (fact, inv_fact) = ([0] * (max_n + 1), [0] * (max_n + 1))
    fact[0] = 1
    for i in range(max_n):
        fact[i + 1] = fact[i] * (i + 1) % mod
    inv_fact[-1] = pow(fact[-1], mod - 2, mod)
    for i in reversed(range(max_n)):
        inv_fact[i] = inv_fact[i + 1] * (i + 1) % mod

    def nCr_mod(n, r):
        res = 1
        while n or r:
            (a, b) = (n % mod, r % mod)
            if a < b:
                return 0
            res = res * fact[a] % mod * inv_fact[b] % mod * inv_fact[a - b] % mod
            n //= mod
            r //= mod
        return res
    return nCr_mod
nCr_mod = make_nCr_mod(max_n=5000)

def main():
    mod = int(1000000000.0 + 7)
    t = int(input())
    allans = []
    for _ in range(t):
        (n, k) = readIntArr()
        c = readIntArr()
        cnts = [0] * n
        max_len = [0] * n
        for i in range(n):
            ccnt = 0
            for j in range(i, -1, -1):
                if c[j] == c[i]:
                    ccnt += 1
                if ccnt == k:
                    j -= 1
                    break
            if ccnt < k:
                cnts[i] = 1
                max_len[i] = 0
                continue
            best_len = 0
            while j >= 0:
                best_len = max(best_len, max_len[j])
                j -= 1
            max_len[i] = best_len + 1
            if best_len == 0:
                ccnt = 0
                for j in range(i):
                    if c[j] == c[i]:
                        ccnt += 1
                cnts[i] = nCr_mod(ccnt, k - 1)
                continue
            assert best_len > 0
            ccnt = 0
            for j in range(i - 1, -1, -1):
                if ccnt >= k - 1 and max_len[j] == best_len:
                    cnts[i] += cnts[j] * nCr_mod(ccnt, k - 1)
                    cnts[i] %= mod
                if c[j] == c[i]:
                    ccnt += 1
        l = 0
        ans = 1
        for i in range(n):
            if max_len[i] > l:
                l = max_len[i]
                ans = cnts[i]
            elif max_len[i] == l:
                ans += cnts[i]
                if ans >= mod:
                    ans -= mod
        if l == 0:
            ans = 1
        allans.append(ans)
    multiLineArrayPrint(allans)
    return
import sys
input = sys.stdin.buffer.readline

def oneLineArrayPrint(arr):
    print(' '.join([str(x) for x in arr]))

def multiLineArrayPrint(arr):
    print('\n'.join([str(x) for x in arr]))

def multiLineArrayOfArraysPrint(arr):
    print('\n'.join([' '.join([str(x) for x in y]) for y in arr]))

def readIntArr():
    return [int(x) for x in input().split()]

def makeArr(defaultValFactory, dimensionArr):
    dv = defaultValFactory
    da = dimensionArr
    if len(da) == 1:
        return [dv() for _ in range(da[0])]
    else:
        return [makeArr(dv, da[1:]) for _ in range(da[0])]

def queryInteractive(a, b):
    print('? {} {}'.format(a, b))
    sys.stdout.flush()
    return int(input())

def answerInteractive(ans):
    print('! {}'.format(ans))
    sys.stdout.flush()
inf = float('inf')
from math import gcd, floor, ceil
import math
for _abc in range(1):
    main()
import collections
import functools
import itertools
import sys
RD = sys.stdin.readline
INT = lambda : int(RD())
INTS = lambda : list(map(int, RD().split()))
mod = 10 ** 9 + 7
_c = [[1], [1, 1]]
for k in range(5000):
    n = len(_c[-1])
    q = [1] * (n + 1)
    for i in range(1, n):
        q[i] = (_c[-1][i] + _c[-1][i - 1]) % mod
    _c.append(q)

def C(n, m):
    assert n >= 0 and 0 <= m <= n
    return _c[n][m]
T = INT()
for _t in range(T):
    (n, k) = INTS()
    c = INTS()
    assert len(c) == n
    maxg = [0] * n
    d = collections.defaultdict(int)
    for i in range(n):
        if i > 0:
            maxg[i] = maxg[i - 1]
        d[c[i]] += 1
        if d[c[i]] >= k:
            maxg[i] += 1
            d.clear()
    if k == 1 or maxg[-1] == 0:
        print(1)
        continue
    di = collections.defaultdict(collections.deque)
    rete = [0] * n
    rett = [0] * n
    for i in range(n):
        q = di[c[i]]
        q.append(i)
        while len(q) >= k:
            b = q[0]
            if (maxg[b - 1] if b > 0 else 0) + 1 < maxg[i]:
                q.popleft()
            else:
                break
        lq = len(q)
        if lq >= k:
            for j in range(lq):
                if lq - j < k:
                    break
                _a = rett[q[j] - 1] if q[j] > 0 and rett[q[j] - 1] > 0 else 1
                _b = C(lq - j - 2, k - 2)
                rete[i] = (rete[i] + _a * _b % mod) % mod
                if rete[i] == 0:
                    rete[i] = mod
        rett[i] = (rett[i - 1] if i > 0 and maxg[i - 1] == maxg[i] else 0) + rete[i]
        if rett[i] > 0 and rett[i] % mod == 0:
            rett[i] = mod
        else:
            rett[i] %= mod
    print(rett[-1] % mod)
from sys import stdin
input = lambda : stdin.readline().strip()
M = 10 ** 9 + 7
cnk = [[0] * (5000 + 1) for i in range(5000 + 1)]

def solve():
    (n, k) = map(int, input().split())
    c = [-1] + [int(x) for x in input().split()]
    if k == 1:
        print(1)
        return
    dp = [[0, 0] for i in range(n + 1)]
    dp[0][0] = 1
    for i in range(1, n + 1):
        sz = 1
        for s in range(i - 1, -1, -1):
            if c[s] == c[i]:
                sz += 1
                if sz == k:
                    dp[i][1] = dp[s - 1][1] + 1
                if sz >= k:
                    if dp[s - 1][1] < dp[i][1] - 1:
                        break
                    dp[i][0] += dp[s - 1][0] * cnk[sz - 2][k - 2] % M
                    dp[i][0] %= M
        if dp[i][1] < dp[i - 1][1]:
            dp[i] = [0, dp[i - 1][1]]
        if dp[i][1] == dp[i - 1][1]:
            dp[i][0] += dp[i - 1][0]
            dp[i][0] %= M
    print(dp[n][0])
for i in range(5000 + 1):
    cnk[i][0] = 1
for i in range(1, 5000 + 1):
    for j in range(1, i + 1):
        cnk[i][j] = (cnk[i - 1][j] + cnk[i - 1][j - 1]) % M
t = int(input())
for _ in range(t):
    solve()