multishoot 1.0.0

//! The periodic direction: collocation grids in theta, the spectral derivative d/dtheta on them, and sines and
//! cosines in double-double, and in octo-double for constants that must be exact in every precision.
//!
//! A field is held on M = 8 B points. Without symmetry ([`Parity::Full`]) the grid is theta_j = 2 pi j / M over the
//! period. A field odd under theta -> theta + pi (odd harmonics only) or even (even harmonics only) is held on the
//! half grid theta_j = pi j / M, j < M, of the 2M-point grid; the other half follows from the symmetry, and the
//! derivative matrix is the full one folded with it. Pointwise products keep the parities (odd times odd is even).
//!
//! The derivative matrix of the N-point grid (N even) is D_ij = (-1)^(i-j) cot((i-j) pi / N) / 2, i != j, and zero
//! on the diagonal: the derivative of the trigonometric interpolant, exact for harmonics below N/2 (Trefethen,
//! Spectral Methods in MATLAB, SIAM 2000, ch. 3).
use crate::lane::DMat;
use multidouble::dd::DD;
use multidouble::od::OD;

/// pi in double-double
pub const PI: DD = DD {
    hi: std::f64::consts::PI,
    lo: 1.2246467991473532e-16,
};

/// pi / 2 in double-double
const PI_2: DD = DD {
    hi: std::f64::consts::FRAC_PI_2,
    lo: 6.123233995736766e-17,
};

/// The symmetry of a field under theta -> theta + pi, which fixes its grid and derivative matrix
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Parity {
    /// odd harmonics only, on the half grid
    Odd,
    /// even harmonics only, on the half grid
    Even,
    /// all harmonics, on the full grid
    Full,
}

/// sin and cos of r by their Taylor series, |r| <= pi / 4 (+ rounding): terms to below 2^-115 of the first
fn sin_cos_series(r: DD) -> (DD, DD) {
    let r2 = r * r;
    let (mut s, mut c) = (r, DD::ONE);
    let (mut ts, mut tc) = (r, DD::ONE);
    for k in 1..30 {
        let kf = k as f64;
        ts = -(ts * r2) / DD::from_f64((2.0 * kf) * (2.0 * kf + 1.0));
        tc = -(tc * r2) / DD::from_f64((2.0 * kf - 1.0) * (2.0 * kf));
        s = s + ts;
        c = c + tc;
        if tc.hi.abs() < 2e-35 && ts.hi.abs() < 2e-35 * r.hi.abs().max(1e-300) {
            break;
        }
    }
    (s, c)
}

/// (sin, cos) of r + n pi / 2 from those of r
fn quadrant(n: i64, (s, c): (DD, DD)) -> (DD, DD) {
    match n.rem_euclid(4) {
        0 => (s, c),
        1 => (c, -s),
        2 => (-s, -c),
        _ => (-c, s),
    }
}

/// sin x and cos x in double-double: x reduced by the nearest multiple of pi / 2 (an absolute error of about
/// |x| 2^-106 from pi's rounding), then the Taylor series
pub fn sin_cos(x: DD) -> (DD, DD) {
    let n = (x.hi / std::f64::consts::FRAC_PI_2).round();
    let r = x - PI_2.mul_f(n);
    quadrant(n as i64, sin_cos_series(r))
}

/// sin and cos of (p / q) pi, q > 0, with the reduction in integers: the remainder (2p - n q) pi / (2q) after the
/// nearest multiple n of pi / 2 is formed from an exact integer, so the result is accurate to about 2^-106 relative
/// for every p
pub fn sin_cos_pi(p: i64, q: i64) -> (DD, DD) {
    assert!(q > 0, "sin_cos_pi: q = {q}");
    // n = round(2p / q), the nearest multiple of pi / 2
    let n = (4 * p + q).div_euclid(2 * q);
    let num = 2 * p - n * q;
    let r = PI.mul_f(num as f64) / DD::from_f64((2 * q) as f64);
    quadrant(n, sin_cos_series(r))
}

/// pi in octo-double: the parts of its greedy round-to-nearest expansion, from 150 decimal digits (checked against
/// `OD::parse` and Machin's formula in the tests)
pub const PI_OD: OD = OD {
    x: [
        std::f64::consts::PI,
        1.2246467991473532e-16,
        -2.9947698097183397e-33,
        1.1124542208633653e-49,
        5.672231979640316e-66,
        1.7449862161352486e-83,
        6.02937273224954e-100,
        1.91012354687999e-116,
    ],
};

/// sin and cos of r in octo-double by their Taylor series, |r| <= pi / 4 (+ rounding): terms to below tol (1e-136
/// for octo-double)
fn sin_cos_series_od(r: OD, tol: f64) -> (OD, OD) {
    let r2 = r * r;
    let (mut s, mut c) = (r, OD::from_f64(1.0));
    let (mut ts, mut tc) = (r, OD::from_f64(1.0));
    for k in 1..80 {
        let kf = k as f64;
        ts = -(ts * r2) / OD::from_f64((2.0 * kf) * (2.0 * kf + 1.0));
        tc = -(tc * r2) / OD::from_f64((2.0 * kf - 1.0) * (2.0 * kf));
        s = s + ts;
        c = c + tc;
        if tc.to_f64().abs() < tol && ts.to_f64().abs() < tol {
            break;
        }
    }
    (s, c)
}

/// sin and cos of (p / q) pi, q > 0, in octo-double, reduced in integers as [`sin_cos_pi`]: accurate to about 2^-424
/// relative for every p, so the constants made from them are exact in every precision of the solve
pub fn sin_cos_pi_od(p: i64, q: i64) -> (OD, OD) {
    assert!(q > 0, "sin_cos_pi_od: q = {q}");
    let n = (4 * p + q).div_euclid(2 * q);
    let num = 2 * p - n * q;
    let r = PI_OD.mul_f(num as f64) / OD::from_f64((2 * q) as f64);
    let (s, c) = sin_cos_series_od(r, 1e-136);
    match n.rem_euclid(4) {
        0 => (s, c),
        1 => (c, -s),
        2 => (-s, -c),
        _ => (-c, s),
    }
}

/// The number of grid points per period: 2M on the half grids, M on the full one
fn period_points(parity: Parity, m: usize) -> i64 {
    match parity {
        Parity::Full => m as i64,
        _ => 2 * m as i64,
    }
}

/// theta_j of the grid with m points: 2 pi j / m on the full grid, pi j / m on the half grids
pub fn theta(parity: Parity, m: usize, j: usize) -> DD {
    PI.mul_f(2.0 * j as f64) / DD::from_f64(period_points(parity, m) as f64)
}

/// sin and cos of k theta_j on the grid with m points, from the integer multiple of pi
pub fn sin_cos_grid(parity: Parity, m: usize, j: usize, k: i64) -> (DD, DD) {
    sin_cos_pi(2 * k * j as i64, period_points(parity, m))
}

/// entry (i, j) of the N-point derivative matrix, i != j mod N: (-1)^(i-j) cot((i-j) pi / N) / 2
fn d_full(d: i64, n: i64) -> DD {
    let (s, c) = sin_cos_pi(d, n);
    let v = (c / s).mul_f(0.5);
    if d.rem_euclid(2) == 0 {
        v
    } else {
        -v
    }
}

/// d/dtheta on the M = 8 B points of the grid of `parity`, in double-double: on the half grids the 2M-point matrix
/// folded with the symmetry, `D_ij = D2M_ij -+ D2M_i,j+M`
pub fn dmat<const B: usize>(parity: Parity) -> DMat<B> {
    let m = 8 * B;
    let n = period_points(parity, m);
    let mut hi = vec![0.0; m * m];
    let mut lo = vec![0.0; m * m];
    for i in 0..m {
        for j in 0..m {
            let mut v = DD::ZERO;
            let images: &[(usize, f64)] = match parity {
                Parity::Full => &[(j, 1.0)],
                Parity::Odd => &[(j, 1.0), (j + m, -1.0)],
                Parity::Even => &[(j, 1.0), (j + m, 1.0)],
            };
            for &(jj, sgn) in images {
                if jj != i {
                    v = v + d_full(i as i64 - jj as i64, n).mul_f(sgn);
                }
            }
            hi[i * m + j] = v.hi;
            lo[i * m + j] = v.lo;
        }
    }
    DMat::<B>::from_rows(&hi, &lo)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::lane::DV;
    use crate::prec::nan_max;
    use multidouble::field::Fl;

    /// sin_cos against sin_cos_pi at multiples of pi / 7 and against f64 at random points; pass: 1e-30 and 1e-15
    #[test]
    fn sin_cos_agree() {
        for p in -30..30 {
            let x = PI.mul_f(p as f64) / DD::from_f64(7.0);
            let (s, c) = sin_cos(x);
            let (s2, c2) = sin_cos_pi(p, 7);
            assert!((s - s2).hi.abs() < 1e-30 && (c - c2).hi.abs() < 1e-30, "p {p}");
            assert!((s.hi - x.hi.sin()).abs() < 1e-15 && (c.hi - x.hi.cos()).abs() < 1e-15);
            // sin^2 + cos^2 = 1
            assert!((s * s + c * c - DD::ONE).hi.abs() < 1e-31);
        }
    }

    /// D applied to sin(k theta) and cos(k theta), k below the Nyquist harmonic, against k cos and -k sin; pass:
    /// 1e-29 (D's entries are O(M), each rounded to 2^-106)
    #[test]
    fn dmat_differentiates_harmonics() {
        const B: usize = 2;
        let m = 8 * B;
        for (parity, ks) in [
            (Parity::Odd, vec![1i64, 3, 5, 15]),
            (Parity::Even, vec![0i64, 2, 8, 14]),
            (Parity::Full, vec![0i64, 1, 4, 7]),
        ] {
            let d = dmat::<B>(parity);
            for &k in &ks {
                let samples = |f: &dyn Fn(usize) -> DD| {
                    let v: Vec<DD> = (0..m).map(f).collect();
                    let (h, l) = multidouble::dd::split(&v);
                    Fl::<DV, B>::load(&h, &l)
                };
                let s = samples(&|j| sin_cos_grid(parity, m, j, k).0);
                let c = samples(&|j| sin_cos_grid(parity, m, j, k).1);
                let ds = d.apply(&s);
                let dc = d.apply(&c);
                for j in 0..m {
                    let kf = k as f64;
                    let e1 = (ds.get(j) - c.get(j).mul_f(kf)).hi.abs();
                    let e2 = (dc.get(j) + s.get(j).mul_f(kf)).hi.abs();
                    assert!(e1 < 1e-29 && e2 < 1e-29, "{parity:?} k {k} j {j}: {e1:e} {e2:e}");
                }
            }
        }
    }

    /// |a - b| in octo-double, as an f64
    fn od_diff(a: OD, b: OD) -> f64 {
        (a - b).to_f64().abs()
    }

    /// arctan(1 / n) in octo-double by its series, n >= 5: terms to below 1e-140
    fn atan_inv(n: f64) -> OD {
        let x = OD::from_f64(1.0) / OD::from_f64(n);
        let x2 = x * x;
        let (mut s, mut p) = (x, x);
        for k in 1..200 {
            p = -(p * x2);
            let t = p / OD::from_f64((2 * k + 1) as f64);
            s = s + t;
            if t.to_f64().abs() < 1e-140 {
                break;
            }
        }
        s
    }

    /// PI_OD against the parse of its 150 decimal digits (every part), and against Machin's formula
    /// pi = 16 arctan(1/5) - 4 arctan(1/239) in octo-double; pass: equal parts, and 1e-125
    #[test]
    fn pi_od_is_pi() {
        let parsed = OD::parse(
            "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813",
        )
        .unwrap();
        assert_eq!(parsed.x, PI_OD.x);
        let machin = atan_inv(5.0).mul_f(16.0) - atan_inv(239.0).mul_f(4.0);
        let e = od_diff(machin, PI_OD);
        println!("|Machin - PI_OD| = {e:.3e}");
        assert!(e < 1e-125, "{e:e}");
    }

    /// The control of `pi_od_is_pi`: pi in quad-double (its first four parts) must fail Machin's bound.
    #[test]
    fn pi_od_control_quad_parts_fail() {
        let mut q = PI_OD;
        q.x[4..].iter_mut().for_each(|v| *v = 0.0);
        let machin = atan_inv(5.0).mul_f(16.0) - atan_inv(239.0).mul_f(4.0);
        assert!(od_diff(machin, q) > 1e-125);
    }

    /// sin_cos_pi_od against sin_cos_pi in double-double (1e-31), sin^2 + cos^2 = 1 (1e-126), and the double angle
    /// sin 2x = 2 sin x cos x against the value at 2p (1e-125), for p / q over three periods and three q
    #[test]
    fn sin_cos_pi_od_agrees() {
        let mut worst = 0.0f64;
        for q in [7i64, 16, 40] {
            for p in -3 * q..3 * q {
                let (s, c) = sin_cos_pi_od(p, q);
                let (sd, cd) = sin_cos_pi(p, q);
                assert!(
                    (s.to_dd() - sd).hi.abs() < 1e-31 && (c.to_dd() - cd).hi.abs() < 1e-31,
                    "{p}/{q}"
                );
                let one = od_diff(s * s + c * c, OD::from_f64(1.0));
                let (s2, _) = sin_cos_pi_od(2 * p, q);
                let dbl = od_diff((s * c).mul_f(2.0), s2);
                worst = worst.max(one).max(dbl);
                assert!(one < 1e-126 && dbl < 1e-125, "{p}/{q}: {one:e} {dbl:e}");
            }
        }
        println!("sin_cos_pi_od: worst identity error {worst:.3e}");
    }

    /// The control of `sin_cos_pi_od_agrees`: the series stopped at terms below 1e-90 must fail sin^2 + cos^2 = 1 at
    /// some angle p pi / 160.
    #[test]
    fn sin_cos_pi_od_control_short_series_fails() {
        let worst = (1..40)
            .map(|p| {
                let (s, c) = sin_cos_series_od(PI_OD.mul_f(p as f64) / OD::from_f64(160.0), 1e-90);
                od_diff(s * s + c * c, OD::from_f64(1.0))
            })
            .fold(0.0f64, f64::max);
        println!("control, short series: sin^2 + cos^2 - 1 up to {worst:.3e}");
        assert!(worst > 1e-126, "{worst:e}");
    }

    /// the control: a harmonic of the wrong parity is not differentiated by the folded matrix (cos 2 theta on the
    /// odd grid)
    #[test]
    fn dmat_control_wrong_parity_fails() {
        const B: usize = 1;
        let m = 8 * B;
        let d = dmat::<B>(Parity::Odd);
        let v: Vec<DD> = (0..m).map(|j| sin_cos_grid(Parity::Odd, m, j, 2).1).collect();
        let (h, l) = multidouble::dd::split(&v);
        let c = Fl::<DV, B>::load(&h, &l);
        let s: Vec<DD> = (0..m).map(|j| sin_cos_grid(Parity::Odd, m, j, 2).0).collect();
        let dc = d.apply(&c);
        let err = (0..m)
            .map(|j| (dc.get(j) + s[j].mul_f(2.0)).hi.abs())
            .fold(0.0, nan_max);
        assert!(err > 1e-3, "{err:e}");
    }
}