Ensure that all quadrature rules store their nodes in sorted order.

CHANGELOG.md

@@ -1,5 +1,13 @@
 This document contains all changes to the crate since version 0.1.8.
 
+# 0.3.0 (unreleased)
+
+- Made all quadrature rules store their nodes sorted in ascending order.
+ This means that all functions that return some view of the nodes (and weights) now return them in sorted order.
+ This affects the Gauss-Legendre, Gauss-Hermite and Gauss-Chebyshev rules.
+ The affected functions are `QuadratureRule::iter()`, `QuadratureRule::into_iter()`, `QuadratureRule::nodes()`, `QuadratureRule::weights()`, `QuadratureRule::as_node_weight_pairs()` and `QuadratureRule::into_node_weight_pairs()`.
+
+
 # 0.2.2
 
 - Add Gauss-Chebyshev quadrature of the first and second kinds.

src/chebyshev/mod.rs

@@ -8,7 +8,9 @@ use core::{f64::consts::PI, fmt};
 use std::backtrace::Backtrace;
 
 #[cfg(feature = "rayon")]
-use rayon::iter::{IntoParallelIterator, IntoParallelRefIterator, ParallelIterator};
+use rayon::iter::{
+    IndexedParallelIterator, IntoParallelIterator, IntoParallelRefIterator, ParallelIterator,
+};
 
 /// A Gauss-Chebyshev quadrature scheme of the first kind.
 ///
@@ -47,6 +49,7 @@ impl GaussChebyshevFirstKind {
 
         Ok(Self {
             node_weight_pairs: (1..degree + 1)
+                .rev()
                 .map(|i| ((PI * (2.0 * (i as f64) - 1.0) / (2.0 * n)).cos(), PI / n))
                 .collect(),
         })
@@ -72,6 +75,7 @@ impl GaussChebyshevFirstKind {
         Ok(Self {
             node_weight_pairs: (1..degree + 1)
                 .into_par_iter()
+                .rev()
                 .map(|i| ((PI * (2.0 * (i as f64) - 1.0) / (2.0 * n)).cos(), PI / n))
                 .collect(),
         })
@@ -156,6 +160,7 @@ impl GaussChebyshevSecondKind {
 
         Ok(Self {
             node_weight_pairs: (1..degree + 1)
+                .rev()
                 .map(|i| {
                     let over_n_plus_1 = 1.0 / (n + 1.0);
                     let sin_val = (PI * i as f64 * over_n_plus_1).sin();
@@ -180,6 +185,7 @@ impl GaussChebyshevSecondKind {
         Ok(Self {
             node_weight_pairs: (1..degree + 1)
                 .into_par_iter()
+                .rev()
                 .map(|i| {
                     let over_n_plus_1 = 1.0 / (n + 1.0);
                     let sin_val = (PI * i as f64 * over_n_plus_1).sin();
@@ -281,15 +287,36 @@ mod test {
         assert!(GaussChebyshevSecondKind::new(1).is_err());
     }
 
+    #[test]
+    fn check_sorted() {
+        for deg in (2..100).step_by(10) {
+            let rule1 = GaussChebyshevFirstKind::new(deg).unwrap();
+            assert!(rule1.as_node_weight_pairs().is_sorted());
+            let rule2 = GaussChebyshevSecondKind::new(deg).unwrap();
+            assert!(rule2.as_node_weight_pairs().is_sorted());
+        }
+    }
+
+    #[cfg(feature = "rayon")]
+    #[test]
+    fn check_par_sorted() {
+        for deg in (2..100).step_by(10) {
+            let rule1 = GaussChebyshevFirstKind::par_new(deg).unwrap();
+            assert!(rule1.as_node_weight_pairs().is_sorted());
+            let rule2 = GaussChebyshevSecondKind::par_new(deg).unwrap();
+            assert!(rule2.as_node_weight_pairs().is_sorted());
+        }
+    }
+
     #[test]
     fn check_chebyshev_1st_deg_5() {
         // Source: https://mathworld.wolfram.com/Chebyshev-GaussQuadrature.html
         let ans = [
-            (0.5 * (0.5 * (5.0 + f64::sqrt(5.0))).sqrt(), PI / 5.0),
-            (0.5 * (0.5 * (5.0 - f64::sqrt(5.0))).sqrt(), PI / 5.0),
-            (0.0, PI / 5.0),
-            (-0.5 * (0.5 * (5.0 - f64::sqrt(5.0))).sqrt(), PI / 5.0),
             (-0.5 * (0.5 * (5.0 + f64::sqrt(5.0))).sqrt(), PI / 5.0),
+            (-0.5 * (0.5 * (5.0 - f64::sqrt(5.0))).sqrt(), PI / 5.0),
+            (0.0, PI / 5.0),
+            (0.5 * (0.5 * (5.0 - f64::sqrt(5.0))).sqrt(), PI / 5.0),
+            (0.5 * (0.5 * (5.0 + f64::sqrt(5.0))).sqrt(), PI / 5.0),
         ];
 
         let rule = GaussChebyshevFirstKind::new(5).unwrap();
@@ -311,7 +338,7 @@ mod test {
 
         for (i, (x, w)) in rule.into_iter().enumerate() {
             // Source: https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature
-            let ii = i as f64 + 1.0;
+            let ii = deg - i as f64;
             let x_should = (ii * PI / (deg + 1.0)).cos();
             let w_should = PI / (deg + 1.0) * (ii * PI / (deg + 1.0)).sin().powi(2);
 

src/hermite/mod.rs

@@ -89,22 +89,26 @@ impl GaussHermite {
         // calculate eigenvalues & vectors
         let eigen = companion_matrix.symmetric_eigen();
 
-        // zip together the iterator over nodes with the one over weights and return as Vec<(f64, f64)>
-        Ok(GaussHermite {
-            node_weight_pairs: eigen
-                .eigenvalues
-                .iter()
-                .copied()
-                .zip(
-                    eigen
-                        .eigenvectors
-                        .row(0)
-                        .map(|x| x * x * PI.sqrt())
-                        .iter()
-                        .copied(),
-                )
-                .collect(),
-        })
+        // zip together the iterator over nodes with the one over weights and collect into a Vec<(f64, f64)>
+        let mut node_weight_pairs: Vec<(Node, Weight)> = eigen
+            .eigenvalues
+            .iter()
+            .copied()
+            .zip(
+                eigen
+                    .eigenvectors
+                    .row(0)
+                    .map(|x| x * x * PI.sqrt())
+                    .iter()
+                    .copied(),
+            )
+            .collect();
+
+        // sort the nodes and weights by the nodes
+        node_weight_pairs
+            .sort_unstable_by(|(node1, _), (node2, _)| node1.partial_cmp(node2).unwrap());
+
+        Ok(GaussHermite { node_weight_pairs })
     }
 
     /// Perform quadrature of e^(-x^2) * `integrand`(x) over the domain (-∞, ∞).
@@ -174,10 +178,18 @@ mod tests {
 
     use super::*;
 
+    #[test]
+    fn check_sorted() {
+        for deg in (2..100).step_by(10) {
+            let rule = GaussHermite::new(deg).unwrap();
+            assert!(rule.as_node_weight_pairs().is_sorted());
+        }
+    }
+
     #[test]
     fn golub_welsch_3() {
         let (x, w): (Vec<_>, Vec<_>) = GaussHermite::new(3).unwrap().into_iter().unzip();
-        let x_should = [1.224_744_871_391_589, 0.0, -1.224_744_871_391_589];
+        let x_should = [-1.224_744_871_391_589, 0.0, 1.224_744_871_391_589];
         let w_should = [
             0.295_408_975_150_919_35,
             1.181_635_900_603_677_4,
@@ -242,7 +254,7 @@ mod tests {
     fn integrate_one() {
         let quad = GaussHermite::new(5).unwrap();
         let integral = quad.integrate(|_x| 1.0);
-        assert_abs_diff_eq!(integral, PI.sqrt(), epsilon = 1e-15);
+        assert_abs_diff_eq!(integral, PI.sqrt(), epsilon = 1e-14);
     }
 
     #[cfg(feature = "rayon")]

src/jacobi/mod.rs

@@ -145,7 +145,8 @@ impl GaussJacobi {
             )
             .collect();
 
-        node_weight_pairs.sort_unstable_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
+        node_weight_pairs
+            .sort_unstable_by(|(node1, _), (node2, _)| node1.partial_cmp(node2).unwrap());
 
         // TO FIX: implement correction
         // eigenvalue algorithm has problem to get the zero eigenvalue for odd degrees
@@ -274,13 +275,8 @@ impl std::error::Error for GaussJacobiError {}
 /// Gauss-Legendre quadrature is equivalent to Gauss-Jacobi quadrature with `alpha` = `beta` = 0.
 impl From<GaussLegendre> for GaussJacobi {
     fn from(value: GaussLegendre) -> Self {
-        let mut node_weight_pairs = value.into_node_weight_pairs();
-        // Gauss-Legendre nodes are generated in reverse sorted order.
-        // This corrects for that since Gauss-Jacobi nodes are currently always sorted
-        // in ascending order.
-        node_weight_pairs.reverse();
         Self {
-            node_weight_pairs,
+            node_weight_pairs: value.into_node_weight_pairs(),
             alpha: 0.0,
             beta: 0.0,
         }
@@ -364,6 +360,26 @@ mod tests {
 
     use super::*;
 
+    #[test]
+    fn check_sort() {
+        // This contains values such that all the possible combinations of two of them
+        // contains the combinations
+        // (0, 0), which is Gauss-Legendre
+        // (-0.5, -0.5), which is Gauss-Chebyshev of the first kind
+        // (0.5, 0.5), which is Gauss-Chebyshev of the second kind
+        // and all the other combinations of the values are Gauss-Jacobi
+        const PARAMS: [f64; 4] = [-0.5, -0.25, 0.0, 0.5];
+
+        for deg in (2..100).step_by(20) {
+            for alpha in PARAMS {
+                for beta in PARAMS {
+                    let rule = GaussJacobi::new(deg, alpha, beta).unwrap();
+                    assert!(rule.as_node_weight_pairs().is_sorted());
+                }
+            }
+        }
+    }
+
     #[test]
     fn sanity_check_chebyshev_delegation() {
         const DEG: usize = 200;
@@ -384,10 +400,7 @@ mod tests {
         let jrule = GaussJacobi::new(DEG, 0.0, 0.0).unwrap();
         let lrule = GaussLegendre::new(DEG).unwrap();
 
-        assert_eq!(
-            jrule.as_node_weight_pairs(),
-            lrule.into_iter().rev().collect::<Vec<_>>(),
-        );
+        assert_eq!(jrule.as_node_weight_pairs(), lrule.as_node_weight_pairs(),);
     }
 
     #[test]

src/laguerre/mod.rs

@@ -112,7 +112,9 @@ impl GaussLaguerre {
                     .copied(),
             )
             .collect();
-        node_weight_pairs.sort_unstable_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
+
+        node_weight_pairs
+            .sort_unstable_by(|(node1, _), (node2, _)| node1.partial_cmp(node2).unwrap());
 
         Ok(GaussLaguerre {
             node_weight_pairs,
@@ -237,6 +239,16 @@ mod tests {
     use approx::assert_abs_diff_eq;
     use core::f64::consts::PI;
 
+    #[test]
+    fn check_sorted() {
+        for deg in (2..100).step_by(10) {
+            for alpha in [-0.9, -0.5, 0.0, 0.5] {
+                let rule = GaussLaguerre::new(deg, alpha).unwrap();
+                assert!(rule.as_node_weight_pairs().is_sorted());
+            }
+        }
+    }
+
     #[test]
     fn golub_welsch_2_alpha_5() {
         let (x, w): (Vec<_>, Vec<_>) = GaussLaguerre::new(2, 5.0).unwrap().into_iter().unzip();

src/legendre/mod.rs

@@ -24,7 +24,9 @@
 //! ```
 
 #[cfg(feature = "rayon")]
-use rayon::prelude::{IntoParallelIterator, IntoParallelRefIterator, ParallelIterator};
+use rayon::prelude::{
+    IndexedParallelIterator, IntoParallelIterator, IntoParallelRefIterator, ParallelIterator,
+};
 
 mod bogaert;
 
@@ -41,6 +43,7 @@ use std::backtrace::Backtrace;
 /// # Examples
 ///
 /// Basic usage:
+///
 /// ```
 /// # use gauss_quad::legendre::{GaussLegendre, GaussLegendreError};
 /// # use approx::assert_abs_diff_eq;
@@ -53,16 +56,18 @@ use std::backtrace::Backtrace;
 /// assert_abs_diff_eq!(integral, 0.0);
 /// # Ok::<(), GaussLegendreError>(())
 /// ```
+///
 /// The nodes and weights are computed in O(n) time,
 /// so large quadrature rules are feasible:
+///
 /// ```
 /// # use gauss_quad::legendre::{GaussLegendre, GaussLegendreError};
 /// # use approx::assert_abs_diff_eq;
 /// let quad = GaussLegendre::new(1_000_000)?;
 ///
 /// let integral = quad.integrate(-3.0, 3.0, |x| x.sin());
 ///
-/// assert_abs_diff_eq!(integral, 0.0);
+/// assert_abs_diff_eq!(integral, 0.0, epsilon = 8e-16);
 /// # Ok::<(), GaussLegendreError>(())
 /// ```
 #[derive(Debug, Clone, PartialEq)]
@@ -89,6 +94,7 @@ impl GaussLegendre {
 
         Ok(Self {
             node_weight_pairs: (1..deg + 1)
+                .rev()
                 .map(|k: usize| NodeWeightPair::new(deg, k).into_tuple())
                 .collect(),
         })
@@ -108,6 +114,7 @@ impl GaussLegendre {
         Ok(Self {
             node_weight_pairs: (1..deg + 1)
                 .into_par_iter()
+                .rev()
                 .map(|k| NodeWeightPair::new(deg, k).into_tuple())
                 .collect(),
         })
@@ -212,11 +219,28 @@ mod tests {
 
     use super::*;
 
+    #[test]
+    fn check_sorted() {
+        for deg in (2..200).step_by(10) {
+            let rule = GaussLegendre::new(deg).unwrap();
+            assert!(rule.as_node_weight_pairs().is_sorted());
+        }
+    }
+
+    #[cfg(feature = "rayon")]
+    #[test]
+    fn check_par_sorted_order() {
+        for deg in (2..100).step_by(10) {
+            let rule = GaussLegendre::par_new(deg).unwrap();
+            assert!(rule.as_node_weight_pairs().is_sorted());
+        }
+    }
+
     #[test]
     fn check_degree_3() {
         let (x, w): (Vec<_>, Vec<_>) = GaussLegendre::new(3).unwrap().into_iter().unzip();
 
-        let x_should = [0.7745966692414834, 0.0000000000000000, -0.7745966692414834];
+        let x_should = [-0.7745966692414834, 0.0000000000000000, 0.7745966692414834];
         let w_should = [0.5555555555555556, 0.8888888888888888, 0.5555555555555556];
         for (i, x_val) in x_should.iter().enumerate() {
             assert_abs_diff_eq!(*x_val, x[i]);
@@ -240,8 +264,8 @@ mod tests {
         #[allow(clippy::excessive_precision)]
         let w_should = [0.0244461801962625182113259,0.0244315690978500450548486,0.0244023556338495820932980,0.0243585572646906258532685,0.0243002001679718653234426,0.0242273192228152481200933,0.0241399579890192849977167,0.0240381686810240526375873,0.0239220121367034556724504,0.0237915577810034006387807,0.0236468835844476151436514,0.0234880760165359131530253,0.0233152299940627601224157,0.0231284488243870278792979,0.0229278441436868469204110,0.0227135358502364613097126,0.0224856520327449668718246,0.0222443288937997651046291,0.0219897106684604914341221,0.0217219495380520753752610,0.0214412055392084601371119,0.0211476464682213485370195,0.0208414477807511491135839,0.0205227924869600694322850,0.0201918710421300411806732,0.0198488812328308622199444,0.0194940280587066028230219,0.0191275236099509454865185,0.0187495869405447086509195,0.0183604439373313432212893,0.0179603271850086859401969,0.0175494758271177046487069,0.0171281354231113768306810,0.0166965578015892045890915,0.0162550009097851870516575,0.0158037286593993468589656,0.0153430107688651440859909,0.0148731226021473142523855,0.0143943450041668461768239,0.0139069641329519852442880,0.0134112712886163323144890,0.0129075627392673472204428,0.0123961395439509229688217,0.0118773073727402795758911,0.0113513763240804166932817,0.0108186607395030762476596,0.0102794790158321571332153,0.0097341534150068058635483,0.0091830098716608743344787,0.0086263777986167497049788,0.0080645898904860579729286,0.0074979819256347286876720,0.0069268925668988135634267,0.0063516631617071887872143,0.0057726375428656985893346,0.0051901618326763302050708,0.0046045842567029551182905,0.0040162549837386423131943,0.0034255260409102157743378,0.0028327514714579910952857,0.0022382884309626187436221,0.0016425030186690295387909,0.0010458126793403487793129,0.0004493809602920903763943];
 
-        for (i, x_val) in x_should.iter().rev().enumerate() {
-            assert_abs_diff_eq!(*x_val, x[i], epsilon = 0.000_000_1);
+        for (i, &x_val) in x_should.iter().rev().enumerate() {
+            assert_abs_diff_eq!(-x_val, x[i], epsilon = 0.000_000_1);
         }
         for (i, w_val) in w_should.iter().rev().enumerate() {
             assert_abs_diff_eq!(*w_val, w[i], epsilon = 0.000_000_1);