Square-free division (hard version)
This is the hard version of the problem. The only difference is that in this version $$$0 \leq k \leq 20$$$.There is an array $$$a_1, a_2, \ldots, a_n$$$ of $$$n$$$ positive integers. You should divide it into a minimal number of continuous segments, such that in each segment there are no two numbers (on different positions), whose product is a perfect square.Moreover, it is allowed to do at most $$$k$$$ such operations before the division: choose a number in the array and change its value to any positive integer.What is the minimum number of continuous segments you should use if you will make changes optimally?
In the first test case it is possible to change the array this way: $$$[\underline{3}, 6, 2, 4, \underline{5}]$$$ (changed elements are underlined). After that the array does not need to be divided, so the answer is $$$1$$$.In the second test case it is possible to change the array this way: $$$[6, 2, \underline{3}, 8, 9, \underline{5}, 3, 6, \underline{10}, \underline{11}, 7]$$$. After that such division is optimal:   $$$[6, 2, 3]$$$  $$$[8, 9, 5, 3, 6, 10, 11, 7]$$$ 
The first line contains a single integer $$$t$$$ $$$(1 \le t \le 1000)$$$  — the number of test cases.The first line of each test case contains two integers $$$n$$$, $$$k$$$ ($$$1 \le n \le 2 \cdot 10^5$$$, $$$0 \leq k \leq 20$$$).The second line of each test case contains $$$n$$$ integers $$$a_1, a_2, \ldots, a_n$$$ ($$$1 \le a_i \le 10^7$$$).It's guaranteed that the sum of $$$n$$$ over all test cases does not exceed $$$2 \cdot 10^5$$$.
For each test case print a single integer  — the answer to the problem.
import sys, os, io
input = io.BytesIO(os.read(0, os.fstat(0).st_size)).readline

def the_sieve_of_eratosthenes(n):
    s = [1] * (n + 1)
    x = []
    for i in range(2, n + 1):
        if s[i]:
            x.append(i)
            for j in range(i, n + 1, i):
                s[j] = 0
    return x

def f(u, v):
    return u * (k + 1) + v
t = int(input())
p = the_sieve_of_eratosthenes(pow(10, 4))
l = pow(10, 7) + 5
x = [0] * (l + 1)
for i in range(l + 1):
    x[i] = i
for i in p:
    j = i * i
    for k in range(1, l // j + 1):
        x[j * k] = x[k]
inf = pow(10, 9) + 1
ans = []
for _ in range(t):
    (n, k) = map(int, input().split())
    a = list(map(int, input().split()))
    (d, y) = (dict(), [])
    for i in range(n):
        j = x[a[i]]
        if j in d:
            y.append(i)
        else:
            d[j] = len(d)
        a[i] = d[j]
    m = len(d)
    a.append(m)
    y.append(n)
    (nx, la) = ([-1] * (n + 1), [-1] * (m + 1))
    ok = [0] * (n + 1)
    for i in range(n + 1):
        if la[a[i]] ^ -1:
            nx[la[a[i]]] = i
        else:
            ok[i] = 1
        la[a[i]] = i
    ok[n] = 0
    dp = [inf] * ((n + 1) * (k + 1))
    dp[0] = 0
    z = []
    k0 = min(len(y), k + 1)
    for i in range(k0):
        ok[y[i]] = 1
        z.append(y[i])
    w = 0
    for i in range(n):
        for j in range(k + 1):
            if dp[f(i, j)] == inf:
                continue
            c = dp[f(i, j)] + 1
            v = j
            for u in z:
                dp[f(u, v)] = min(dp[f(u, v)], c)
                v += 1
                if v > k:
                    break
        if nx[i] ^ -1:
            ok[nx[i]] = 1
        for j in range(len(z)):
            if nx[i] <= z[j]:
                if nx[i] == z[j]:
                    z.pop(j)
                break
        while w <= n and ok[w]:
            w += 1
        if w <= n and len(z) ^ k0:
            z.append(w)
            w += 1
    ans0 = inf
    for i in range(f(n, 0), f(n + 1, 0)):
        ans0 = min(ans0, dp[i])
    ans.append(ans0)
sys.stdout.write('\n'.join(map(str, ans)))
def nr():
    return nrs()[0]

def nrs():
    return [int(i) for i in input().split()]
n = 10 ** 7
squares = [i * i for i in range(1, 3162)]
p = list(range(n + 1))
for i in range(1, n + 1):
    if p[i] == i:
        for sq in squares:
            if i * sq > n:
                break
            p[i * sq] = i
for _ in range(nr()):
    (N, K) = nrs()
    A = [p[a] for a in nrs()]
    dp = [N] * (K + 1)
    dp[0] = 1
    used = [{}] * (K + 1)
    for a in A:
        for j in range(K, -1, -1):
            if dp[j] == N:
                continue
            if a in used[j]:
                if j < K and dp[j + 1] > dp[j]:
                    dp[j + 1] = dp[j]
                    used[j + 1] = used[j]
                dp[j] += 1
                used[j] = {}
            used[j][a] = 1
    print(min(dp))
def nr():
    return nrs()[0]

def nrs():
    return [int(i) for i in input().split()]
n = 10 ** 7
squares = [i * i for i in range(1, 3162)]
p = list(range(n + 1))
for i in range(1, n + 1):
    if p[i] == i:
        for sq in squares:
            if i * sq > n:
                break
            p[i * sq] = i
for _ in range(nr()):
    (N, K) = nrs()
    A = [p[a] for a in nrs()]
    dp = [N] * (K + 1)
    dp[0] = 1
    used = [{}] * (K + 1)
    for a in A:
        for j in range(K, -1, -1):
            if dp[j] == N:
                continue
            if a in used[j]:
                if j < K and dp[j + 1] > dp[j]:
                    dp[j + 1] = dp[j]
                    used[j + 1] = used[j]
                dp[j] += 1
                used[j] = {}
            used[j][a] = 1
    print(min(dp))