Resolve issue 5

[terasakisatoshi]

This PR closes #5
Note that this PR breaks test_leg5, which does not appear in the original C++ implementation.

 cargo test
   Compiling xprec v0.2.0 (/Users/atelierarith/work/atelierarith/markus/xprec-rs)
    Finished `test` profile [unoptimized + debuginfo] target(s) in 0.20s
     Running unittests src/lib.rs (target/debug/deps/xprec-229d85e47e0949a4)

running 33 tests
test checks::test::test_class ... ok
test checks::test::test_isclose ... ok
test consts::test::test_consts ... ok
test convert::test::test_from_primitive_edge_cases ... ok
test convert::test::test_from_primitive_f32 ... ok
test convert::test::test_from_primitive_f64 ... ok
test convert::test::test_from_primitive_integers ... ok
test exp::test::test_expm1_kernel ... ok
test funcs::test::test_fract ... ok
test gauss::test::test_leg16 ... ok
test gauss::test::test_leg5 ... FAILED
test gauss::test::test_leg7 ... ok
test exp::test::test_expm1 ... ok
test exp::test::test_log1p ... ok
test hyperbolic::test::test_cosh ... ok
test circular::test::test_circ ... ok
test exp::test::test_exp ... ok
test hyperbolic::test::test_sinh ... ok
test round::test::test_ceil ... ok
test round::test::test_round ... ok
test round::test::test_trunc ... ok
test traits::test::test_traits ... ok
test hyperbolic::test::test_tanh ... ok
test circular::test::test_kernels ... ok
test hyperbolic::test::test_arc ... ok
test arith::test::test_arith_d ... ok
test exp::test::test_log ... ok
test roots::test::test_roots_q ... ok
test roots::test::test_hypot ... ok
test arith::test::test_arith_q ... ok
test arith::test::test_arith_dd ... ok
test arith::test::test_arith_qq ... ok
test arith::test::test_arith_qd ... ok

failures:

---- gauss::test::test_leg5 stdout ----

thread 'gauss::test::test_leg5' panicked at src/gauss.rs:101:13:
assert_abs_diff_eq!(w[i], W5[i], epsilon = 2.0*Df64::EPSILON)

    left  = Compensated { hi: 0.47862867049936647, lo: -2.8616217915820644e-18 }
    right = Compensated { hi: 0.47862867049936647, lo: -2.8616217915821202e-18 }


note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace


failures:
    gauss::test::test_leg5

test result: FAILED. 32 passed; 1 failed; 0 ignored; 0 measured; 0 filtered out; finished in 4.16s

error: test failed, to rerun pass `--lib`

[terasakisatoshi]

I've found that the following C++ code works in the libxprec project, so test_leg5 in Rust should pass.

TEST_CASE("leg-cmp-5", "[gauss]") {
    const static DDouble x_ref[5] = {
        {-0.906179845938664, -2.909730553174891e-17},
        {-0.5384693101056831, 1.6678154894696646e-17},
        {0.0, 0.0},
        {0.5384693101056831, -1.6678154894696646e-17},
        {0.906179845938664, 2.909730553174891e-17}
    };
    const static DDouble w_ref[5] = {
        {0.23692688505618908, 2.6149055638876413e-18},
        {0.47862867049936647, -2.8616217915821202e-18},
        {0.5688888888888889, 4.934324553889585e-19},
        {0.47862867049936647, -2.8616217915821202e-18},
        {0.23692688505618908, 2.6149055638876413e-18},
    };

    std::vector<DDouble> x(5), w(5);
    gauss_legendre(5, x.data(), w.data());
    for (int i = 0; i < 5; ++i) {
        REQUIRE_THAT(x[i], WithinAbs(x_ref[i], std::numeric_limits<DDouble>::epsilon()));
        REQUIRE_THAT(w[i], WithinAbs(w_ref[i], 0.2 * std::numeric_limits<DDouble>::epsilon()));
    }
}

However, I have no idea how to resolve 🫠

[shinaoka]

I think we should carefully check circular.rs, arith.rs before testing gauss.rs since gauss.rs depends on them.

Many tests are still missing in the Rust library:
https://github.com/tuwien-cms/libxprec/blob/mainline/test/circular.cpp

Let us create a dependency graph and check leaves first.

[terasakisatoshi]

The following function should match the original C++ implementation and pass test_leg5, test_leg7 and test_leg16.
Remark: Cursor (AI agent) created this.

pub fn gauss_legendre(x: &mut [Df64], w: &mut [Df64])
{
    let n = x.len();
    gauss_chebyshev_theta(x);
    // Convert theta to x before Newton iteration (matching C++ implementation)
    for i in 0..n {
        x[i] = cos(x[i]);
    }
    // Perform Newton iteration in x space (matching C++ implementation)
    for _iter in 0..10 {
        legendre_x_newton(n as i64, x, w);
    }
    // Compute weights using stored dPn values (matching C++ implementation)
    // Match C++: PowerOfTwo(2.0) * reciprocal((1-x^2) * dPn^2)
    for i in 0..n {
        // Match C++: ExDouble(1.0).add_small(-x[i] * x[i])
        // Use addfast_dq to match C++: add_small(DDouble) -> addfast_dq(double, DDouble)
        // Match C++ calculation order: -x[i] * x[i] = neg(mul(x[i], x[i]))
        let neg_x_sq = neg_q(mul_qq(x[i], x[i]));
        let one_minus_x2 = addfast_dq(1.0, neg_x_sq);
        // w = 2 / ((1-x^2) * dPn^2)
        // Match C++: PowerOfTwo(2.0) * reciprocal(one_minus_x2 * w[i] * w[i])
        // Note: w[i] contains dPn from Newton iteration
        #[allow(non_snake_case)]
        let dPn_sq = square_q(w[i]);
        w[i] = mul_pow2(reciprocal_q(one_minus_x2 * dPn_sq), 2.0);
    }
}

fn gauss_chebyshev_theta(θ: &mut [Df64])
{
    let n = θ.len();
    // Match C++: fact = pi / (1.0 * n)
    let fact = consts::PI / Df64::new(1.0 * n as f64);
    for i in 0..n {
        // Match C++: x[i] = cos((n - i - 0.5) * fact)
        // Store theta first, then convert to x
        θ[i] = (Df64::new((n - i) as f64) - Df64::new(0.5)) * fact;
    }
}

fn legendre_x_newton(n: i64, x: &mut [Df64], w: &mut [Df64])
{
    // Newton iteration in x space (matching C++ implementation)
    // Match C++: leg_deriv(n, x[i], Pn, dPn); dx = -Pn / dPn; x[i] += dx; w[i] = dPn;
    for i in 0..x.len() {
        #[allow(non_snake_case)]
        let (pn, dPn) = leg_deriv(n, x[i]);
        // Match C++: dx = -Pn / dPn
        // Use explicit division to match C++ order
        let neg_pn = neg_q(pn);
        let dx = div_qq(neg_pn, dPn);
        // Match C++: x[i] += dx
        // Use addfast_qq to match C++ '+' operator
        x[i] = addfast_qq(x[i], dx);
        // Store dPn for weight computation (matching C++ implementation)
        w[i] = dPn;
    }
}

// Compute Legendre polynomial and its derivative using recursion (matching C++ leg_deriv)
fn leg_deriv(n: i64, x: Df64) -> (Df64, Df64)
{
    assert!(n >= 1);

    let mut pn_1 = Df64::ONE;
    let mut dpn_1 = Df64::ZERO;
    let mut pn = x;
    let mut dpn = Df64::ONE;

    for k in 1..n {
        // Match C++: ((2 * n + 1.0) * x * Pn - n * Pn_1) / (n + 1.0)
        let two_k_plus_one = Df64::new((2 * k + 1) as f64);
        let k_dd = Df64::new(k as f64);
        let k_plus_one = (k + 1) as f64; // Match C++: (n + 1.0) is double
        // Match C++ calculation order exactly
        // Match C++: operator- uses add_qq(x, -y)
        let two_k_plus_one_x_pn = two_k_plus_one * x * pn;
        let k_dd_pn_1 = k_dd * pn_1;
        let pnext = div_qd(add_qq(two_k_plus_one_x_pn, neg_q(k_dd_pn_1)), k_plus_one);
        // Match C++: ((2 * n + 1.0) * (x * dPn + Pn) - n * dPn_1) / (n + 1.0)
        // Compute x * dpn first, then add pn to match C++ order
        // Use add_qq to match C++ '+' operator (xprec_add_qq)
        let x_dpn = x * dpn;
        let x_dpn_plus_pn = add_qq(x_dpn, pn);
        let two_k_plus_one_x_dpn_plus_pn = two_k_plus_one * x_dpn_plus_pn;
        let k_dd_dpn_1 = k_dd * dpn_1;
        // Match C++: operator- uses add_qq(x, -y)
        let dpnext = div_qd(add_qq(two_k_plus_one_x_dpn_plus_pn, neg_q(k_dd_dpn_1)), k_plus_one);

        // shift terms by one
        pn_1 = pn;
        pn = pnext;
        dpn_1 = dpn;
        dpn = dpnext;
    }
    return (pn, dpn);
}

[shinaoka]

@mwallerb

Do you have any idea? Why did you switch to a different algorithm?

xprec-rs/src/gauss.rs

Lines 35 to 48 in e9ce611

	fn legendre_theta_newton(n: i64, θ: &mut [Df64], w: &mut [Df64])
	{
	// Newton iteration for theta rather than x
	// SIAM J. SCI. COMPUT., Vol. 35, No. 2, p. A652
	#[allow(non_snake_case)]
	for i in 0..θ.len() {
	let (s, c) = sincos(θ[i]);
	let (pn_1, pn) = plx(n, c);
	let pn_θ = (n as f64) * (c * pn - pn_1) / s;
	let Δθ = pn / pn_θ;
	θ[i] -= Δθ;
	w[i] = 2.0 / square_q(pn_θ);
	}
	}

[mwallerb]

Thanks a lot for the added unit tests!

I suspect the test failures are due to numerical noise. Instead of changing the algorithm etc., let's just relax the tolerances slightly.

[mwallerb]

@mwallerb
Do you have any idea? Why did you switch to a different algorithm?

The other algorithm has better error properties in general, especially for higher Gauss order. The tolerances in the original implementation were tighter than what theory would suggest, so (depending on the exact algorithm) tests sometimes fail or they don't.

If we want to squeeze out the very last bits of precision, then we will eventually have to switch to a sextuple-precision algorithm. We could also add a lookup table with the exact values for some low-order quadratures we use, say 8 and 16. I doubt that it is really necessary though.

[mwallerb]

Great job!~