import numpy as np


def fourier_modes(w, modes):
    """Circular-kernel coefficients sum_j w[j] exp(-2 pi i k j/N)."""
    w = np.asarray(w)
    modes = np.asarray(modes)
    return np.fft.fft(w)[modes % len(w)]


def newton_root(what, a=1.0, tau_r=1.0, tau_d=0.5, iterations=12, init=None):
    """Solve tau_r*lambda = -1 + a*what*exp(-lambda*tau_d)."""
    lam = complex((-1.0 + a * what) / tau_r if init is None else init)
    for _ in range(iterations):
        e = np.exp(-lam * tau_d)
        f = tau_r * lam + 1.0 - a * what * e
        df = tau_r + a * what * tau_d * e
        if abs(df) < 1e-12:
            break
        lam -= f / df
    return lam


def mode_regularizer(w, modes=(-3, -2, -1, 1, 2, 3), a=1.0,
                     tau_r=1.0, tau_d=0.5, epsilon=0.02, eta=0.0,
                     target_velocity=0.0, T=1.0):
    """Detached-root spectral penalty for selected Fourier modes."""
    modes = tuple(modes)
    hats = fourier_modes(w, modes)
    roots = np.array([newton_root(h, a, tau_r, tau_d) for h in hats])
    growth = roots.real
    stable = np.logaddexp(0.0, growth + epsilon) ** 2
    velocity = np.array([T * z.imag / (2.0 * np.pi * k)
                         for z, k in zip(roots, modes)])
    total = np.sum(stable + eta * (velocity - target_velocity) ** 2)
    return float(total), roots, hats