// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package integrate

import (
	"math"

	"sourcedock.dev/petrbalvin/tensor/internal/base"
)

// Oscillatory quadrature: the integral of a smooth amplitude against a
// sine or cosine of a high frequency, the shape every spectral
// reduction produces and one a plain adaptive rule pays for double: it
// must resolve the carrier, not the amplitude, so the evaluation count
// grows with the frequency and the per-panel rules start aliasing.
//
// The scheme is Filon-type. The interval splits into equal panels, the
// amplitude f is interpolated on each panel by a polynomial through
// Gauss-Legendre nodes, and the product of that polynomial with the
// oscillatory kernel is carried out exactly through per-panel weights.
// The error therefore tracks the smoothness of f alone and falls like
// the panel width to the interpolation order, no matter how large the
// frequency grows, while the plain adaptive rule must spend roughly
// twenty evaluations per carrier wavelength to see it at all.

// FilonOptions tunes IntegrateFilon. Nodes ≤ 0 means 16, the
// polynomial degree of the amplitude interpolant per panel is Nodes−1.
// Panels ≤ 0 means automatic: the count that keeps each panel at most
// about Nodes half-wavelengths of the carrier, the range where the
// moment construction below is exact to the rounding floor.
type FilonOptions struct {
	Panels int
	Nodes  int
}

// filonAlphaCap bounds the forced-panel moment phase: a panel may
// carry at most this many half-wavelengths of the carrier before the
// auxiliary rule that builds the weights would have to grow without
// bound. The automatic panel count never reaches it.
const filonAlphaCap = 4096.0

// IntegrateFilon returns the two definite integrals
//
//	cosIntegral = ∫ f(x)·cos(kx) dx,  sinIntegral = ∫ f(x)·sin(kx) dx
//
// over [a, b], the real and imaginary parts of ∫ f(x)·e^{ikx} dx. A
// reversed interval integrates in the negative direction and k = 0
// degenerates to the plain integral of f with a zero sine part. The
// construction is exact whenever f is a polynomial of degree below
// Nodes, so on smooth amplitudes the answer sits at the rounding floor
// even for frequencies whose carrier a sampled rule cannot see.
//
// Errors: NaN or infinite bounds, an infinite frequency, a NaN
// frequency, Nodes outside [2, 32], a forced Panels whose panels would
// carry more than filonAlphaCap half-wavelengths of the carrier, a span
// that overflows the float64 range, a frequency whose span product
// leaves no representable panel count, and an f that fails or returns
// a non-finite value.
func IntegrateFilon(f func(x float64) (float64, error), a, b, k float64, opts FilonOptions) (float64, float64, error) {
	if opts.Nodes <= 0 {
		opts.Nodes = 16
	}
	if opts.Nodes < 2 || opts.Nodes > 32 {
		return 0, 0, base.Errf("IntegrateFilon: Nodes must be between 2 and 32, got %d", opts.Nodes)
	}
	if math.IsNaN(a) || math.IsNaN(b) || math.IsNaN(k) {
		return 0, 0, base.Errf("IntegrateFilon: bounds and frequency must not be NaN")
	}
	if math.IsInf(a, 0) || math.IsInf(b, 0) {
		return 0, 0, base.Errf("IntegrateFilon: bounds must be finite, got [%g, %g]", a, b)
	}
	if math.IsInf(k, 0) {
		return 0, 0, base.Errf("IntegrateFilon: the frequency must be finite, got %g", k)
	}
	sign := 1.0
	if b < a {
		a, b = b, a
		sign = -1
	}
	if a == b {
		return 0, 0, nil
	}
	// Two finite bounds can still sit so far apart that their span
	// overflows: the panel width would be infinite and the carrier's
	// phase at the panel centre 0·Inf or k·Inf, a quiet NaN pair.
	if span := b - a; math.IsInf(span, 0) {
		return 0, 0, base.Errf("IntegrateFilon: the span from %g to %g overflows, leaving no representable panel width", a, b)
	}
	if opts.Panels > 0 {
		if alpha := math.Abs(k) * (b - a) / (2 * float64(opts.Panels)); alpha > filonAlphaCap {
			return 0, 0, base.Errf("IntegrateFilon: %d panels leave %g half-wavelengths of the carrier per panel, above the %g the weights can be built within; raise Panels or leave them automatic",
				opts.Panels, alpha, filonAlphaCap)
		}
	}
	panels := opts.Panels
	if panels <= 0 {
		panels = 1
		if k != 0 {
			// A panel of h carries |k|h/2 half-wavelengths; the cap at
			// Nodes keeps the moment construction in its exact range
			// and the interpolation error far under the floor. The
			// estimate can also leave the int range while still
			// finite, and the conversion of such a ceiling is
			// implementation-dependent garbage: on saturation it asks
			// for an unending loop, elsewhere it wraps negative and
			// the empty loop reports a quiet zero. Refuse anything
			// the platform's int cannot represent.
			est := math.Abs(k) * (b - a) / (2 * float64(opts.Nodes))
			if est >= math.MaxInt {
				return 0, 0, base.Errf("IntegrateFilon: the frequency %g over the span %g leaves no representable panel count", k, b-a)
			}
			panels = int(math.Ceil(est))
		}
	}
	h := (b - a) / float64(panels)
	alpha := k * h / 2

	nodes, _, err := GaussLegendreNodes(opts.Nodes)
	if err != nil {
		return 0, 0, err
	}
	wCos, wSin, err := filonWeights(nodes, alpha)
	if err != nil {
		return 0, 0, err
	}

	// One sweep over the panels: sample the amplitude at the nodes,
	// contract it with the weights into the panel's two amplitudes C
	// and S, and rotate them into place by the carrier's phase at the
	// panel centre.
	var cosTotal, sinTotal float64
	for p := range panels {
		centre := a + (float64(p)+0.5)*h
		half := h / 2
		var c, s float64
		for i := range nodes {
			fx, ferr := f(centre + half*nodes[i])
			if ferr != nil {
				return 0, 0, base.Errf("IntegrateFilon: %w", ferr)
			}
			if math.IsNaN(fx) || math.IsInf(fx, 0) {
				return 0, 0, base.Errf("IntegrateFilon: the amplitude returned the non-finite value %g on panel %d", fx, p)
			}
			c += wCos[i] * fx
			s += wSin[i] * fx
		}
		phase := k * centre
		cosP, sinP := math.Cos(phase), math.Sin(phase)
		cosTotal += cosP*c - sinP*s
		sinTotal += sinP*c + cosP*s
	}
	return sign * cosTotal * h / 2, sign * sinTotal * h / 2, nil
}

// filonWeights returns, for the Gauss-Legendre nodes of [-1, 1], the
// Filon weights: the exact integrals of each Lagrange basis polynomial
// against cos(αy) and sin(αy). With these the panel integral of the
// interpolating polynomial times the carrier is one dot product per
// part, and every trace of the carrier's phase lives in the weights,
// built once, never per panel.
//
// The basis moments come from a composite 32-point Gauss-Legendre rule
// whose subinterval count follows α, so the auxiliary rule resolves
// the carrier the amplitude is multiplied by; the automatic panel cap
// keeps that cost at one subinterval and the rule at the rounding
// floor.
func filonWeights(nodes []float64, alpha float64) (wCos, wSin []float64, err error) {
	m := len(nodes)
	// Barycentric weights of the interpolation nodes.
	bw := make([]float64, m)
	for i := range m {
		p := 1.0
		for j := range m {
			if j != i {
				p *= nodes[i] - nodes[j]
			}
		}
		if p == 0 {
			return nil, nil, base.Errf("IntegrateFilon: repeated interpolation nodes")
		}
		bw[i] = 1 / p
	}
	// The auxiliary rule: 32-point Gauss-Legendre over enough equal
	// subintervals of [-1, 1] that each carries at most 16
	// half-wavelengths of e^{iαy}.
	subs := 1
	if a := math.Abs(alpha); a > 16 {
		subs = int(math.Ceil(a / 16))
	}
	auxNodes, auxWeights, err := GaussLegendreNodes(32)
	if err != nil {
		return nil, nil, err
	}
	wCos = make([]float64, m)
	wSin = make([]float64, m)
	span := 2.0 / float64(subs)
	for s := range subs {
		lo := -1 + float64(s)*span
		for t := range auxNodes {
			// The aux nodes live on [-1, 1]; map them into the
			// subinterval [lo, lo+span] with the half-span the affine
			// change of variables carries.
			y := lo + span*0.5*(auxNodes[t]+1)
			// Barycentric evaluation of every basis polynomial at y,
			// with the exact hit a node coincidence asks for.
			den := 0.0
			exact := -1
			for i := range m {
				d := y - nodes[i]
				if d == 0 {
					exact = i
					break
				}
				den += bw[i] / d
			}
			cy, sy := math.Cos(alpha*y), math.Sin(alpha*y)
			w := span * 0.5 * auxWeights[t]
			for i := range m {
				var li float64
				if exact >= 0 {
					if i == exact {
						li = 1
					}
				} else {
					li = bw[i] / (y - nodes[i]) / den
				}
				wCos[i] += w * li * cy
				wSin[i] += w * li * sy
			}
		}
	}
	return wCos, wSin, nil
}