sparse_ir_basis/
polyfourier.rs

1//! Piecewise Legendre polynomial Fourier transform implementation for SparseIR
2//!
3//! This module provides Fourier transform functionality for piecewise Legendre
4//! polynomials, enabling evaluation in Matsubara frequency domain.
5
6use num_complex::Complex64;
7use std::f64::consts::PI;
8
9use crate::error::Error;
10use crate::freq::MatsubaraFreq;
11use crate::poly::{PiecewiseLegendrePoly, PiecewiseLegendrePolyVector};
12use crate::special_functions::spherical_bessel_j;
13use crate::traits::{Bosonic, Fermionic, Statistics, StatisticsType};
14
15/// Power model for asymptotic behavior
16#[derive(Debug, Clone)]
17pub struct PowerModel {
18    pub(crate) moments: Vec<f64>,
19}
20
21impl PowerModel {
22    /// Create a new power model with given moments
23    pub fn new(moments: Vec<f64>) -> Self {
24        Self { moments }
25    }
26}
27
28/// Piecewise Legendre polynomial with Fourier transform capability
29///
30/// This represents a piecewise Legendre polynomial that can be evaluated
31/// in the Matsubara frequency domain using Fourier transform.
32#[derive(Debug, Clone)]
33pub struct PiecewiseLegendreFT<S: StatisticsType> {
34    /// The underlying piecewise Legendre polynomial
35    pub(crate) poly: PiecewiseLegendrePoly,
36    /// Asymptotic cutoff frequency index
37    pub(crate) n_asymp: f64,
38    /// Power model for asymptotic behavior
39    pub(crate) model: PowerModel,
40    _phantom: std::marker::PhantomData<S>,
41}
42
43// Type aliases for convenience
44pub type FermionicPiecewiseLegendreFT = PiecewiseLegendreFT<Fermionic>;
45pub type BosonicPiecewiseLegendreFT = PiecewiseLegendreFT<Bosonic>;
46
47impl<S: StatisticsType> PiecewiseLegendreFT<S> {
48    /// Create a new PiecewiseLegendreFT from a polynomial and statistics
49    ///
50    /// # Arguments
51    /// * `poly` - The underlying piecewise Legendre polynomial
52    /// * `stat` - Statistics type (Fermionic or Bosonic)
53    /// * `n_asymp` - Asymptotic cutoff frequency index (default: infinity)
54    ///
55    /// # Errors
56    ///
57    /// * [`Error::InvalidParameter`] if the domain of `poly` is not [-1, 1]
58    ///   within 1e-12, or `n_asymp` is NaN or negative (infinity, the
59    ///   default, disables the asymptotic model)
60    pub fn new(poly: PiecewiseLegendrePoly, _stat: S, n_asymp: Option<f64>) -> Result<Self, Error> {
61        // Validate domain
62        if !((poly.xmin + 1.0).abs() <= 1e-12 && (poly.xmax - 1.0).abs() <= 1e-12) {
63            return Err(Error::InvalidParameter {
64                name: "poly",
65                value: format!("a polynomial on [{:?}, {:?}]", poly.xmin, poly.xmax),
66                reason: "must be defined on [-1, 1] (within 1e-12)".to_string(),
67            });
68        }
69        if let Some(n) = n_asymp.filter(|n| !(*n >= 0.0)) {
70            return Err(Error::InvalidParameter {
71                name: "n_asymp",
72                value: format!("{n:?}"),
73                reason: "must be non-negative (infinity disables the asymptotic model)".to_string(),
74            });
75        }
76
77        let n_asymp = n_asymp.unwrap_or(f64::INFINITY);
78        let model = Self::power_model(&poly);
79
80        Ok(Self {
81            poly,
82            n_asymp,
83            model,
84            _phantom: std::marker::PhantomData,
85        })
86    }
87
88    /// Get the asymptotic cutoff frequency index
89    pub fn get_n_asymp(&self) -> f64 {
90        self.n_asymp
91    }
92
93    /// Get the statistics type
94    pub fn get_statistics(&self) -> Statistics {
95        S::STATISTICS
96    }
97
98    /// Get the zeta value for this statistics type
99    pub fn zeta(&self) -> i64 {
100        match S::STATISTICS {
101            Statistics::Fermionic => 1,
102            Statistics::Bosonic => 0,
103        }
104    }
105
106    /// Get a reference to the underlying polynomial
107    pub fn get_poly(&self) -> &PiecewiseLegendrePoly {
108        &self.poly
109    }
110
111    /// Evaluate the Fourier transform at a Matsubara frequency
112    ///
113    /// # Arguments
114    /// * `omega` - Matsubara frequency
115    ///
116    /// # Returns
117    /// The complex Fourier transform value
118    pub fn evaluate(&self, omega: &MatsubaraFreq<S>) -> Complex64 {
119        let n = omega.get_n();
120        if (n as f64).abs() < self.n_asymp {
121            self.compute_unl_inner(&self.poly, n)
122        } else {
123            self.giw(n)
124        }
125    }
126
127    /// Evaluate at integer Matsubara index
128    ///
129    /// # Errors
130    ///
131    /// [`Error::InvalidMatsubaraIndex`] if `n` has the wrong parity for the
132    /// statistics (it returned 0 before)
133    pub fn evaluate_at_n(&self, n: i64) -> Result<Complex64, Error> {
134        Ok(self.evaluate(&MatsubaraFreq::<S>::new(n)?))
135    }
136
137    /// Evaluate at multiple Matsubara indices
138    ///
139    /// # Errors
140    ///
141    /// The error of [`Self::evaluate_at_n`] for the first invalid index
142    pub fn evaluate_at_ns(&self, ns: &[i64]) -> Result<Vec<Complex64>, Error> {
143        ns.iter().map(|&n| self.evaluate_at_n(n)).collect()
144    }
145
146    /// Create power model for asymptotic behavior
147    fn power_model(poly: &PiecewiseLegendrePoly) -> PowerModel {
148        // Evaluate at the right end of the domain: `new` accepts an xmax
149        // within 1e-12 of 1, where derivs(1.0) would be outside the domain.
150        let deriv_x1 = poly.derivs(poly.xmax);
151        let moments = Self::power_moments(&deriv_x1, poly.l);
152        PowerModel::new(moments)
153    }
154
155    /// Compute power moments for asymptotic expansion
156    ///
157    /// Integration by parts gives moments that combine the derivatives of the
158    /// function at x = 1 and at x = -1. Only those at x = 1 (`deriv_x1`) are
159    /// used; the ones at x = -1 follow from the parity `(-1)^l`. `l` must
160    /// therefore be the index of the singular function in its SVE result (see
161    /// `PiecewiseLegendrePoly::l`), not its index within an even or odd block.
162    fn power_moments(deriv_x1: &[f64], l: i32) -> Vec<f64> {
163        let statsign = match S::STATISTICS {
164            Statistics::Fermionic => -1.0,
165            Statistics::Bosonic => 1.0,
166        };
167
168        let mut moments = deriv_x1.to_vec();
169        for (m, moment) in moments.iter_mut().enumerate() {
170            let m_f64 = (m + 1) as f64; // Julia uses 1-based indexing
171            *moment *=
172                -(statsign * (-1.0_f64).powi(m_f64 as i32) + (-1.0_f64).powi(l)) / 2.0_f64.sqrt();
173        }
174        moments
175    }
176
177    /// Compute the inner Fourier transform (for small frequencies)
178    fn compute_unl_inner(&self, poly: &PiecewiseLegendrePoly, wn: i64) -> Complex64 {
179        let wred = PI / 4.0 * wn as f64;
180        let phase_wi = Self::phase_stable(poly, wn);
181        let mut res = Complex64::new(0.0, 0.0);
182
183        let order_max = poly.polyorder;
184        let segment_count = poly.knots.len() - 1;
185
186        for order in 0..order_max {
187            for j in 0..segment_count {
188                let data_oj = poly.data[[order, j]];
189                let tnl = Self::get_tnl(order as i32, wred * poly.delta_x[j]);
190                res += data_oj * tnl * phase_wi[j] / poly.norms[j];
191            }
192        }
193
194        res / 2.0_f64.sqrt()
195    }
196
197    /// Compute asymptotic behavior (for large frequencies)
198    fn giw(&self, wn: i64) -> Complex64 {
199        let iw = Complex64::new(0.0, PI / 2.0 * wn as f64);
200        if wn == 0 {
201            return Complex64::new(0.0, 0.0);
202        }
203
204        let inv_iw = 1.0 / iw;
205
206        inv_iw * Self::evalpoly(inv_iw, &self.model.moments)
207    }
208
209    /// Evaluate polynomial at complex point (Horner's method)
210    fn evalpoly(x: Complex64, coeffs: &[f64]) -> Complex64 {
211        let mut result = Complex64::new(0.0, 0.0);
212        for i in (0..coeffs.len()).rev() {
213            result = result * x + Complex64::new(coeffs[i], 0.0);
214        }
215        result
216    }
217
218    /// Compute midpoint relative to nearest integer
219    ///
220    /// Returns (xmid_diff, extra_shift) where:
221    /// - xmid_diff: midpoint values for numerical stability
222    /// - extra_shift: nearest integer shift (-1, 0, or 1)
223    fn shift_xmid(knots: &[f64], delta_x: &[f64]) -> (Vec<f64>, Vec<i32>) {
224        let n_segments = delta_x.len();
225        let delta_x_half: Vec<f64> = delta_x.iter().map(|&dx| dx / 2.0).collect();
226
227        // xmid_m1: cumsum(Δx) - Δx_half
228        let mut xmid_m1 = Vec::with_capacity(n_segments);
229        let mut cumsum = 0.0;
230        for i in 0..n_segments {
231            cumsum += delta_x[i];
232            xmid_m1.push(cumsum - delta_x_half[i]);
233        }
234
235        // xmid_p1: -reverse(cumsum(reverse(Δx))) + Δx_half
236        let mut xmid_p1 = Vec::with_capacity(n_segments);
237        let mut cumsum_rev = 0.0;
238        for i in (0..n_segments).rev() {
239            cumsum_rev += delta_x[i];
240            xmid_p1.insert(0, -cumsum_rev + delta_x_half[i]);
241        }
242
243        // xmid_0: knots[1:] - Δx_half
244        let xmid_0: Vec<f64> = (0..n_segments)
245            .map(|i| knots[i + 1] - delta_x_half[i])
246            .collect();
247
248        // Determine shift and diff
249        let mut xmid_diff = Vec::with_capacity(n_segments);
250        let mut extra_shift = Vec::with_capacity(n_segments);
251
252        for i in 0..n_segments {
253            let shift = xmid_0[i].round() as i32;
254            extra_shift.push(shift);
255
256            // Choose appropriate xmid based on shift
257            let diff = match shift {
258                -1 => xmid_m1[i],
259                0 => xmid_0[i],
260                1 => xmid_p1[i],
261                _ => xmid_0[i], // Fallback
262            };
263            xmid_diff.push(diff);
264        }
265
266        (xmid_diff, extra_shift)
267    }
268
269    /// Compute stable phase factors
270    ///
271    /// Computes: im^mod(wn * (extra_shift + 1), 4) * cispi(wn * xmid_diff / 2)
272    /// where cispi(x) = exp(i*π*x)
273    fn phase_stable(poly: &PiecewiseLegendrePoly, wn: i64) -> Vec<Complex64> {
274        let (xmid_diff, extra_shift) = Self::shift_xmid(&poly.knots, &poly.delta_x);
275        let mut phase_wi = Vec::with_capacity(xmid_diff.len());
276
277        let im_unit = Complex64::new(0.0, 1.0);
278
279        for j in 0..xmid_diff.len() {
280            // Compute im^mod(wn * (extra_shift[j] + 1), 4)
281            let power = ((wn * (extra_shift[j] as i64 + 1)) % 4 + 4) % 4; // Ensure positive mod
282            let im_power = im_unit.powi(power as i32);
283
284            // Compute cispi(wn * xmid_diff[j] / 2) = exp(i*π*wn*xmid_diff/2)
285            let arg = PI * wn as f64 * xmid_diff[j] / 2.0;
286            let cispi = Complex64::new(arg.cos(), arg.sin());
287
288            phase_wi.push(im_power * cispi);
289        }
290
291        phase_wi
292    }
293
294    /// Get T_nl coefficient (special function)
295    ///
296    /// This implements the T_nl function which is related to spherical Bessel functions:
297    /// T_nl(w) = 2 * i^l * j_l(|w|) * (w < 0 ? conj : identity)
298    /// where j_l is the spherical Bessel function of the first kind.
299    pub fn get_tnl(l: i32, w: f64) -> Complex64 {
300        let abs_w = w.abs();
301
302        // Use the high-precision spherical Bessel function from special_functions
303        let sph_bessel = spherical_bessel_j(l, abs_w);
304
305        // Compute 2 * i^l
306        let im_unit = Complex64::new(0.0, 1.0);
307        let im_power = im_unit.powi(l);
308        let result = 2.0 * im_power * sph_bessel;
309
310        // Apply conjugation for negative w
311        if w < 0.0 { result.conj() } else { result }
312    }
313}
314
315/// Vector of PiecewiseLegendreFT polynomials
316#[derive(Debug, Clone)]
317pub struct PiecewiseLegendreFTVector<S: StatisticsType> {
318    pub(crate) polyvec: Vec<PiecewiseLegendreFT<S>>,
319    _phantom: std::marker::PhantomData<S>,
320}
321
322// Type aliases for convenience
323pub type FermionicPiecewiseLegendreFTVector = PiecewiseLegendreFTVector<Fermionic>;
324pub type BosonicPiecewiseLegendreFTVector = PiecewiseLegendreFTVector<Bosonic>;
325
326impl<S: StatisticsType> PiecewiseLegendreFTVector<S> {
327    /// Create an empty vector
328    pub fn new() -> Self {
329        Self {
330            polyvec: Vec::new(),
331            _phantom: std::marker::PhantomData,
332        }
333    }
334
335    /// Create from a vector of PiecewiseLegendreFT
336    pub fn from_vector(polyvec: Vec<PiecewiseLegendreFT<S>>) -> Self {
337        Self {
338            polyvec,
339            _phantom: std::marker::PhantomData,
340        }
341    }
342
343    /// Get the number of polynomials in the vector
344    pub fn len(&self) -> usize {
345        self.polyvec.len()
346    }
347
348    /// Check if the vector is empty
349    pub fn is_empty(&self) -> bool {
350        self.polyvec.is_empty()
351    }
352
353    /// Create from PiecewiseLegendrePolyVector and statistics
354    ///
355    /// # Errors
356    ///
357    /// The errors of [`PiecewiseLegendreFT::new`] for the first polynomial
358    /// that it rejects
359    pub fn from_poly_vector(
360        polys: &PiecewiseLegendrePolyVector,
361        _stat: S,
362        n_asymp: Option<f64>,
363    ) -> Result<Self, Error> {
364        let mut polyvec = Vec::with_capacity(polys.size());
365
366        for poly in polys.get_polys() {
367            let ft_poly = PiecewiseLegendreFT::new(poly.clone(), _stat, n_asymp)?;
368            polyvec.push(ft_poly);
369        }
370
371        Ok(Self {
372            polyvec,
373            _phantom: std::marker::PhantomData,
374        })
375    }
376
377    /// Get the size of the vector
378    pub fn size(&self) -> usize {
379        self.polyvec.len()
380    }
381
382    /// The transformed polynomials, in order
383    pub fn get_polys(&self) -> &[PiecewiseLegendreFT<S>] {
384        &self.polyvec
385    }
386
387    /// Get element by index (immutable)
388    pub fn get(&self, index: usize) -> Option<&PiecewiseLegendreFT<S>> {
389        self.polyvec.get(index)
390    }
391
392    /// Get element by index (mutable)
393    pub fn get_mut(&mut self, index: usize) -> Option<&mut PiecewiseLegendreFT<S>> {
394        self.polyvec.get_mut(index)
395    }
396
397    /// Set element at index
398    ///
399    /// # Errors
400    /// [`Error::InvalidParameter`] if `index` is not less than the length;
401    /// the vector is then unchanged
402    pub fn set(&mut self, index: usize, poly: PiecewiseLegendreFT<S>) -> Result<(), Error> {
403        let len = self.polyvec.len();
404        if index >= len {
405            return Err(Error::InvalidParameter {
406                name: "index",
407                value: index.to_string(),
408                reason: format!("must be less than the length {len}"),
409            });
410        }
411        self.polyvec[index] = poly;
412        Ok(())
413    }
414
415    /// Create a similar empty vector
416    pub fn similar(&self) -> Self {
417        Self::new()
418    }
419
420    /// Get n_asymp from the first element (if any)
421    pub fn n_asymp(&self) -> f64 {
422        self.polyvec.first().map_or(f64::INFINITY, |p| p.n_asymp)
423    }
424
425    /// Evaluate all polynomials at a Matsubara frequency
426    pub fn evaluate_at(&self, omega: &MatsubaraFreq<S>) -> Vec<Complex64> {
427        self.polyvec
428            .iter()
429            .map(|poly| poly.evaluate(omega))
430            .collect()
431    }
432
433    /// Evaluate all polynomials at multiple Matsubara frequencies
434    pub fn evaluate_at_many(&self, omegas: &[MatsubaraFreq<S>]) -> Vec<Vec<Complex64>> {
435        omegas.iter().map(|omega| self.evaluate_at(omega)).collect()
436    }
437}
438
439// Indexing operators
440impl<S: StatisticsType> std::ops::Index<usize> for PiecewiseLegendreFTVector<S> {
441    type Output = PiecewiseLegendreFT<S>;
442
443    fn index(&self, index: usize) -> &Self::Output {
444        &self.polyvec[index]
445    }
446}
447
448impl<S: StatisticsType> std::ops::IndexMut<usize> for PiecewiseLegendreFTVector<S> {
449    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
450        &mut self.polyvec[index]
451    }
452}
453
454// Default implementations
455impl<S: StatisticsType> Default for PiecewiseLegendreFTVector<S> {
456    fn default() -> Self {
457        Self::new()
458    }
459}
460
461// ===== Matsubara sampling point selection =====
462
463/// Default grid for finding extrema/sign changes
464/// Matches C++ DEFAULT_GRID: [0:2^6-1] followed by exponential spacing up to 2^25
465/// Generated from Julia: [range(0; length=2^6); trunc.(Int, exp2.(range(6, 25; length=32 * (25 - 6) + 1)))]
466const DEFAULT_GRID: &[i64] = &[
467    0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
468    26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
469    50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 68, 69, 71, 72, 74, 76, 77,
470    79, 81, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100, 103, 105, 107, 109, 112, 114, 117, 119, 122,
471    125, 128, 130, 133, 136, 139, 142, 145, 148, 152, 155, 158, 162, 165, 169, 173, 177, 181, 184,
472    189, 193, 197, 201, 206, 210, 215, 219, 224, 229, 234, 239, 245, 250, 256, 261, 267, 273, 279,
473    285, 291, 297, 304, 311, 317, 324, 331, 339, 346, 354, 362, 369, 378, 386, 394, 403, 412, 421,
474    430, 439, 449, 459, 469, 479, 490, 501, 512, 523, 534, 546, 558, 570, 583, 595, 608, 622, 635,
475    649, 663, 678, 693, 708, 724, 739, 756, 772, 789, 806, 824, 842, 861, 879, 899, 918, 939, 959,
476    980, 1002, 1024, 1046, 1069, 1092, 1116, 1141, 1166, 1191, 1217, 1244, 1271, 1299, 1327, 1357,
477    1386, 1417, 1448, 1479, 1512, 1545, 1579, 1613, 1649, 1685, 1722, 1759, 1798, 1837, 1878, 1919,
478    1961, 2004, 2048, 2092, 2138, 2185, 2233, 2282, 2332, 2383, 2435, 2488, 2543, 2599, 2655, 2714,
479    2773, 2834, 2896, 2959, 3024, 3090, 3158, 3227, 3298, 3370, 3444, 3519, 3596, 3675, 3756, 3838,
480    3922, 4008, 4096, 4185, 4277, 4371, 4466, 4564, 4664, 4766, 4870, 4977, 5086, 5198, 5311, 5428,
481    5547, 5668, 5792, 5919, 6049, 6181, 6316, 6455, 6596, 6741, 6888, 7039, 7193, 7351, 7512, 7676,
482    7844, 8016, 8192, 8371, 8554, 8742, 8933, 9129, 9328, 9533, 9741, 9955, 10173, 10396, 10623,
483    10856, 11094, 11336, 11585, 11838, 12098, 12363, 12633, 12910, 13193, 13482, 13777, 14078,
484    14387, 14702, 15024, 15353, 15689, 16032, 16384, 16742, 17109, 17484, 17866, 18258, 18657,
485    19066, 19483, 19910, 20346, 20792, 21247, 21712, 22188, 22673, 23170, 23677, 24196, 24726,
486    25267, 25820, 26386, 26964, 27554, 28157, 28774, 29404, 30048, 30706, 31378, 32065, 32768,
487    33485, 34218, 34968, 35733, 36516, 37315, 38132, 38967, 39821, 40693, 41584, 42494, 43425,
488    44376, 45347, 46340, 47355, 48392, 49452, 50535, 51641, 52772, 53928, 55108, 56315, 57548,
489    58809, 60096, 61412, 62757, 64131, 65536, 66971, 68437, 69936, 71467, 73032, 74631, 76265,
490    77935, 79642, 81386, 83168, 84989, 86850, 88752, 90695, 92681, 94711, 96785, 98904, 101070,
491    103283, 105545, 107856, 110217, 112631, 115097, 117618, 120193, 122825, 125514, 128263, 131072,
492    133942, 136875, 139872, 142935, 146064, 149263, 152531, 155871, 159284, 162772, 166337, 169979,
493    173701, 177504, 181391, 185363, 189422, 193570, 197809, 202140, 206566, 211090, 215712, 220435,
494    225262, 230195, 235236, 240387, 245650, 251029, 256526, 262144, 267884, 273750, 279744, 285870,
495    292129, 298526, 305063, 311743, 318569, 325545, 332674, 339958, 347402, 355009, 362783, 370727,
496    378845, 387141, 395618, 404281, 413133, 422180, 431424, 440871, 450525, 460390, 470472, 480774,
497    491301, 502059, 513053, 524288, 535768, 547500, 559488, 571740, 584259, 597053, 610126, 623487,
498    637139, 651091, 665348, 679917, 694805, 710019, 725567, 741455, 757690, 774282, 791236, 808562,
499    826267, 844360, 862849, 881743, 901051, 920781, 940944, 961548, 982603, 1004119, 1026107,
500    1048576, 1071536, 1095000, 1118977, 1143480, 1168519, 1194106, 1220253, 1246974, 1274279,
501    1302182, 1330696, 1359834, 1389611, 1420039, 1451134, 1482910, 1515381, 1548564, 1582473,
502    1617125, 1652535, 1688721, 1725699, 1763487, 1802102, 1841563, 1881888, 1923096, 1965207,
503    2008239, 2052214, 2097152, 2143073, 2190000, 2237955, 2286960, 2337038, 2388212, 2440507,
504    2493948, 2548558, 2604364, 2661392, 2719669, 2779222, 2840079, 2902269, 2965820, 3030763,
505    3097128, 3164947, 3234250, 3305071, 3377443, 3451399, 3526975, 3604205, 3683127, 3763777,
506    3846193, 3930414, 4016479, 4104428, 4194304, 4286147, 4380001, 4475911, 4573920, 4674076,
507    4776425, 4881015, 4987896, 5097116, 5208729, 5322785, 5439339, 5558445, 5680159, 5804538,
508    5931641, 6061527, 6194257, 6329894, 6468501, 6610142, 6754886, 6902798, 7053950, 7208411,
509    7366255, 7527555, 7692387, 7860828, 8032958, 8208857, 8388608, 8572294, 8760003, 8951822,
510    9147841, 9348153, 9552851, 9762031, 9975792, 10194233, 10417458, 10645571, 10878678, 11116890,
511    11360318, 11609077, 11863283, 12123055, 12388515, 12659788, 12937002, 13220285, 13509772,
512    13805597, 14107900, 14416823, 14732510, 15055110, 15384774, 15721657, 16065917, 16417714,
513    16777216, 17144589, 17520006, 17903645, 18295683, 18696307, 19105702, 19524063, 19951584,
514    20388467, 20834916, 21291142, 21757357, 22233781, 22720637, 23218155, 23726566, 24246110,
515    24777031, 25319577, 25874004, 26440571, 27019544, 27611195, 28215801, 28833647, 29465021,
516    30110221, 30769549, 31443315, 32131834, 32835429, 33554432,
517];
518
519/// Find sign changes of a Matsubara basis function
520///
521/// Returns Matsubara frequencies where the function changes sign.
522///
523/// # Errors
524///
525/// [`Error::NotSupported`] if `u_hat` has no definite parity (symm = 0, #183)
526pub(crate) fn sign_changes<S: StatisticsType + 'static>(
527    u_hat: &PiecewiseLegendreFT<S>,
528    positive_only: bool,
529) -> Result<Vec<MatsubaraFreq<S>>, Error> {
530    let f = func_for_part(u_hat)?;
531    let x0 = find_all(&f, DEFAULT_GRID);
532
533    // Convert to Matsubara indices: n = 2*x + zeta
534    let mut indices: Vec<i64> = x0.iter().map(|&x| 2 * x + u_hat.zeta()).collect();
535
536    if !positive_only {
537        symmetrize_matsubara_inplace(&mut indices);
538    }
539
540    Ok(indices
541        .iter()
542        .filter_map(|&n| MatsubaraFreq::<S>::new(n).ok())
543        .collect())
544}
545
546/// Find extrema of a Matsubara basis function
547///
548/// Returns Matsubara frequencies where the function has local extrema.
549///
550/// # Errors
551///
552/// [`Error::NotSupported`] if `u_hat` has no definite parity (symm = 0, #183)
553pub(crate) fn find_extrema<S: StatisticsType + 'static>(
554    u_hat: &PiecewiseLegendreFT<S>,
555    positive_only: bool,
556) -> Result<Vec<MatsubaraFreq<S>>, Error> {
557    let f = func_for_part(u_hat)?;
558    let x0 = discrete_extrema(&f, DEFAULT_GRID);
559
560    // Convert to Matsubara indices: n = 2*x + zeta
561    let mut indices: Vec<i64> = x0.iter().map(|&x| 2 * x + u_hat.zeta()).collect();
562
563    if !positive_only {
564        symmetrize_matsubara_inplace(&mut indices);
565    }
566
567    Ok(indices
568        .iter()
569        .filter_map(|&n| MatsubaraFreq::<S>::new(n).ok())
570        .collect())
571}
572
573/// The part of a Matsubara basis function whose sign changes choose the
574/// default sampling points
575#[derive(Clone, Copy)]
576enum Part {
577    Real,
578    Imaginary,
579}
580
581/// The part of `û(iν)` that carries its sign changes, from the parity `symm`
582/// of `u` and the statistics: for an even `u` (symm = 1) the fermionic `û`
583/// is imaginary and the bosonic `û` real, for an odd `u` (symm = -1) the
584/// other way round
585///
586/// # Errors
587///
588/// [`Error::NotSupported`] if `symm` is not ±1: the singular functions of an
589/// SVE that is not centrosymmetric have no parity (#183).
590fn part_with_the_sign_changes<S: StatisticsType>(symm: i32) -> Result<Part, Error> {
591    match (symm, S::STATISTICS) {
592        (1, Statistics::Bosonic) | (-1, Statistics::Fermionic) => Ok(Part::Real),
593        (1, Statistics::Fermionic) | (-1, Statistics::Bosonic) => Ok(Part::Imaginary),
594        _ => Err(Error::NotSupported {
595            what: format!(
596                "default Matsubara sampling points, which need basis functions of definite \
597                 parity (symm = ±1, from a centrosymmetric SVE); got symm = {symm}"
598            ),
599        }),
600    }
601}
602
603impl Part {
604    fn of(self, value: Complex64) -> f64 {
605        match self {
606            Part::Real => value.re,
607            Part::Imaginary => value.im,
608        }
609    }
610}
611
612/// Create a function that extracts the part of `poly_ft` that carries its
613/// sign changes (see [`part_with_the_sign_changes`])
614fn func_for_part<S: StatisticsType + 'static>(
615    poly_ft: &PiecewiseLegendreFT<S>,
616) -> Result<Box<dyn Fn(i64) -> f64>, Error> {
617    let part = part_with_the_sign_changes::<S>(poly_ft.poly.symm())?;
618    let zeta = poly_ft.zeta();
619    let poly_ft_clone = poly_ft.clone();
620
621    Ok(Box::new(move |n: i64| {
622        let omega = MatsubaraFreq::<S>::new(2 * n + zeta).unwrap();
623        part.of(poly_ft_clone.evaluate(&omega))
624    }))
625}
626
627/// Integer bisection: find the zero crossing of f in [a, b] where f(a) and f(b)
628/// have different signs. Returns the largest integer x such that f(x) has the
629/// same sign as f(a). Matches Julia's `bisect` implementation.
630fn bisect(f: &dyn Fn(i64) -> f64, a: i64, b: i64, fa: f64) -> i64 {
631    let mut lo = a;
632    let mut hi = b;
633    let mut flo = fa;
634
635    while (hi - lo).abs() > 1 {
636        let mid = lo + (hi - lo) / 2;
637        let fmid = f(mid);
638        if flo.signum() != fmid.signum() {
639            hi = mid;
640        } else {
641            lo = mid;
642            flo = fmid;
643        }
644    }
645    lo
646}
647
648/// Integer bisection for extrema: find the discrete extremum of f in [a, b].
649/// Returns the integer x in [a, b] that maximizes |f(x)|.
650/// Matches Julia's `bisect_discr_extremum` implementation.
651fn bisect_discr_extremum(f: &dyn Fn(i64) -> f64, a: i64, b: i64) -> i64 {
652    let mut lo = a;
653    let mut hi = b;
654
655    while hi - lo > 2 {
656        let mid1 = lo + (hi - lo) / 3;
657        let mid2 = hi - (hi - lo) / 3;
658        if f(mid1).abs() < f(mid2).abs() {
659            lo = mid1;
660        } else {
661            hi = mid2;
662        }
663    }
664
665    // Check all remaining candidates
666    let mut best = lo;
667    let mut best_val = f(lo).abs();
668    for x in (lo + 1)..=hi {
669        let val = f(x).abs();
670        if val > best_val {
671            best = x;
672            best_val = val;
673        }
674    }
675    best
676}
677
678/// Find all sign changes of a function on a grid, using bisection to locate
679/// the precise integer position within each grid interval.
680/// Matches Julia's `find_all` with `bisect`.
681fn find_all(f: &dyn Fn(i64) -> f64, xgrid: &[i64]) -> Vec<i64> {
682    if xgrid.is_empty() {
683        return Vec::new();
684    }
685
686    let mut results = Vec::new();
687    let mut prev_val = f(xgrid[0]);
688
689    for i in 1..xgrid.len() {
690        let val = f(xgrid[i]);
691        // Detect sign change (both values must be nonzero)
692        if prev_val.signum() != val.signum() && prev_val != 0.0 && val != 0.0 {
693            // Bisect to find precise integer location of the zero
694            let root = bisect(f, xgrid[i - 1], xgrid[i], prev_val);
695            results.push(root);
696        }
697        prev_val = val;
698    }
699
700    results
701}
702
703/// Find discrete extrema of a function on a grid, using bisection to locate
704/// the precise integer position within each grid interval.
705/// Also checks boundary extrema. Matches Julia v1 / Python 1.x `discrete_extrema`.
706fn discrete_extrema(f: &dyn Fn(i64) -> f64, xgrid: &[i64]) -> Vec<i64> {
707    if xgrid.len() < 3 {
708        return Vec::new();
709    }
710
711    let fx: Vec<f64> = xgrid.iter().map(|&x| f(x)).collect();
712    let absfx: Vec<f64> = fx.iter().map(|v| v.abs()).collect();
713
714    // Forward differences: gx[i] = fx[i+1] - fx[i]
715    // signdfdx[i] = signbit(gx[i])
716    // A derivative sign change at index i means the extremum is strictly
717    // between xgrid[i] and xgrid[i+2].
718    let signdfdx: Vec<bool> = fx
719        .windows(2)
720        .map(|w| (w[1] - w[0]).is_sign_negative())
721        .collect();
722
723    let mut results = Vec::new();
724
725    for i in 0..signdfdx.len() - 1 {
726        if signdfdx[i] != signdfdx[i + 1] {
727            // Extremum between xgrid[i] and xgrid[i+2]
728            let refined = bisect_discr_extremum(f, xgrid[i], xgrid[i + 2]);
729            results.push(refined);
730        }
731    }
732
733    // Boundary: first point is extremum if |f| decreases or sign changes inwards
734    let sfx: Vec<bool> = fx.iter().map(|v| v.is_sign_negative()).collect();
735    if absfx[0] > absfx[1] || sfx[0] != sfx[1] {
736        results.insert(0, xgrid[0]);
737    }
738
739    // Boundary: last point is extremum if |f| decreases or sign changes inwards
740    let n = fx.len();
741    if absfx[n - 1] > absfx[n - 2] || sfx[n - 1] != sfx[n - 2] {
742        results.push(xgrid[n - 1]);
743    }
744
745    results
746}
747
748/// Symmetrize Matsubara indices by adding negative frequencies
749fn symmetrize_matsubara_inplace(xs: &mut Vec<i64>) {
750    if xs.is_empty() {
751        return;
752    }
753
754    // Remove zero if present
755    xs.retain(|&x| x != 0);
756
757    // Add negative frequencies
758    let positives: Vec<i64> = xs.iter().filter(|&&x| x > 0).copied().collect();
759    let mut negatives: Vec<i64> = positives.iter().map(|&x| -x).collect();
760
761    xs.append(&mut negatives);
762    xs.sort();
763    xs.dedup();
764}
765
766#[cfg(test)]
767#[path = "polyfourier_tests.rs"]
768mod polyfourier_tests;
sparse_ir_basis/polyfourier.rs

sparse_ir_basis/
polyfourier.rs
