I think the project should be called "NIntegrate".

BTW, that is not a serious suggestion; it is just that Wolfram Language (aka Matheamtica) has both `Integrate` and `NIntegrate` for symbolic and numeric integration, respectively.

Gaussian quadrature points are typically solved numerically. There's a good chance these ultimately came from a table.

Additionally, compile time floating-point evaluation is limited. When I looked around recently, I didn't see a rust equivalent of gcem; any kind of transcendental function evaluation (which finding Gaussian quadrature points absolutely would require) would not allow compile-time evaluation.

this was mainly to use an Iteration free method in this paper: https://www.cfm.brown.edu/faculty/gk/APMA2560/Handouts/GL_qu...

this method is much faster and simpler.

Probably but that would slow down compilation a lot.

It seems like it is lacking the functionality R's integrate has for handling infinite boundaries, but I suppose you could implement that yourself on the outside.

For what it's worth,

    use integrate::adaptive_quadrature::simpson::adaptive_simpson_method;
    use statrs::distribution::{Continuous, Normal};

    fn dnorm(x: f64) -> f64 {
        Normal::new(0.0, 1.0).unwrap().pdf(x)
    }
    
    fn main() {
        let result = adaptive_simpson_method(dnorm, -100.0, 100.0, 1e-2, 1e-8);
        println!("Result: {:?}", result);
    }

It does seem to be a usability hazard that the function being integrated is defined as a fn, rather than a Fn, as you can't pass closures that capture variables, requiring the weird dnorm definition

You will be completely blown away, then, from what Wolfram Language (aka Mathematica) can do. (When it comes to numerical integration.)

https://reference.wolfram.com/language/tutorial/NIntegrateOv...

for ]-inf, inf[ integrals, you can use Gauss Hermite method, just keep in mind to multiply your function with exp(x^2).

    use integrate::{
        gauss_quadrature::hermite::gauss_hermite_rule,
    };
    use statrs::distribution::{Continuous, Normal};

    fn dnorm(x: f64) -> f64 {
        Normal::new(0.0, 1.0).unwrap().pdf(x)* x.powi(2).exp()
    }

    fn main() {
        let n: usize = 170;
        let result = gauss_hermite_rule(dnorm, n);
        println!("Result: {:?}", result);
    }

How many evaluations of the underlying function does it make? (Hoping someone will fire up their R interpreter and find out.)

Or, probably, dnorm is a probability distribution which includes a likeliness function, and a cumulative likeliness function, etc. I bet it doesn't work on arbitrary functions.

Thanks! That’s awesome to hear—I’d love to see how your Raku numerical integration project turns out!

You can email me if you want to, I'll be happy to help.