import math

def is_prime(n: int) -> bool:
    """Deterministic Miller-Rabin for 64‑bit integers, fallback for larger."""
    if n < 2:
        return False
    # small primes
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
    for p in small_primes:
        if n % p == 0:
            return n == p
    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1
    # test these bases – enough for n < 2^64
    for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
        if a % n == 0:
            continue
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = (x * x) % n
            if x == n - 1:
                break
        else:
            return False
    return True

def factor_cct(c: int):
    """
    Factor a semiprime c = p * q (both primes) using the CCT‑ODE approach.
    
    The method solves s^2 - d^2 = 4c by searching over s = ceil(sqrt(4c)) upward.
    When s^2 - 4c is a perfect square, we obtain p = (s + d)//2, q = (s - d)//2.
    This is equivalent to the iterative probabilistic difference algorithm but
    deterministic. It is most efficient when the two primes are close.
    
    Returns:
        (p, q) such that p <= q and p * q == c.
        
    Raises:
        ValueError: if no pair found (c not semiprime or unbalanced primes too far).
    """
    if c % 2 == 0:
        # trivial if even
        p = 2
        q = c // 2
        if p * q == c and is_prime(p) and is_prime(q):
            return p, q
        raise ValueError(f"{c} is not a product of two large primes (even but not 2*prime)")

    # start from sqrt(4c)
    s = math.isqrt(4 * c) + 1
    # optional: small bound to avoid infinite loops – can be increased if needed
    max_iter = 10 ** 7   # safety cap, adjust based on expected |p-q|
    for _ in range(max_iter):
        diff_sq = s * s - 4 * c
        if diff_sq < 0:
            s += 1
            continue
        d = math.isqrt(diff_sq)
        if d * d == diff_sq:
            # candidate pair (a, b) = (s-d, s+d)
            a = s - d
            b = s + d
            if a % 2 == 0 and b % 2 == 0:
                p = a // 2
                q = b // 2
                if p > q:
                    p, q = q, p
                if p * q == c and is_prime(p) and is_prime(q):
                    return p, q
        s += 1

    raise ValueError(f"No factor pair found for {c} within allowed iterations. "
                     f"The primes might be too far apart for this method.")


if __name__ == "__main__":
    # Example from the markdown
    import numpy as np
    c1 = np.random.randint(10000,20000)
    p1, q1 = factor_cct(c1)
    print(f"c = {c1} → {p1} * {q1} = {p1 * q1}")

    # A larger semiprime with close primes (typical RSA‑like)
    c2 = 101 * 103   # = 10403, primes are close
    p2, q2 = factor_cct(c2)
    print(f"c = {c2} → {p2} * {q2} = {p2 * q2}")

    # Another larger example (balanced, but not tiny)
    c3 = 1000003 * 1000033   # difference = 30
    p3, q3 = factor_cct(c3)
    print(f"c = {c3} → {p3} * {q3} = {p3 * q3}")