add hermite and laguerre tables

src/numerical_derivative/derivator.rs

@@ -3,10 +3,20 @@ use crate::utils::error_codes::ErrorCode;
 
 pub trait DerivatorSingleVariable: Default + Clone + Copy
 {
+    ///generic n-th derivative for a single variable function
     fn get<T: ComplexFloat>(&self, order: usize, func: &dyn Fn(T) -> T, point: T) -> Result<T, ErrorCode>;
 }
 
 pub trait DerivatorMultiVariable: Default + Clone + Copy
 {
+    ///specialized version for a single partial derivative
+    /// Specialized versions for most-used cases make for a smoother experience with API
+    fn get_single_partial<T: ComplexFloat, const NUM_VARS: usize>(&self, func: &dyn Fn(&[T; NUM_VARS]) -> T, idx_to_derivate: usize, point: &[T; NUM_VARS]) -> Result<T, ErrorCode>;
+
+    ///specialized version for a double partial derivative
+    ///Specialized versions for most-used cases make for a smoother experience with API
+    fn get_double_partial<T: ComplexFloat, const NUM_VARS: usize>(&self, func: &dyn Fn(&[T; NUM_VARS]) -> T, idx_to_derivate: &[usize; 2], point: &[T; NUM_VARS]) -> Result<T, ErrorCode>;
+
+    ///generic n-th derivative for a multivariable function
     fn get<T: ComplexFloat, const NUM_VARS: usize, const NUM_ORDER: usize>(&self, order: usize, func: &dyn Fn(&[T; NUM_VARS]) -> T, idx_to_derivate: &[usize; NUM_ORDER], point: &[T; NUM_VARS]) -> Result<T, ErrorCode>;
 }

src/numerical_derivative/finite_difference.rs

@@ -287,4 +287,34 @@ impl DerivatorMultiVariable for MultiVariableSolver
             mode::FiniteDifferenceMode::Central => return Ok(self.get_central_difference_multi_variable(order, func, idx_to_derivate, point)) 
         }    
     }
+
+    fn get_single_partial<T: ComplexFloat, const NUM_VARS: usize>(&self, func: &dyn Fn(&[T; NUM_VARS]) -> T, idx_to_derivate: usize, point: &[T; NUM_VARS]) -> Result<T, ErrorCode> 
+    {
+        if self.step_size == 0.0
+        {
+            return Err(ErrorCode::NumberOfStepsCannotBeZero);
+        }
+
+        match self.method
+        {
+            mode::FiniteDifferenceMode::Forward => return Ok(self.get_forward_difference_multi_variable(1, func, &[idx_to_derivate], point)),
+            mode::FiniteDifferenceMode::Backward => return Ok(self.get_backward_difference_multi_variable(1, func, &[idx_to_derivate], point)),
+            mode::FiniteDifferenceMode::Central => return Ok(self.get_central_difference_multi_variable(1, func, &[idx_to_derivate], point)) 
+        }    
+    }
+
+    fn get_double_partial<T: ComplexFloat, const NUM_VARS: usize>(&self, func: &dyn Fn(&[T; NUM_VARS]) -> T, idx_to_derivate: &[usize; 2], point: &[T; NUM_VARS]) -> Result<T, ErrorCode> 
+    {
+        if self.step_size == 0.0
+        {
+            return Err(ErrorCode::NumberOfStepsCannotBeZero);
+        }
+
+        match self.method
+        {
+            mode::FiniteDifferenceMode::Forward => return Ok(self.get_forward_difference_multi_variable(2, func, idx_to_derivate, point)),
+            mode::FiniteDifferenceMode::Backward => return Ok(self.get_backward_difference_multi_variable(2, func, idx_to_derivate, point)),
+            mode::FiniteDifferenceMode::Central => return Ok(self.get_central_difference_multi_variable(2, func, idx_to_derivate, point)) 
+        }  
+    }
 }

src/numerical_derivative/hessian.rs

@@ -71,7 +71,7 @@ impl<D: DerivatorMultiVariable> Hessian<D>
             {
                 if result[row_index][col_index].is_nan()
                 {
-                    result[row_index][col_index] = self.derivator.get(2, function, &[row_index, col_index], vector_of_points)?;
+                    result[row_index][col_index] = self.derivator.get_double_partial(function, &[row_index, col_index], vector_of_points)?;
 
                     result[col_index][row_index] = result[row_index][col_index]; //exploit the fact that a hessian is a symmetric matrix
                 }

src/numerical_derivative/jacobian.rs

@@ -94,7 +94,7 @@ impl<D: DerivatorMultiVariable> Jacobian<D>
         {
             for col_index in 0..NUM_VARS
             {
-                result[row_index][col_index] = self.derivator.get(1, &function_matrix[row_index], &[col_index], vector_of_points)?;
+                result[row_index][col_index] = self.derivator.get_single_partial(&function_matrix[row_index], col_index, vector_of_points)?;
             }
         }
 
@@ -132,7 +132,7 @@ impl<D: DerivatorMultiVariable> Jacobian<D>
             let mut cur_row: Vec<T> = Vec::new();
             for col_index in 0..NUM_VARS
             {
-                cur_row.push(self.derivator.get(1,&function_matrix[row_index], &[col_index], vector_of_points)?);
+                cur_row.push(self.derivator.get_single_partial(&function_matrix[row_index], col_index, vector_of_points)?);
             }
 
             result.push(cur_row);

src/numerical_derivative/test.rs

@@ -75,7 +75,7 @@ fn test_single_derivative_partial_1()
     derivator.set_method(FiniteDifferenceMode::Central);
 
     //partial derivate for (x, y) = (1.0, 3.0), partial derivative for x is known to be 6*x + 2*y
-    let val = derivator.get(1, &func, &[0], &point).unwrap();
+    let val = derivator.get_single_partial(&func, 0, &point).unwrap();
     assert!(f64::abs(val - 12.0) < 0.000001);
 
     //partial derivate for (x, y) = (1.0, 3.0), partial derivative for y is known to be 2.0*x
@@ -521,4 +521,4 @@ fn test_hessian_1()
             assert!(f64::abs(result[i][j] - expected_result[i][j]) < 0.01);
         }
     }
-}
+}

src/numerical_integration/gaussian_integration.rs

@@ -4,13 +4,13 @@ use num_complex::ComplexFloat;
 use crate::utils::gl_table as gl_table;
 use crate::numerical_integration::integrator::Integrator;
 
-pub struct Gaussian
+pub struct GaussianQuadrature
 {
     order: usize,
     method_type: mode::GaussianMethod
 }
 
-impl Gaussian
+impl GaussianQuadrature
 {
     pub fn get_order(&self) -> usize
     {
@@ -34,7 +34,7 @@ impl Gaussian
 
     pub fn with_parameters(order: usize, integration_method: mode::GaussianMethod) -> Self 
     {
-        Gaussian
+        GaussianQuadrature
         {
             order: order,
             method_type: integration_method
@@ -43,9 +43,9 @@ impl Gaussian
 
     fn check_for_errors<T: ComplexFloat>(&self, integration_limit: &[T; 2]) -> Result<(), ErrorCode> 
     {
-        if !(2..=gl_table::MAX_GL_ORDER).contains(&self.order)
+        if !(1..=gl_table::MAX_GL_ORDER).contains(&self.order)
         {
-            return Err(ErrorCode::GaussLegendreOrderOutOfRange);
+            return Err(ErrorCode::GaussianQuadratureOrderOutOfRange);
         }
 
         if integration_limit[0].abs() >= integration_limit[1].abs()
@@ -100,7 +100,7 @@ impl Gaussian
     }
 }
 
-impl Integrator for Gaussian
+impl Integrator for GaussianQuadrature
 {
     fn get_single_total<T: ComplexFloat, const NUM_VARS: usize>(&self, func: &dyn Fn(&[T; NUM_VARS]) -> T, integration_limit: &[T; 2]) -> Result<T, ErrorCode>
     {

src/numerical_integration/test.rs

@@ -1,8 +1,10 @@
+use std::println;
+
 use crate::numerical_integration::mode::*;
 use crate::utils::error_codes::ErrorCode;
 use crate::numerical_integration::integrator::Integrator;
 use crate::numerical_integration::iterative_integration::Iterative;
-use crate::numerical_integration::gaussian_integration::Gaussian;
+use crate::numerical_integration::gaussian_integration::GaussianQuadrature;
 
 #[test]
 fn test_booles_integration_1()
@@ -103,11 +105,11 @@ fn test_gauss_legendre_quadrature_integration_1()
 
     let integration_limit = [0.0, 2.0];
 
-    let integrator = Gaussian::with_parameters(2, GaussianMethod::GaussLegendre);
+    let integrator = GaussianQuadrature::with_parameters(4, GaussianMethod::GaussLegendre);
 
     //simple integration for x, known to be x^4 - x^3, expect a value of ~8.00
     let val = integrator.get_single_total(&func, &integration_limit).unwrap();
-    assert!(f64::abs(val - 8.0) < 0.00001);
+    assert!(f64::abs(val - 8.0) < 1e-14);
 }
 
 #[test] 
@@ -122,25 +124,25 @@ fn test_gauss_legendre_quadrature_integration_2()
     let integration_limit = [0.0, 1.0];
     let point = [1.0, 2.0, 3.0];
 
-    let integrator = Gaussian::with_parameters(2, GaussianMethod::GaussLegendre);
+    let integrator = GaussianQuadrature::with_parameters(2, GaussianMethod::GaussLegendre);
 
     //partial integration for x, known to be x*x + x*y*z, expect a value of ~7.00
     let val = integrator.get_single_partial(&func, 0, &integration_limit, &point).unwrap();
-    assert!(f64::abs(val - 7.0) < 0.00001);
+    assert!(f64::abs(val - 7.0) < 1e-14);
 
 
     let integration_limit = [0.0, 2.0];
 
     //partial integration for y, known to be 2.0*x*y + y*y*z/2.0, expect a value of ~10.00 
     let val = integrator.get_single_partial(&func, 1, &integration_limit, &point).unwrap();
-    assert!(f64::abs(val - 10.0) < 0.00001);
+    assert!(f64::abs(val - 10.0) < 1e-14);
 
 
     let integration_limit = [0.0, 3.0];
 
     //partial integration for z, known to be 2.0*x*z + y*z*z/2.0, expect a value of ~15.0 
     let val = integrator.get_single_partial(&func, 2, &integration_limit, &point).unwrap();
-    assert!(f64::abs(val - 15.0) < 0.00001);
+    assert!(f64::abs(val - 15.0) < 1e-14);
 }
 
 #[test]
@@ -153,11 +155,11 @@ fn test_gauss_legendre_quadrature_integration_3()
     };
 
     let integration_limits = [[0.0, 2.0], [0.0, 2.0]];
-    let integrator = Gaussian::with_parameters(2, GaussianMethod::GaussLegendre);
+    let integrator = GaussianQuadrature::with_parameters(2, GaussianMethod::GaussLegendre);
 
     //simple double integration for 6*x, expect a value of ~24.00
     let val = integrator.get_double_total(&func, &integration_limits).unwrap();
-    assert!(f64::abs(val - 24.0) < 0.00001);
+    assert!(f64::abs(val - 24.0) < 1e-14);
 }
 
 #[test]
@@ -384,11 +386,11 @@ fn test_error_checking_3()
 
     let integration_limit = [0.0, 2.0];
 
-    //Gauss Legendre not valid for n < 2
-    let integrator = Gaussian::with_parameters(1, GaussianMethod::GaussLegendre);
+    //Gauss Legendre not valid for n < 1
+    let integrator = GaussianQuadrature::with_parameters(0, GaussianMethod::GaussLegendre);
     let result = integrator.get_single_total(&func, &integration_limit);
     assert!(result.is_err());
-    assert!(result.unwrap_err() == ErrorCode::GaussLegendreOrderOutOfRange);
+    assert!(result.unwrap_err() == ErrorCode::GaussianQuadratureOrderOutOfRange);
 }
 
 #[test]
@@ -403,8 +405,8 @@ fn test_error_checking_4()
     let integration_limit = [0.0, 2.0];
 
     //Gauss Legendre not valid for n > 20
-    let integrator = Gaussian::with_parameters(21, GaussianMethod::GaussLegendre);
+    let integrator = GaussianQuadrature::with_parameters(21, GaussianMethod::GaussLegendre);
     let result = integrator.get_single_total(&func, &integration_limit);
     assert!(result.is_err());
-    assert!(result.unwrap_err() == ErrorCode::GaussLegendreOrderOutOfRange);
+    assert!(result.unwrap_err() == ErrorCode::GaussianQuadratureOrderOutOfRange);
 }

src/utils/error_codes.rs

@@ -26,5 +26,5 @@ pub enum ErrorCode
 
     //Can be returned by single_integration and double_integration if using the Gauss Legendre integration method
     //Returned if requested order of integration is < 2 or > 15
-    GaussLegendreOrderOutOfRange
+    GaussianQuadratureOrderOutOfRange
 }

src/utils/gauss_laguerre_table.rs

@@ -0,0 +1,83 @@
+pub const MAX_GH_ORDER: usize = 15;
+use crate::utils::error_codes::ErrorCode;
+
+pub fn get_gauss_laguerre_weights_and_abscissae(order: usize, index: usize) -> Result<(f64, f64), ErrorCode>
+{
+    let ref_abs: f64 = match order 
+    {
+        1 => LAGUERRE_ABSCISSA_1[index],
+        2 => LAGUERRE_ABSCISSA_2[index],
+        3 => LAGUERRE_ABSCISSA_3[index],
+        4 => LAGUERRE_ABSCISSA_4[index],
+        5 => LAGUERRE_ABSCISSA_5[index],
+        6 => LAGUERRE_ABSCISSA_6[index],
+        7 => LAGUERRE_ABSCISSA_7[index],
+        8 => LAGUERRE_ABSCISSA_8[index],
+        9 => LAGUERRE_ABSCISSA_9[index],
+        10 => LAGUERRE_ABSCISSA_10[index],
+        11 => LAGUERRE_ABSCISSA_11[index],
+        12 => LAGUERRE_ABSCISSA_12[index],
+        13 => LAGUERRE_ABSCISSA_13[index],
+        14 => LAGUERRE_ABSCISSA_14[index],
+        15 => LAGUERRE_ABSCISSA_15[index],
+        _ => return Err(ErrorCode::GaussianQuadratureOrderOutOfRange),
+    };
+
+    let ref_weight: f64 = match order 
+    {
+        1 => LAGUERRE_WEIGHT_1[index],
+        2 => LAGUERRE_WEIGHT_2[index],
+        3 => LAGUERRE_WEIGHT_3[index],
+        4 => LAGUERRE_WEIGHT_4[index],
+        5 => LAGUERRE_WEIGHT_5[index],
+        6 => LAGUERRE_WEIGHT_6[index],
+        7 => LAGUERRE_WEIGHT_7[index],
+        8 => LAGUERRE_WEIGHT_8[index],
+        9 => LAGUERRE_WEIGHT_9[index],
+        10 => LAGUERRE_WEIGHT_10[index],
+        11 => LAGUERRE_WEIGHT_11[index],
+        12 => LAGUERRE_WEIGHT_12[index],
+        13 => LAGUERRE_WEIGHT_13[index],
+        14 => LAGUERRE_WEIGHT_14[index],
+        15 => LAGUERRE_WEIGHT_15[index],
+        _ => return Err(ErrorCode::GaussianQuadratureOrderOutOfRange),
+    };
+
+    return Ok((ref_abs, ref_weight));
+}
+
+// =============================================================================
+// Table generated using Python's numpy.polynomial.laguerre
+// =============================================================================
+const LAGUERRE_ABSCISSA_1: [f64; 1] = [1.0];
+const LAGUERRE_ABSCISSA_2: [f64; 2] = [0.58578644, 3.41421356];
+const LAGUERRE_ABSCISSA_3: [f64; 3] = [0.41577456, 2.29428036, 6.28994508];
+const LAGUERRE_ABSCISSA_4: [f64; 4] = [0.32254769, 1.7457611, 4.5366203, 9.39507091];
+const LAGUERRE_ABSCISSA_5: [f64; 5] = [0.26356032, 1.41340306, 3.59642577, 7.08581001, 12.64080084];
+const LAGUERRE_ABSCISSA_6: [f64; 6] = [0.2228466, 1.1889321, 2.99273633, 5.77514357, 9.83746742, 15.98287398];
+const LAGUERRE_ABSCISSA_7: [f64; 7] = [0.19304368, 1.0266649, 2.56787674, 4.90035308, 8.18215344, 12.73418029, 19.39572786];
+const LAGUERRE_ABSCISSA_8: [f64; 8] = [0.17027963, 0.90370178, 2.25108663, 4.26670017, 7.0459054, 10.75851601, 15.74067864, 22.86313174];
+const LAGUERRE_ABSCISSA_9: [f64; 9] = [0.15232223,  0.80722002,  2.00513516,  3.78347397,  6.20495678, 9.37298525, 13.46623691, 18.83359779, 26.37407189];
+const LAGUERRE_ABSCISSA_10: [f64; 10] = [0.13779347, 0.72945455, 1.8083429, 3.4014337, 5.55249614, 8.33015275, 11.84378584, 16.27925783, 21.99658581, 29.92069701];
+const LAGUERRE_ABSCISSA_11: [f64; 11] = [ 0.12579644, 0.66541826, 1.64715055, 3.09113814, 5.0292844, 7.50988786, 10.605951, 14.43161376, 19.1788574, 25.21770934, 33.49719285];
+const LAGUERRE_ABSCISSA_12: [f64; 12] = [ 0.11572212, 0.61175748, 1.51261027, 2.83375134, 4.59922764, 6.84452545, 9.62131684, 13.00605499, 17.11685519, 22.15109038, 28.48796725, 37.09912104];
+const LAGUERRE_ABSCISSA_13: [f64; 13] = [0.10714239, 0.5661319, 1.39856434, 2.61659711, 4.23884593, 6.29225627, 8.81500194, 11.86140359, 15.51076204, 19.88463566, 25.18526386, 31.8003863, 40.72300867];
+const LAGUERRE_ABSCISSA_14: [f64; 14] = [0.09974751, 0.52685765, 1.30062912, 2.43080108, 3.93210282, 5.82553622, 8.14024014, 10.91649951, 14.21080501, 18.10489222, 22.72338163, 28.27298172, 35.14944366, 44.36608171];
+const LAGUERRE_ABSCISSA_15: [f64; 15] = [0.09330781, 0.49269174, 1.21559541, 2.26994953, 3.66762272, 5.42533663,  7.56591623, 10.12022857, 13.13028248, 16.65440771, 20.7764789, 25.62389423, 31.40751917, 38.53068331, 48.02608557];
+
+
+const LAGUERRE_WEIGHT_1: [f64; 1] = [1.0];
+const LAGUERRE_WEIGHT_2: [f64; 2] = [0.85355339, 0.14644661];
+const LAGUERRE_WEIGHT_3: [f64; 3] = [0.71109301, 0.27851773, 0.01038926];
+const LAGUERRE_WEIGHT_4: [f64; 4] = [6.03154104e-01, 3.57418692e-01, 3.88879085e-02, 5.39294706e-04];
+const LAGUERRE_WEIGHT_5: [f64; 5] = [5.21755611e-01, 3.98666811e-01, 7.59424497e-02, 3.61175868e-03, 2.33699724e-05];
+const LAGUERRE_WEIGHT_6: [f64; 6] = [4.58964674e-01, 4.17000831e-01, 1.13373382e-01, 1.03991975e-02, 2.61017203e-04, 8.98547906e-07];
+const LAGUERRE_WEIGHT_7: [f64; 7] = [4.09318952e-01, 4.21831278e-01, 1.47126349e-01, 2.06335145e-02, 1.07401014e-03, 1.58654643e-05, 3.17031548e-08];
+const LAGUERRE_WEIGHT_8: [f64; 8] = [3.69188589e-01, 4.18786781e-01, 1.75794987e-01, 3.33434923e-02, 2.79453624e-03, 9.07650877e-05, 8.48574672e-07, 1.04800117e-09];
+const LAGUERRE_WEIGHT_9: [f64; 9] = [3.36126422e-01, 4.11213980e-01, 1.99287525e-01, 4.74605628e-02, 5.59962661e-03, 3.05249767e-04, 6.59212303e-06, 4.11076933e-08, 3.29087403e-11];
+const LAGUERRE_WEIGHT_10: [f64; 10] = [3.08441116e-01, 4.01119929e-01, 2.18068288e-01, 6.20874561e-02, 9.50151698e-03, 7.53008389e-04, 2.82592335e-05, 4.24931398e-07, 1.83956482e-09, 9.91182722e-13];
+const LAGUERRE_WEIGHT_11: [f64; 11] = [2.84933213e-01, 3.89720890e-01, 2.32781832e-01, 7.65644535e-02, 1.43932828e-02, 1.51888085e-03, 8.51312244e-05, 2.29240388e-06, 2.48635370e-08, 7.71262693e-11, 2.88377587e-14];
+const LAGUERRE_WEIGHT_12: [f64; 12] = [2.64731371e-01, 3.77759276e-01, 2.44082011e-01, 9.04492222e-02, 2.01023812e-02, 2.66397354e-03, 2.03231593e-04, 8.36505586e-06, 1.66849388e-07, 1.34239103e-09, 3.06160164e-12, 8.14807747e-16];
+const LAGUERRE_WEIGHT_13: [f64; 13] = [2.47188708e-01, 3.65688823e-01, 2.52562420e-01, 1.03470758e-01, 2.64327544e-02, 4.22039604e-03, 4.11881770e-04, 2.35154740e-05, 7.31731162e-07, 1.10884163e-08, 6.77082669e-11, 1.15997996e-13, 2.24509320e-17];
+const LAGUERRE_WEIGHT_14: [f64; 14] = [2.31815577e-01, 3.53784692e-01, 2.58734610e-01, 1.15482894e-01, 3.31920922e-02, 6.19286944e-03, 7.39890378e-04, 5.49071947e-05, 2.40958576e-06, 5.80154398e-08, 6.81931469e-10, 3.22120775e-12, 4.22135244e-15, 6.05237502e-19];
+const LAGUERRE_WEIGHT_15: [f64; 15] = [2.18234886e-01, 3.42210178e-01, 2.63027578e-01, 1.26425818e-01, 4.02068649e-02, 8.56387780e-03, 1.21243615e-03, 1.11674392e-04, 6.45992676e-06, 2.22631691e-07, 4.22743038e-09, 3.92189727e-11, 1.45651526e-13, 1.48302705e-16, 1.60059491e-20];