"""
FourierKunzVector — the normalised Kunz step function on S^1 with Fourier analysis.
"""
__all__ = ['FourierKunzVector']
import math
from ..numerical_semigroup import NumericalSemigroup
from ._vector import KunzVector
[docs]
class FourierKunzVector:
"""
The normalized Kunz function viewed as a step function on the circle S^1.
Starting from a KunzVector v = (w_1, ..., w_{m-1}) we build:
f : S^1 → [0, 1]
f(i/m) = w_i / max_j(w_j) for i = 1, ..., m-1
f(0) = 0 (convention: 0 residue class maps to 0)
The domain is discretised as {0, 1/m, 2/m, ..., (m-1)/m} ⊂ [0,1)
and f is extended to all of [0,1) as a right-continuous step function
(piecewise-constant on each interval [i/m, (i+1)/m)).
Parameters
----------
source : KunzVector | NumericalSemigroup
Either a KunzVector or a NumericalSemigroup (converted automatically).
"""
def __init__(self, source):
if isinstance(source, NumericalSemigroup):
source = KunzVector(source)
if not isinstance(source, KunzVector):
raise TypeError("source must be a KunzVector or NumericalSemigroup.")
self._kunz = source
m = source.multiplicity
if m == 1:
# The semigroup generated by 1 has an empty Kunz vector, so there
# is nothing to normalise. Keep only the residue-0 grid value.
raw_max = 0.0
self._grid: tuple[float, ...] = (0.0,)
else:
raw_max = max(source) # max over w_1 ... w_{m-1}
# normalised values at each grid point i/m, i = 0, 1, ..., m-1
# index 0 (residue 0) -> 0 by convention
self._grid: tuple[float, ...] = (0.0,) + tuple(
source[i] / raw_max for i in range(m - 1)
)
self._m = m
self._max_raw = raw_max
# --- properties ---
@property
def kunz_vector(self) -> KunzVector:
return self._kunz
@property
def multiplicity(self) -> int:
return self._m
@property
def grid_points(self) -> tuple[float, ...]:
"""The m normalised domain points {0, 1/m, ..., (m-1)/m}."""
return tuple(i / self._m for i in range(self._m))
@property
def grid_values(self) -> tuple[float, ...]:
"""Normalised function values at each grid point."""
return self._grid
# --- evaluation ---
def __call__(self, x: float) -> float:
"""
Evaluate f at x ∈ [0, 1).
f is the right-continuous step function that equals w_i / max(w)
on the half-open interval [i/m, (i+1)/m) for i = 0, 1, ..., m-1,
with the convention that f on [0, 1/m) equals 0 (residue 0 class).
Parameters
----------
x : float
A value in [0, 1). x = 1.0 is also accepted and mapped to 0
(periodicity).
Returns
-------
float
f(x) ∈ [0, 1].
"""
x = float(x) % 1.0 # enforce periodicity
m = self._m
# which bin does x fall into? bin i covers [i/m, (i+1)/m)
i = int(x * m)
if i >= m: # numerical edge case
i = m - 1
return self._grid[i]
[docs]
def evaluate(self, x: float) -> float:
"""Alias for __call__."""
return self(x)
# --- Fourier coefficients ---
[docs]
def fourier_coefficient(self, n: int) -> complex:
"""
Compute the n-th Fourier coefficient of f.
Because f is piecewise-constant on each half-open interval
``[i/m, (i+1)/m)``, the integral
.. math::
c_n = \\int_0^1 f(x)\\, e^{-2\\pi i n x}\\, dx
reduces **exactly** to a finite sum. Integrating ``f(i/m) * e^{-2πinx}``
over ``[i/m, (i+1)/m)`` gives ``(1/m) * f(i/m) * e^{-2πi n (i/m)}``, so:
.. math::
c_n = \\frac{1}{m} \\sum_{i=0}^{m-1} f\\!\\left(\\frac{i}{m}\\right)
e^{-2\\pi i n i/m}
This is the standard DFT of the grid values, scaled by ``1/m``.
Parameters
----------
n : int
Fourier mode index.
Returns
-------
complex
The n-th Fourier coefficient c_n.
"""
m = self._m
omega = 2 * math.pi / m
total = sum(
self._grid[i] * complex(math.cos(omega * n * i),
-math.sin(omega * n * i))
for i in range(m)
)
return total / m
[docs]
def fourier_coefficients(self, n_max: int) -> dict[int, complex]:
"""
Compute Fourier coefficients c_n for n = -n_max, ..., n_max.
Uses numpy's FFT when available (O(m log m)), otherwise falls back
to the per-coefficient DFT loop (O(m * n_max)).
Returns
-------
dict mapping int -> complex
"""
try:
import numpy as np
m = self._m
# np.fft.fft gives F[k] = sum_{j} grid[j] * e^{-2pi i j k / m}
# our convention is c_n = (1/m) * F[n]
full = np.fft.fft(self._grid) / m
return {n: complex(full[n % m]) for n in range(-n_max, n_max + 1)}
except ImportError:
return {n: self.fourier_coefficient(n) for n in range(-n_max, n_max + 1)}
[docs]
def partial_sum(self, x: float, n_max: int) -> float:
"""
Evaluate the Fourier partial sum
:math:`S_{n_{\\max}}(x) = \\sum_{|n| \\le n_{\\max}} c_n e^{2\\pi i n x}`.
Useful for visualising how well the Fourier series reconstructs f.
Parameters
----------
x : float
Point in [0, 1).
n_max : int
Number of modes on each side.
Returns
-------
float
Real part of the partial sum (imaginary part is ~0 for real f).
"""
coeffs = self.fourier_coefficients(n_max)
total = sum(
coeffs[n] * complex(math.cos(2 * math.pi * n * x),
math.sin(2 * math.pi * n * x))
for n in range(-n_max, n_max + 1)
)
return total.real
# --- L^p distance ---
[docs]
def distance(self, other: "FourierKunzVector", norm: str = "L2") -> float:
"""
Compute the distance between this FourierKunzVector and another.
Both functions are evaluated on a common grid of size lcm(m_self, m_other),
which is the finest grid on which both step functions are simultaneously
constant.
Parameters
----------
other : FourierKunzVector
norm : {"L1", "L2", "Linf"}
Which norm to use:
- ``"L1"`` — :math:`\\int_0^1 |f(x) - g(x)| \\, dx`
- ``"L2"`` — :math:`\\left(\\int_0^1 |f(x) - g(x)|^2 \\, dx\\right)^{1/2}` *(default)*
- ``"Linf"`` — :math:`\\sup_{x \\in [0,1)} |f(x) - g(x)|`
Returns
-------
float
The distance ≥ 0.
Raises
------
ValueError
If ``norm`` is not one of the supported values.
TypeError
If ``other`` is not a FourierKunzVector, or if ``norm`` is not a string.
"""
if not isinstance(other, FourierKunzVector):
raise TypeError("other must be a FourierKunzVector.")
if not isinstance(norm, str):
raise TypeError(f"norm must be a string ('L1', 'L2', or 'Linf'), got {type(norm).__name__}.")
norm = norm.upper()
if norm not in ("L1", "L2", "LINF"):
raise ValueError(f"norm must be 'L1', 'L2', or 'Linf', got {norm!r}.")
N = math.lcm(self._m, other._m)
if norm == "L1":
total = 0.0
for k in range(N):
total += abs(self(k / N) - other(k / N))
return total / N
elif norm == "L2":
total_sq = 0.0
for k in range(N):
d = abs(self(k / N) - other(k / N))
total_sq += d * d
return math.sqrt(total_sq / N)
else: # LINF
max_diff = 0.0
for k in range(N):
d = abs(self(k / N) - other(k / N))
if d > max_diff:
max_diff = d
return max_diff
# --- display ---
def __repr__(self) -> str:
pts = ", ".join(f"{x:.3f}" for x in self._grid)
return (f"FourierKunzVector(m={self._m}, "
f"grid_values=({pts}))")