from math import gcd, isqrt
from collections import defaultdict, Counter

def liouville(m):
    n = m
    f = []
    for d in [2, 3, 5]:
        while n % d == 0:
            f.append(d)
            n //= d
    inc = (4, 2, 4, 2, 4, 6, 2, 6)
    i, d = 0, 7
    while d * d <= n:
        while n % d == 0:
            n //= d
            f.append(d)

        d += inc[i]
        i = (i + 1) % 8
    if n > 1:
        f.append(n)
    cnt = Counter(f)
    return (-1)**sum(cnt.values())

print([liouville(n) for n in range(1,11)])
def solve(n):
    r = 0
    for x in range(1, n + 1):
        for y in range(1, n + 1):
            g = gcd(x, y)
            fs = isqrt(g)
            if fs**2 == g:
                r += 1
    return r

def sum_ent(n):
    return sum(n // k for k in range(1, n + 1))
def nums(n):
    s = 0
    d = n
    while d > 0:
        k = n // d
        r = (n // k - n // (k + 1))
        s += r * k

        d = n // (k + 1)
        
    return s

def num_ineff(n):
    s = 0
    for d in range(1, n +  1):
        k = n // d
        s += k**2 * liouville(d)
    return s
    
def num_squares(n):
    s = 0
    d = n
 
    while d > 0:
        k = n // d
        s += k**2 * sum(liouville(x) for x in range(n//(k+1)+1, n//k + 1))
        d = n // (k + 1)


    return s

def count_square_gcd_pairs(n):
    # --- Linear sieve for Liouville function ---
    lam = [0] * (n + 1)
    lam[1] = 1
    primes = []
    is_composite = bytearray(n + 1)
    min_prime = [0] * (n + 1)   # smallest prime factor

    for i in range(2, n + 1):
        if not is_composite[i]:
            primes.append(i)
            lam[i] = -1
            min_prime[i] = i
        for p in primes:
            if i * p > n:
                break
            is_composite[i * p] = 1
            min_prime[i * p] = p
            if i % p == 0:
                lam[i * p] = -lam[i]   # one more factor of p
                break
            lam[i * p] = -lam[i]       # p is a new prime factor

    # Liouville prefix sums L[i] = sum_{d<=i} lambda(d)
    L = [0] * (n + 2)
    for i in range(1, n + 1):
        L[i] = L[i - 1] + lam[i]
    print(lam)
    print(L)
    # --- O(√n) block summation ---
    result = 0
    d = 1
    while d <= n:
        q = n // d
        next_d = n // q          # last index with floor(n/d) == q
        result += (L[next_d] - L[d - 1]) * q * q
        d = next_d + 1

    return result


n = 30
r = solve(n)
t = num_ineff(n)
s = num_squares(n)
print("s :", s, t, r, count_square_gcd_pairs(n))