import cmath

def phase_distribution(bits, theta):
    if type(bits) is not int or bits < 1 or not 0.0 <= theta < 1.0:
        raise ValueError("bits must be positive and phase must be in [0, 1)")
    size = 2 ** bits
    probabilities = []
    for outcome in range(size):
        amplitude = sum(cmath.exp(2j * cmath.pi * k * (theta - outcome / size))
                        for k in range(size)) / size
        probabilities.append(abs(amplitude) ** 2)
    return probabilities

def estimate(bits, theta):
    probabilities = phase_distribution(bits, theta)
    peak = max(range(len(probabilities)), key=probabilities.__getitem__)
    return {"bits": bits, "peak": peak, "phase": peak / (2 ** bits),
            "peak_probability": probabilities[peak], "mass": sum(probabilities)}

exact = estimate(3, 3 / 8)
coarse = estimate(3, 1 / 3)
fine = estimate(5, 1 / 3)
invalid_rejected = False
try:
    estimate(3, 1.0)
except ValueError:
    invalid_rejected = True
assert exact["peak"] == 3 and abs(exact["mass"] - 1.0) < 1e-12
assert abs(fine["phase"] - 1 / 3) < abs(coarse["phase"] - 1 / 3)
assert invalid_rejected and exact["peak_probability"] > 1.0 - 1e-12
print(f"PASS: 32 phase estimation exact={exact['peak']:03b} coarse_error={abs(coarse['phase']-1/3):.6f} fine_error={abs(fine['phase']-1/3):.6f}")