One way to think about is that these techniques work by integrating exactly the polynomial that interpolates the function where you're sampling it, so you need to resolve the features of the function to get good accuracy.

The Nyquist sampling theorem is of course proved by considering a Fourier transform, which is given by an integral, so the relation to integration in general should not be surprising.

Well in both cases it comes down to aliasing I think (high frequency wave presents as low frequency wave with too big step size)

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.

herbie recommends `fma(h, (i + 0.5), a)`, but also doesn't report any accuracy problems with the original either

In my experience Rust is very good about using simd for loading and not great at using it automatically for math. This is from some experimentation at work and checking disassembly so ymmv

There are common library extensions for that.

It's a minimal implementation (small, lightweight) of well-known and established algorithms (fast).

Those are adjectives not comparatives.

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.

Mayve because they aren't guaranteed to be actual functions (in the mathematical sense) and could return random values