[Discussion] Adding an `init` function that could modify static variables directly

[canndrew]

If we want this feature badly enough then a better way to do it might be to allow main to take an extra argument of type StaticMut<T> where T is some type that can be const-initialized and StaticMut<T> a ZST that can be dereffed to a &'static mut T with the address known at compile time. You'd also have a Static<T> type which is the non-mut version of StaticMut<T> and a way to convert from one to the other once you've finished initializing the data.

Here's what I mean:

struct MyStaticData {
    foo: u32,
}

const fn the_actual_data_that_gets_put_in_the_binary() -> MyStaticData {
    MyStaticData {
        foo: 0,
    }
}

fn intialize_static_data_at_runtime(
    my_static_data: StaticMut<MyStaticData>,
) -> Static<MyStaticData> {
    my_static_data.foo = 123;
    my_static_data.i_dont_need_this_to_be_mut_anymore()
}

fn uses_static_data_without_really_taking_an_argument_at_runtime(
    my_static_data: Static<MyStaticData>,
) {
    println!("foo is {}", my_static_data.foo);
}

#[static_init(my_static_data = the_actual_data_that_gets_put_in_the_binary())]
fn main(my_static_data: StaticMut<MyStaticData>) {
    let my_static_data = initialize_static_data_at_runtime(my_static_data);
    uses_static_data_without_really_taking_an_argument_at_runtime(my_static_data);
}

This is zero-cost, safe, non-racey, and doesn't introduce any life-before-main. I don't think it's actually worth adding to the language though.

[Nemo157]

Seems plausible to build as a library around a static UnsafeCell<MaybeUninit<T>> and ZST tokens referencing it.

[trentj]

@NoamB

Isn't the crate graph non-cyclic? This provides a well defined init order, provided initializers are per-crate.

An acyclic graph doesn't have to have just one "proper" topological order. If Z depends on both X and Y that doesn't impose an ordering between X and Y.

[bjorn3]

@NoamB

Isn't the crate graph non-cyclic? This provides a well defined init order, provided initializers are per-crate.

Using #[no_mangle] and extern {} it is possible to call a function defined by one of the crates that depends on you. The standard library handles a lot of cycles this way.

[NoamB]

We could sort the rest alphabetically, or with any other consistent ordering. I don't really like it, but it should work.

Regarding extern and #[no_mangle] - This would mean that any crate that exports a symbol this way can't transitively rely on initializers. Maybe this can be enforced?

[farnz]

The

@NoamB

Regarding extern and #[no_mangle] - This would mean that any crate that exports a symbol this way can't transitively rely on initializers. Maybe this can be enforced?

The problem case is crates importing symbols this way, not exporting. If I import a symbol by using extern and #[no_mangle]

@Neutron3529

I strongly doubt about that.

do_maths(params, &data_table); should be slower since it passes one extra parameter(&data_table). If multiple parameters are passed, I can't believe that is no cost.

Just for fun, I changed my benchmark to show a case where parameter passing is faster than a global:

lib.rs:

use std::ops::Range;

pub const SIZE: usize = 16 * 1024;

pub static mut DATA_TABLE: [usize; SIZE] = [0; SIZE];

/// # Safety
///  This initializes DATA_TABLE
pub unsafe fn init_data_table() {
    for (index, elem) in unsafe { DATA_TABLE.iter_mut().enumerate() } {
        *elem = index * index;
    }
}

fn sum_static() -> usize {
    unsafe { DATA_TABLE }.iter().sum()
}

pub fn use_data_static(range: Range<usize>) -> usize {
    range.fold(0, |acc, elem| acc + elem * sum_static())
}

fn sum_param(data: &[usize; SIZE]) -> usize {
    data.iter().sum()
}

pub fn use_data_param(data: &[usize; SIZE], range: Range<usize>) -> usize {
    range.fold(0, |acc, elem| acc + elem * sum_param(data))
}

pub fn use_data_mixed(range: Range<usize>) -> usize {
    let data = unsafe { &DATA_TABLE };
    range.fold(0, |acc, elem| acc + elem * sum_param(data))
}

benchmark:

use benchmark_static::*;
use criterion::{black_box, criterion_group, criterion_main, Criterion};

pub fn bench_global(c: &mut Criterion) {
    unsafe {
        init_data_table();
    }

    c.bench_function("global", |b| b.iter(|| use_data_static(black_box(0..100))));
}

pub fn bench_param(c: &mut Criterion) {
    unsafe {
        init_data_table();
    }

    let data = unsafe { &DATA_TABLE };

    c.bench_function("param", |b| {
        b.iter(|| use_data_param(data, black_box(0..100)))
    });
}

pub fn bench_mixed(c: &mut Criterion) {
    unsafe {
        init_data_table();
    }

    c.bench_function("mixed", |b| b.iter(|| use_data_mixed(black_box(0..100))));
}

criterion_group!(global, bench_global);
criterion_group!(param, bench_param);
criterion_group!(mixed, bench_mixed);

criterion_main!(global, param, mixed);

And results:

     Running benches/compare.rs (target/release/deps/compare-2f28991b28549b29)
Gnuplot not found, using plotters backend
global                  time:   [401.26 µs 403.95 µs 408.56 µs]
Found 1 outliers among 100 measurements (1.00%)
  1 (1.00%) high severe

param                   time:   [147.75 µs 149.04 µs 150.49 µs]

mixed                   time:   [169.84 µs 170.34 µs 170.95 µs]
Found 5 outliers among 100 measurements (5.00%)
  1 (1.00%) high mild
  4 (4.00%) high severe

The new benchmark mixed gives you a strong hint about why passing a parameter is faster here: under Rust's aliasing rules, for as long as a borrow of DATA_TABLE exists, Rust can assume that there's no external mutation, but in theory, an external mutator could slip in between two calls to sum_static and change DATA_TABLE, and the compiler is not allowed to assume that this doesn't happen.

In contrast, use_data_mixed guarantees to the compiler that no external mutators slip in between let data = unsafe { &DATA_TABLE }; and the end of the function (that's the promise we have to uphold to make that use of unsafe sound). It can thus optimize knowing that every call to sum_param gets the same data table.

[Neutron3529]

@farnz

Just for fun, I changed my benchmark to show a case where parameter passing is faster than a global:

I firstly tried several program to explain why global parameter is slower, but as soon as I notice this benchmark, I realize that, static is not equals to static mut.

What's more, your example behave surprisingly. This is the result when I comment the whole body of init_data_table and change static mut to static

global                  time:   [232.48 µs 232.49 µs 232.49 µs]
                        change: [-64.535% -64.506% -64.482%] (p = 0.00 < 0.05)
                        Performance has improved.
Found 10 outliers among 100 measurements (10.00%)
  1 (1.00%) low severe
  3 (3.00%) high mild
  6 (6.00%) high severe

param                   time:   [232.58 µs 232.62 µs 232.67 µs]
                        change: [+17.898% +17.988% +18.068%] (p = 0.00 < 0.05)
                        Performance has regressed.
Found 9 outliers among 100 measurements (9.00%)
  7 (7.00%) high mild
  2 (2.00%) high severe

mixed                   time:   [232.49 µs 232.50 µs 232.51 µs]
                        change: [+18.066% +18.096% +18.118%] (p = 0.00 < 0.05)
                        Performance has regressed.
Found 5 outliers among 100 measurements (5.00%)
  1 (1.00%) low severe
  2 (2.00%) low mild
  2 (2.00%) high mild

Since I comment the init function, those loops should run faster, and change static mut to static is almost a no-op. But they runs slower, which is highly unexpected.

Thus, it is not my fault that make code runs slower. Global initializers should run faster.

[farnz]

@Neutron3529

Global initializers should run faster.

Why? I cannot see how the compiler will generate codepaths involving significantly fewer instructions for a global variable as opposed to a parameter - and optimization efforts in compilers tend to ignore global variables.

@Neutron3529

What's more, your example behave surprisingly. This is the result when I comment the whole body of init_data_table and change static mut to static

When I do the same (using taskset to isolate the benchmarking to a single P core on my Intel system, and with CPU clock speed locked on the test core), I get:

     Running benches/compare.rs (target/release/deps/compare-2f28991b28549b29)
Gnuplot not found, using plotters backend
global                  time:   [293.51 µs 294.19 µs 294.73 µs]
                        change: [-5.6473% -5.2971% -4.8791%] (p = 0.00 < 0.05)
                        Performance has improved.

param                   time:   [164.15 µs 164.18 µs 164.22 µs]
                        change: [-0.2689% -0.1116% +0.0205%] (p = 0.14 > 0.05)
                        No change in performance detected.
Found 9 outliers among 100 measurements (9.00%)
  4 (4.00%) high mild
  5 (5.00%) high severe

mixed                   time:   [176.41 µs 176.82 µs 177.40 µs]
                        change: [-18.967% -18.782% -18.551%] (p = 0.00 < 0.05)
                        Performance has improved.
Found 16 outliers among 100 measurements (16.00%)
  1 (1.00%) low severe
  6 (6.00%) high mild
  9 (9.00%) high severe

Because the runtime is so low (around 10,000 CPU clock cycles per iteration for the slowest of the 3), it's quite hard to isolate all the things that can affect the final result. But even so, passing a parameter around remains the fastest, in large part because the static cannot be optimized as well as a parameter passed by reference, because the aliasing (and hence external mutation) rules are different for statics and for shared references.

[NoamB]

@farnz

The problem case is crates importing symbols this way, not exporting. If I import a symbol by using extern and #[no_mangle]

If the exported symbol doesn't need any initializer to run before it is called, then it can be called at any point. This also allows its dependencies to call into it.

What limitation can be placed on the importing crate?

If we assume that we don't need to impose any limitation on the exporting crate, the following becomes possible:

If we have crates A, B, C, D where:

• A has a static initializer, and exports symbol fn needs_init(). needs_init needs the initizlizer to run for it to be sound to call (Allowed by our initial assumption).
• B links to A::needs_init and calls it from fn calls_dependent()
• C depends on B, and in its initializer, calls B::calls_dependent()

If C's initializer is called before A's (for example, if A depends on C), then needs_init is called before its static initializer.

[farnz]

@NoamB

What limitation can be placed on the importing crate?

Restricting the importing crate is hard - effectively, what it's doing is going around the normal linking mechanisms to import a symbol, bypassing the normal mechanisms that are supposed to make sure that A's initializers run before C can call B::calls_dependent().

There is a certain extent to which this sort of trickery is "play silly games, win stupid prizes", but it still needs considering.

[NoamB]

I think that restricting the importing crate is intractable - that's why I suggested restricting the exporting crate. This would allow low level crates which know about each other to do the cyclic dependency trickery.

More flexibility can be gained by allowing a crate which exports a symbol to explicitly allow some of its dependencies to have initializers.

[farnz]

@NoamB

I think that restricting the importing crate is intractable - that's why I suggested restricting the exporting crate. This would allow low level crates which know about each other to do the cyclic dependency trickery.

That's the underlying challenge, though - it's the importing crates that have all the control and power, right up until they call code in the exporting crate. Nothing stops C from calling A::needs_init or B::calls_dependent before A's static initialization code runs, especially once you consider environments like WebAssembly.

This is what lazy initialization works around, at the expense of 3 instructions per global access - it ensures that when C calls A::needs_init or B::calls_dependent, A's initializer runs, because it's run lazily on demand. And if having the first call to A::needs_init or B::calls_dependent run slow is not always acceptable, A can export a eager_init function that does all the lazy init work up-front.

Part of the question, though, is how much work we actually ought to do for this case - if you can pass a global context by reference, instead of having it found lurking in the environment, you get performance benefits where it enables extra optimizations, and you get to have multiple global contexts in one program. If you can't pass a reference to the context, but can add a parameter, you get the ZST init token trick. If you can afford 3 instructions when looking up your context, you can use OnceLock or LazyLock.

And there are at least two directions that I'm aware of that people have thought seriously about which would be more general. One is just doing a better job of VRA, so that the 3 instruction penalty for OnceLock and friends goes away if it's possible for the compiler to see that it's been initialized already, just as the 3 instruction penalty for slice indexing goes away if the compiler can see that the slice is clearly large enough; the other is capabilities/contexts which would provide a way to pass the global context around without requiring it to be a parameter.

[Neutron3529]

@farnz

because the aliasing (and hence external mutation) rules are different for statics and for shared references

have you comment initialization part like that?

pub static DATA_TABLE: [usize; SIZE] = [0; SIZE];

/// # Safety
///  This initializes DATA_TABLE
pub unsafe fn init_data_table() {
//    for (index, elem) in unsafe { DATA_TABLE.iter_mut().enumerate() } {
//        *elem = index * index;
//    }
}

I tried comment init_data_table but keep the static mut unchanged, then I could got a huge speedup which lowering time to <100 us, which shouldn't be regard as "No change in performance detected"

@NoamB

If C's initializer is called before A's (for example, if A depends on C), then needs_init is called before its static initializer.

@NoamB

I think that restricting the importing crate is intractable

@farnz

Nothing stops C from calling A::needs_init or B::calls_dependent before A's static initialization code runs, especially once you consider environments like WebAssembly.

restricting the importing crate is tractable since rust does not allow cycle dependencies.

In your case, B dependent on A, and C dependent on B, thus A's initializer runs first, B the second, and C the last.

---

What's more, an intresting thing is found:

@farnz

fn sum_static() -> usize {
    unsafe { DATA_TABLE }.iter().sum()
}

@farnz

let data = unsafe { &DATA_TABLE };

You have actually found what affect the calculating speed, but you omit that and simply said it is static variables that runs slower.

The fact is behind the reference.

unsafe {DATA_TABLE} // Copy
unsafe {&DATA_TABLE} // Reference only

you create copies, thus sum_static is slower.

The change could be simple:

fn sum_static() -> usize {
    unsafe { &DATA_TABLE }.iter().sum()
}

using such method we could calculate sum_static and sum_param at the same speed.

[farnz]

@Neutron3529

@farnz

because the aliasing (and hence external mutation) rules are different for statics and for shared references

have you comment initialization part like that?

pub static DATA_TABLE: [usize; SIZE] = [0; SIZE];

/// # Safety
///  This initializes DATA_TABLE
pub unsafe fn init_data_table() {
//    for (index, elem) in unsafe { DATA_TABLE.iter_mut().enumerate() } {
//        *elem = index * index;
//    }
}

I tried comment init_data_table but keep the static mut unchanged, then I could got a huge speedup which lowering time to <100 us, which shouldn't be regard as "No change in performance detected"

I didn't just comment out the body of init_data_table - I removed the entire function. I suspect the speed-up you're seeing is from failing to lock the CPU clock for the core you're testing on, and not a consequence of code changes.

@Neutron3529

@NoamB

If C's initializer is called before A's (for example, if A depends on C), then needs_init is called before its static initializer.

@NoamB

I think that restricting the importing crate is intractable

@farnz

Nothing stops C from calling A::needs_init or B::calls_dependent before A's static initialization code runs, especially once you consider environments like WebAssembly.

restricting the importing crate is tractable since rust does not allow cycle dependencies.

In your case, B dependent on A, and C dependent on B, thus A's initializer runs first, B the second, and C the last.

This is for the case where C has been built with disabled running of static initializers, but A and B were not built that way - and we got at the symbols in A and B via a route like extern, instead of a proper dependency build.

This is an edge case, but we need to thing about it.

@Neutron3529

The fact is behind the reference.

unsafe {DATA_TABLE} // Copy
unsafe {&DATA_TABLE} // Reference only

you create copies, thus sum_static is slower.

This is not true in the static case - see Compiler Explorer where the generated code is identical for both versions.

It is true in the static mut case, because that's how Rust handles the risk of external changes. But even with this change to the static mut case, I still see the parameter passing case being faster than the global case - it just becomes the same timings as the static version, where the global version is about 100 µs slower than parameter passing.

[Neutron3529]

@farnz

But even with this change to the static mut case, I still see the parameter passing case being faster than the global case

This because something is out of controlled.

I could easily wrote code that shows parameter passing is slower, but that is not "one method just faster than other", it is just some noise that we couldn't control.

code:

pub static mut coef0:f64=-1.3;
pub static mut coef1:f64=0.5;
pub static mut coef2:f64=0.3;
pub static mut coef3:f64=0.2;
pub static mut coef4:f64=0.1;
pub static mut coef5:f64=0.05;
pub static mut coef6:f64=0.03;
pub static mut coef7:f64=0.02;
pub static mut coef8:f64=0.01;
pub static mut coef9:f64=0.001;


#[inline(never)] // it is faster enforce not inline, don't know why.
fn fn_static(x:f64) -> f64 {
    unsafe{
        (((((((((coef9*x+coef8)*x+coef7)*x+coef6)*x+coef5)*x+coef4)*x+coef3)*x+coef2)*x+coef1)*x+coef0)
    }
}
#[inline(never)]
pub fn use_fn_static(mut x:f64, rep:usize) -> f64 {
    for _ in 0..rep {
        x=fn_static(x).exp();
    }
    x
}
// #[inline(never)]
fn fn_param(x:f64,c0:f64,c1:f64,c2:f64,c3:f64,c4:f64,c5:f64,c6:f64,c7:f64,c8:f64,c9:f64) -> f64 {
    (((((((((c9*x+c8)*x+c7)*x+c6)*x+c5)*x+c4)*x+c3)*x+c2)*x+c1)*x+c0)
}

// #[inline(never)]
pub fn use_fn_param(mut x:f64,rep:usize,c0:f64,c1:f64,c2:f64,c3:f64,c4:f64,c5:f64,c6:f64,c7:f64,c8:f64,c9:f64) -> f64 {
    for _ in 0..rep {
        x=fn_param(x,c0,c1,c2,c3,c4,c5,c6,c7,c8,c9).exp();
    }
    x
}

bench:

use benchmark_static::*;
use criterion::{black_box, criterion_group, criterion_main, Criterion};


pub fn bench_global(c: &mut Criterion) {
    c.bench_function("global", |b| b.iter(|| use_fn_static(1.0,black_box(100))));
}
pub fn bench_param(c: &mut Criterion) {
    let [c0,c1,c2,c3,c4,c5,c6,c7,c8,c9]=unsafe{[coef0,coef1,coef2,coef3,coef4,coef5,coef6,coef7,coef8,coef9]};
    c.bench_function("param", |b| b.iter(|| use_fn_param(1.0,black_box(100),c0,c1,c2,c3,c4,c5,c6,c7,c8,c9)));
}

criterion_group!(global, bench_global);
criterion_group!(param, bench_param);

criterion_main!(global, param);

They should run at the same speed, but if we put #[inline(never)] correctly (to be precise, on fn_static), use_fn_static would got ~1% performance gain.

global                  time:   [4.5881 µs 4.5887 µs 4.5895 µs]
Found 3 outliers among 100 measurements (3.00%)
  1 (1.00%) high mild
  2 (2.00%) high severe

param                   time:   [4.6124 µs 4.6126 µs 4.6128 µs]
Found 4 outliers among 100 measurements (4.00%)
  1 (1.00%) low mild
  1 (1.00%) high mild
  2 (2.00%) high severe

Actually, benchmark could say nothing.

[farnz]

This is getting away from the original point - my point is that given two variants, where one gets the global context as a parameter, and the other gets the global context via a static, the version with the static will never be faster than the version with the parameter. I've shown via benchmarks that this is true - I've also shown a benchmark where the version with the parameter is faster, because (looking at the generated code on my machine) Rust has been able to inline in both cases, and then determine that sum_param is a constant value after inlining, whereas sum_static cannot be because of the need to reload &DATA_TABLE each time it does the addition.

I've repeated your test on my machine, where the clock speed is locked down, and the test has the same core each time (and that core's sibling is offlined), and I see them both taking almost exactly the same amount of time - with param usually being slightly faster.

Can you show a case where the two versions differ only in that one gets passed a reference to the global, and other finds it via a static, and passing via a static is consistently much faster?

[Neutron3529]

@farnz

Rust has been able to inline in both cases

Will, this is not true, since we could have REALLY COMPLEX function that could not be inlined.

For simple example, I have a program that computing bivariate normal probabilities, which is firstly wrote in fortran (thus source code is available in orig.rs ), It uses massive of consts which need no initialize, but fortunately, we could modify those constant for benchmark reason.

Full program is provided as follows:

src/lib.rs

pub mod orig;
pub mod r#static;
pub mod param;

src/orig.rs

/// Rust version mvbvn
/// A function for computing bivariate normal probabilities
pub fn mvbvn(mut li: f64, mut lj: f64, mut ui: f64, mut uj: f64, mut corr: f64) -> f64 {
    if li >= ui || lj >= uj {
        0.0
    } else {
        if li.abs() > ui.abs() {
            (li, ui, corr) = (-ui, -li, -corr)
        }
        if lj.abs() > uj.abs() {
            (lj, uj, corr) = (-uj, -lj, -corr)
        }
        //         if li==f64::NEG_INFINITY {MVPHI(uj)-MVPHI(lj)}
        mvbvu(li, lj, corr)
            - if uj != f64::INFINITY {
                mvbvu(li, uj, corr)
            } else {
                0.0
            }
            - if ui != f64::INFINITY {
                mvbvu(ui, lj, corr)
            } else {
                0.0
            }
            + if ui + uj != f64::INFINITY {
                mvbvu(ui, uj, corr)
            } else {
                0.0
            }
    }
}

// // benchmark only
// fn garbage(li:f64,lj:f64,ui:f64)->f64{
//     let mut ans=0.;
//     for _ in 0..1000 {
//         ans-=li;
//         ans=ans.exp()+lj;
//         ans=ans.ln()*ui;
//     }
//     ans
// }

// #[inline(always)]
// fn mvbvu(sh:f64,sk:f64,r:f64)->f64{
//     unsafe {MVBVU(&sh,&sk,&r)}
// }

///       DOUBLE PRECISION FUNCTION MVBVU( SH, SK, R )
/// *
/// *     A function for computing bivariate normal probabilities;
/// *       developed using
/// *         Drezner, Z. and Wesolowsky, G. O. (1989),
/// *         On the Computation of the Bivariate Normal Integral,
/// *         J. Stat. Comput. Simul.. 35 pp. 101-107.
/// *       with extensive modications for double precisions by
/// *         Alan Genz and Yihong Ge
/// *         Department of Mathematics
/// *         Washington State University
/// *         Pullman, WA 99164-3113
/// *         Email : alangenz@wsu.edu
/// *
/// * BVN - calculate the probability that X is larger than SH and Y is
/// *       larger than SK.
/// *
/// * Parameters
/// *
/// *   SH  REAL, integration limit
/// *   SK  REAL, integration limit
/// *   R   REAL, correlation coefficient
/// *   LG  INTEGER, number of Gauss Rule Points and Weights
/// *
fn mvbvu(sh: f64, sk: f64, r: f64) -> f64 {
    //      DOUBLE PRECISION BVN, SH, SK, R, ZERO, TWOPI
    //     let mut bvn=0f64;
    //      INTEGER I, LG, NG
    let lg;
    let ng;
    //      PARAMETER ( ZERO = 0, TWOPI = 6.283185307179586D0 )
    const TWOPI: f64 = 6.283185307179586;
    //      DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS
    //      DOUBLE PRECISION MVPHI, SN, ASR, H, K, BS, HS, HK
    //      SAVE X, W
    // *     Gauss Legendre Points and Weights, N =  6
    //       DATA ( W(I,1), X(I,1), I = 1, 3 ) /
    //      *  0.1713244923791705D+00,-0.9324695142031522D+00,
    //      *  0.3607615730481384D+00,-0.6612093864662647D+00,
    //      *  0.4679139345726904D+00,-0.2386191860831970D+00/
    // *     Gauss Legendre Points and Weights, N = 12
    //       DATA ( W(I,2), X(I,2), I = 1, 6 ) /
    //      *  0.4717533638651177D-01,-0.9815606342467191D+00,
    //      *  0.1069393259953183D+00,-0.9041172563704750D+00,
    //      *  0.1600783285433464D+00,-0.7699026741943050D+00,
    //      *  0.2031674267230659D+00,-0.5873179542866171D+00,
    //      *  0.2334925365383547D+00,-0.3678314989981802D+00,
    //      *  0.2491470458134029D+00,-0.1252334085114692D+00/
    // *     Gauss Legendre Points and Weights, N = 20
    //       DATA ( W(I,3), X(I,3), I = 1, 10 ) /
    //      *  0.1761400713915212D-01,-0.9931285991850949D+00,
    //      *  0.4060142980038694D-01,-0.9639719272779138D+00,
    //      *  0.6267204833410906D-01,-0.9122344282513259D+00,
    //      *  0.8327674157670475D-01,-0.8391169718222188D+00,
    //      *  0.1019301198172404D+00,-0.7463319064601508D+00,
    //      *  0.1181945319615184D+00,-0.6360536807265150D+00,
    //      *  0.1316886384491766D+00,-0.5108670019508271D+00,
    //      *  0.1420961093183821D+00,-0.3737060887154196D+00,
    //      *  0.1491729864726037D+00,-0.2277858511416451D+00,
    //      *  0.1527533871307259D+00,-0.7652652113349733D-01/

    const W: [[f64; 10]; 3] = [
        [
            0.1713244923791705,
            0.3607615730481384,
            0.4679139345726904,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
        ],
        [
            0.04717533638651177,
            0.1069393259953183,
            0.1600783285433464,
            0.2031674267230659,
            0.2334925365383547,
            0.2491470458134029,
            0.0,
            0.0,
            0.0,
            0.0,
        ],
        [
            0.01761400713915212,
            0.04060142980038694,
            0.06267204833410905,
            0.08327674157670475,
            0.1019301198172404,
            0.1181945319615184,
            0.1316886384491766,
            0.1420961093183821,
            0.1491729864726037,
            0.1527533871307259,
        ],
    ]; // python script (ignore the leading //): [[float(i[0]) for i in i]+[0.0]*(10-len(i)) for i in [[i.strip().replace('D','e').split(',') for i in i.split('*')[1:] for i in i.strip().split('\n')] for i in "\n".join([a for a in a.split('\n') if 'D+' in a]).split('/')[:3]]]
    const X: [[f64; 10]; 3] = [
        [
            -0.9324695142031522,
            -0.6612093864662647,
            -0.238619186083197,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
        ],
        [
            -0.9815606342467191,
            -0.904117256370475,
            -0.769902674194305,
            -0.5873179542866171,
            -0.3678314989981802,
            -0.1252334085114692,
            0.0,
            0.0,
            0.0,
            0.0,
        ],
        [
            -0.9931285991850949,
            -0.9639719272779138,
            -0.912234428251326,
            -0.8391169718222188,
            -0.7463319064601508,
            -0.636053680726515,
            -0.5108670019508271,
            -0.3737060887154196,
            -0.2277858511416451,
            -0.07652652113349732,
        ],
    ]; // python script (ignore the leading //): [[float(i[1]) for i in i]+[0.0]*(10-len(i)) for i in [[i.strip().replace('D','e').split(',') for i in i.split('*')[1:] for i in i.strip().split('\n')] for i in "\n".join([a for a in a.split('\n') if 'D+' in a]).split('/')[:3]]]
    //       IF ( ABS(R) .LT. 0.3 ) THEN
    //          NG = 1
    //          LG = 3
    //       ELSE IF ( ABS(R) .LT. 0.75 ) THEN
    //          NG = 2
    //          LG = 6
    //       ELSE
    //          NG = 3
    //          LG = 10
    //       ENDIF
    if r.abs() < 0.3 {
        ng = 0;
        lg = 3;
    } else if r.abs() < 0.75 {
        ng = 1;
        lg = 6;
    } else {
        ng = 2;
        lg = 10;
    }
    //       H = SH
    //       K = SK
    let h = sh;
    let mut k = sk;
    //       HK = H*K
    let mut hk = h * k;
    //       BVN = 0
    //       IF ( ABS(R) .LT. 0.925 ) THEN
    if r.abs() < 0.925 {
        //          HS = ( H*H + K*K )/2
        let hs = (h * h + k * k) / 2.0;
        //          ASR = ASIN(R)
        let asr = r.asin();
        //          DO I = 1, LG
        //         for i in 0..lg{
        // //             SN = SIN(ASR*( X(I,NG)+1 )/2)
        //             sn=(asr*(( X[ng][i]+1f64 )/2f64)).sin();
        // //             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
        //             bvn+=W[ng][i]*((sn*hk-hs)/(1-sn*sn)).exp()
        // //             SN = SIN(ASR*(-X(I,NG)+1 )/2)
        // //             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
        // //          END DO
        //         }
        //          BVN = BVN*ASR/(2*TWOPI) + MVPHI(-H)*MVPHI(-K)
        (0..lg)
            .map(|i| {
                let sn1 = (asr * ((X[ng][i] + 1.0) / 2.0)).sin();
                let sn2 = (asr * ((-X[ng][i] + 1.0) / 2.0)).sin();
                W[ng][i]
                    * (((sn1 * hk - hs) / (1.0 - sn1 * sn1)).exp()
                        + ((sn2 * hk - hs) / (1.0 - sn2 * sn2)).exp())
            })
            .sum::<f64>()
            * asr
            / (TWOPI * 2.0)
            + mvphi(-h) * mvphi(-k)
    } else {
        //       ELSE
        //          IF ( R .LT. 0 ) THEN
        //             K = -K
        //             HK = -HK
        //          ENDIF
        if r < 0.0 {
            k = -k;
            hk = -hk;
        }
        //          IF ( ABS(R) .LT. 1 ) THEN
        let bvn = if r.abs() < 1.0 {
            //             AS = ( 1 - R )*( 1 + R )
            let r#as = (1.0 - r) * (1.0 + r);
            //             A = SQRT(AS)
            let a = r#as.sqrt();
            //             BS = ( H - K )**2
            let bs = (h - k) * (h - k);
            //             C = ( 4 - HK )/8
            //             D = ( 12 - HK )/16
            let c = (4.0 - hk) / 8.0;
            let d = (12.0 - hk) / 16.0;
            //             BVN = A*EXP( -(BS/AS + HK)/2 )
            //      +             *( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 )
            let mut bvn = a
                * (-(bs / r#as + hk) / 2.0).exp()
                * (1.0 - c * (bs - r#as) * (1.0 - d * bs / 5.0) / 3.0 + c * d * r#as * r#as / 5.0);
            //             IF ( HK .GT. -160 ) THEN
            if hk > -160.0 {
                //                B = SQRT(BS)
                let b = bs.sqrt();
                //                BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*MVPHI(-B/A)*B
                //      +                    *( 1 - C*BS*( 1 - D*BS/5 )/3 )
                bvn -= (-hk / 2.0).exp()
                    * TWOPI.sqrt()
                    * mvphi(-b / a)
                    * b
                    * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0);
                //             ENDIF
            }
            //             A = A/2
            let a = a / 2.0;
            //             DO I = 1, LG
            //                XS = ( A*(X(I,NG)+1) )**2
            //                RS = SQRT( 1 - XS )
            //                BVN = BVN + A*W(I,NG)*
            //      +              ( EXP( -BS/(2*XS) - HK/(1+RS) )/RS
            //      +              - EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) )
            //                XS = AS*(-X(I,NG)+1)**2/4
            //                RS = SQRT( 1 - XS )
            //                BVN = BVN + A*W(I,NG)*EXP( -(BS/XS + HK)/2 )
            //      +                    *( EXP( -HK*(1-RS)/(2*(1+RS)) )/RS
            //      +                       - ( 1 + C*XS*( 1 + D*XS ) ) )
            //             END DO
            -(bvn
                + (0..lg)
                    .map(|i| {
                        let xs = (a * (X[ng][i as usize] + 1.0)) * (a * (X[ng][i as usize] + 1.0));
                        let rs = (1.0 - xs).sqrt();
                        let bvn = a
                            * W[ng][i as usize]
                            * ((-bs / (2.0 * xs) - hk / (1.0 + rs)).exp() / rs
                                - (-(bs / xs + hk) / 2.0).exp() * (1.0 + c * xs * (1.0 + d * xs)));
                        let xs =
                            r#as * (-X[ng][i as usize] + 1.0) * (-X[ng][i as usize] + 1.0) / 4.0;
                        let rs = (1.0 - xs).sqrt();
                        bvn + a
                            * W[ng][i as usize]
                            * (-(bs / xs + hk) / 2.0).exp()
                            * ((-hk * (1.0 - rs) / (2.0 * (1.0 + rs))).exp() / rs
                                - (1.0 + c * xs * (1.0 + d * xs)))
                    })
                    .sum::<f64>())
                / TWOPI
        //             BVN = -BVN/TWOPI

        //          ENDIF
        } else {
            0.0
        };
        //          IF ( R .GT. 0 ) BVN =  BVN + MVPHI( -MAX( H, K ) )
        //          IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, MVPHI(-H) - MVPHI(-K) )
        if r > 0.0 {
            bvn + mvphi(-(h.max(k)))
        } else {
            -bvn + 0f64.max(mvphi(-h) - mvphi(-k))
        }
        //       ENDIF
    }
    //       MVBVU = BVN
    //       END
}

///       DOUBLE PRECISION FUNCTION MVPHI(Z)
/// *
/// *     Normal distribution probabilities accurate to 1d-15.
/// *     Reference: J.L. Schonfelder, Math Comp 32(1978), pp 1232-1240.
/// *
fn mvphi(z: f64) -> f64 {
    //       INTEGER I, IM
    //       DOUBLE PRECISION A(0:43), BM, B, BP, P, RTWO, T, XA, Z
    //       PARAMETER( RTWO = 1.414213562373095048801688724209D0, IM = 24 )
    const IM: i32 = 24;
    const RTWO: f64 = 1.414213562373095048801688724209;
    //       SAVE A
    //       DATA ( A(I), I = 0, 43 )/
    //      &    6.10143081923200417926465815756D-1,
    //      &   -4.34841272712577471828182820888D-1,
    //      &    1.76351193643605501125840298123D-1,
    //      &   -6.0710795609249414860051215825D-2,
    //      &    1.7712068995694114486147141191D-2,
    //      &   -4.321119385567293818599864968D-3,
    //      &    8.54216676887098678819832055D-4,
    //      &   -1.27155090609162742628893940D-4,
    //      &    1.1248167243671189468847072D-5, 3.13063885421820972630152D-7,
    //      &   -2.70988068537762022009086D-7, 3.0737622701407688440959D-8,
    //      &    2.515620384817622937314D-9, -1.028929921320319127590D-9,
    //      &    2.9944052119949939363D-11, 2.6051789687266936290D-11,
    //      &   -2.634839924171969386D-12, -6.43404509890636443D-13,
    //      &    1.12457401801663447D-13, 1.7281533389986098D-14,
    //      &   -4.264101694942375D-15, -5.45371977880191D-16,
    //      &    1.58697607761671D-16, 2.0899837844334D-17,
    //      &   -5.900526869409D-18, -9.41893387554D-19, 2.14977356470D-19,
    //      &    4.6660985008D-20, -7.243011862D-21, -2.387966824D-21,
    //      &    1.91177535D-22, 1.20482568D-22, -6.72377D-25, -5.747997D-24,
    //      &   -4.28493D-25, 2.44856D-25, 4.3793D-26, -8.151D-27, -3.089D-27,
    //      &    9.3D-29, 1.74D-28, 1.6D-29, -8.0D-30, -2.0D-30 /
    // *
    const A: [f64; 44] = [
        6.10143081923200417926465815756e-1,
        -4.34841272712577471828182820888e-1,
        1.76351193643605501125840298123e-1,
        -6.0710795609249414860051215825e-2,
        1.7712068995694114486147141191e-2,
        -4.321119385567293818599864968e-3,
        8.54216676887098678819832055e-4,
        -1.27155090609162742628893940e-4,
        1.1248167243671189468847072e-5,
        3.13063885421820972630152e-7,
        -2.70988068537762022009086e-7,
        3.0737622701407688440959e-8,
        2.515620384817622937314e-9,
        -1.028929921320319127590e-9,
        2.9944052119949939363e-11,
        2.6051789687266936290e-11,
        -2.634839924171969386e-12,
        -6.43404509890636443e-13,
        1.12457401801663447e-13,
        1.7281533389986098e-14,
        -4.264101694942375e-15,
        -5.45371977880191e-16,
        1.58697607761671e-16,
        2.0899837844334e-17,
        -5.900526869409e-18,
        -9.41893387554e-19,
        2.14977356470e-19,
        4.6660985008e-20,
        -7.243011862e-21,
        -2.387966824e-21,
        1.91177535e-22,
        1.20482568e-22,
        -6.72377e-25,
        -5.747997e-24,
        -4.28493e-25,
        2.44856e-25,
        4.3793e-26,
        -8.151e-27,
        -3.089e-27,
        9.3e-29,
        1.74e-28,
        1.6e-29,
        -8.0e-30,
        -2.0e-30,
    ];
    //       XA = ABS(Z)/RTWO
    let xa = z.abs() / RTWO;
    //       IF ( XA .GT. 100 ) THEN
    let p = if xa >= 100.0 {
        0.0
    //          P = 0
    //       ELSE
    } else {
        //          T = ( 8*XA - 30 ) / ( 4*XA + 15 )
        let t = (8.0 * xa - 30.0) / (4.0 * xa + 15.0);
        let mut bm = 0.0;
        let mut b = 0.0;
        let mut bp = 0.0;
        //          BM = 0
        //          B  = 0
        //          DO I = IM, 0, -1
        //             BP = B
        //             B  = BM
        //             BM = T*B - BP  + A(I)
        //          END DO
        for i in (0..=IM).rev() {
            bp = b;
            b = bm;
            bm = t * b - bp + A[i as usize]
        }
        //          P = EXP( -XA*XA )*( BM - BP )/4
        (-xa * xa).exp() * (bm - bp) / 4.0
        //       END IF
    };
    //       IF ( Z .GT. 0 ) P = 1 - P
    //       MVPHI = P
    if z > 0.0 {
        1.0 - p
    } else {
        p
    }
    //       END
}

src/static.rs

static W0: [f64; 3] =
        [
            0.1713244923791705,
            0.3607615730481384,
            0.4679139345726904,
        ];
static W1:[f64;6]=
        [
            0.04717533638651177,
            0.1069393259953183,
            0.1600783285433464,
            0.2031674267230659,
            0.2334925365383547,
            0.2491470458134029,
        ];
static W2:[f64;10]=
        [
            0.01761400713915212,
            0.04060142980038694,
            0.06267204833410905,
            0.08327674157670475,
            0.1019301198172404,
            0.1181945319615184,
            0.1316886384491766,
            0.1420961093183821,
            0.1491729864726037,
            0.1527533871307259,
        ];
static X0: [f64; 3] =
        [
            -0.9324695142031522,
            -0.6612093864662647,
            -0.238619186083197,
        ];
static X1:[f64; 6] = [
            -0.9815606342467191,
            -0.904117256370475,
            -0.769902674194305,
            -0.5873179542866171,
            -0.3678314989981802,
            -0.1252334085114692,
        ];
static X2:[f64;10]=[
            -0.9931285991850949,
            -0.9639719272779138,
            -0.912234428251326,
            -0.8391169718222188,
            -0.7463319064601508,
            -0.636053680726515,
            -0.5108670019508271,
            -0.3737060887154196,
            -0.2277858511416451,
            -0.07652652113349732,
        ];

static A: [f64; 44] = [
    6.10143081923200417926465815756e-1,
    -4.34841272712577471828182820888e-1,
    1.76351193643605501125840298123e-1,
    -6.0710795609249414860051215825e-2,
    1.7712068995694114486147141191e-2,
    -4.321119385567293818599864968e-3,
    8.54216676887098678819832055e-4,
    -1.27155090609162742628893940e-4,
    1.1248167243671189468847072e-5,
    3.13063885421820972630152e-7,
    -2.70988068537762022009086e-7,
    3.0737622701407688440959e-8,
    2.515620384817622937314e-9,
    -1.028929921320319127590e-9,
    2.9944052119949939363e-11,
    2.6051789687266936290e-11,
    -2.634839924171969386e-12,
    -6.43404509890636443e-13,
    1.12457401801663447e-13,
    1.7281533389986098e-14,
    -4.264101694942375e-15,
    -5.45371977880191e-16,
    1.58697607761671e-16,
    2.0899837844334e-17,
    -5.900526869409e-18,
    -9.41893387554e-19,
    2.14977356470e-19,
    4.6660985008e-20,
    -7.243011862e-21,
    -2.387966824e-21,
    1.91177535e-22,
    1.20482568e-22,
    -6.72377e-25,
    -5.747997e-24,
    -4.28493e-25,
    2.44856e-25,
    4.3793e-26,
    -8.151e-27,
    -3.089e-27,
    9.3e-29,
    1.74e-28,
    1.6e-29,
    -8.0e-30,
    -2.0e-30,
];

pub fn mvbvn(mut li: f64, mut lj: f64, mut ui: f64, mut uj: f64, mut corr: f64) -> f64 {
    if li >= ui || lj >= uj {
        0.0
    } else {
        if li.abs() > ui.abs() {
            (li, ui, corr) = (-ui, -li, -corr)
        }
        if lj.abs() > uj.abs() {
            (lj, uj, corr) = (-uj, -lj, -corr)
        }
        //         if li==f64::NEG_INFINITY {MVPHI(uj)-MVPHI(lj)}
        mvbvu(li, lj, corr)
            - if uj != f64::INFINITY {
                mvbvu(li, uj, corr)
            } else {
                0.0
            }
            - if ui != f64::INFINITY {
                mvbvu(ui, lj, corr)
            } else {
                0.0
            }
            + if ui + uj != f64::INFINITY {
                mvbvu(ui, uj, corr)
            } else {
                0.0
            }
    }
}

fn mvbvu(sh: f64, sk: f64, r: f64) -> f64 {
    //      DOUBLE PRECISION BVN, SH, SK, R, ZERO, TWOPI
    //     let mut bvn=0f64;
    //      INTEGER I, LG, NG
    //      PARAMETER ( ZERO = 0, TWOPI = 6.283185307179586D0 )
    const TWOPI: f64 = 6.283185307179586;
    //      DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS
    //      DOUBLE PRECISION MVPHI, SN, ASR, H, K, BS, HS, HK
    //      SAVE X, W
    // *     Gauss Legendre Points and Weights, N =  6
    //       DATA ( W(I,1), X(I,1), I = 1, 3 ) /
    //      *  0.1713244923791705D+00,-0.9324695142031522D+00,
    //      *  0.3607615730481384D+00,-0.6612093864662647D+00,
    //      *  0.4679139345726904D+00,-0.2386191860831970D+00/
    // *     Gauss Legendre Points and Weights, N = 12
    //       DATA ( W(I,2), X(I,2), I = 1, 6 ) /
    //      *  0.4717533638651177D-01,-0.9815606342467191D+00,
    //      *  0.1069393259953183D+00,-0.9041172563704750D+00,
    //      *  0.1600783285433464D+00,-0.7699026741943050D+00,
    //      *  0.2031674267230659D+00,-0.5873179542866171D+00,
    //      *  0.2334925365383547D+00,-0.3678314989981802D+00,
    //      *  0.2491470458134029D+00,-0.1252334085114692D+00/
    // *     Gauss Legendre Points and Weights, N = 20
    //       DATA ( W(I,3), X(I,3), I = 1, 10 ) /
    //      *  0.1761400713915212D-01,-0.9931285991850949D+00,
    //      *  0.4060142980038694D-01,-0.9639719272779138D+00,
    //      *  0.6267204833410906D-01,-0.9122344282513259D+00,
    //      *  0.8327674157670475D-01,-0.8391169718222188D+00,
    //      *  0.1019301198172404D+00,-0.7463319064601508D+00,
    //      *  0.1181945319615184D+00,-0.6360536807265150D+00,
    //      *  0.1316886384491766D+00,-0.5108670019508271D+00,
    //      *  0.1420961093183821D+00,-0.3737060887154196D+00,
    //      *  0.1491729864726037D+00,-0.2277858511416451D+00,
    //      *  0.1527533871307259D+00,-0.7652652113349733D-01/

    //       IF ( ABS(R) .LT. 0.3 ) THEN
    //          NG = 1
    //          LG = 3
    //       ELSE IF ( ABS(R) .LT. 0.75 ) THEN
    //          NG = 2
    //          LG = 6
    //       ELSE
    //          NG = 3
    //          LG = 10
    //       ENDIF
    //       H = SH
    //       K = SK
    let h = sh;
    let mut k = sk;
    //       HK = H*K
    let mut hk = h * k;
    //       BVN = 0
    //       IF ( ABS(R) .LT. 0.925 ) THEN
    if r.abs() < 0.925 {
        //          HS = ( H*H + K*K )/2
        let hs = (h * h + k * k) / 2.0;
        //          ASR = ASIN(R)
        let asr = r.asin();
        //          DO I = 1, LG
        //         for i in 0..lg{
        // //             SN = SIN(ASR*( X(I,NG)+1 )/2)
        //             sn=(asr*(( X[ng][i]+1f64 )/2f64)).sin();
        // //             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
        //             bvn+=W[ng][i]*((sn*hk-hs)/(1-sn*sn)).exp()
        // //             SN = SIN(ASR*(-X(I,NG)+1 )/2)
        // //             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
        // //          END DO
        //         }
        //          BVN = BVN*ASR/(2*TWOPI) + MVPHI(-H)*MVPHI(-K)
        (if r.abs() < 0.3 {
            X0.iter().zip(W0.iter())
        } else if r.abs() < 0.75 {
            X1.iter().zip(W1.iter())
        } else {
            X2.iter().zip(W2.iter())
        }).map(|(x,w)| {
                let sn1 = (asr * ((x + 1.0) / 2.0)).sin();
                let sn2 = (asr * ((-x + 1.0) / 2.0)).sin();
                w
                    * (((sn1 * hk - hs) / (1.0 - sn1 * sn1)).exp()
                        + ((sn2 * hk - hs) / (1.0 - sn2 * sn2)).exp())
            })
            .sum::<f64>()
            * asr
            / (TWOPI * 2.0)
            + mvphi(-h) * mvphi(-k)
    } else {
        //       ELSE
        //          IF ( R .LT. 0 ) THEN
        //             K = -K
        //             HK = -HK
        //          ENDIF
        if r < 0.0 {
            k = -k;
            hk = -hk;
        }
        //          IF ( ABS(R) .LT. 1 ) THEN
        let bvn = if r.abs() < 1.0 {
            //             AS = ( 1 - R )*( 1 + R )
            let r#as = (1.0 - r) * (1.0 + r);
            //             A = SQRT(AS)
            let a = r#as.sqrt();
            //             BS = ( H - K )**2
            let bs = (h - k) * (h - k);
            //             C = ( 4 - HK )/8
            //             D = ( 12 - HK )/16
            let c = (4.0 - hk) / 8.0;
            let d = (12.0 - hk) / 16.0;
            //             BVN = A*EXP( -(BS/AS + HK)/2 )
            //      +             *( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 )
            let mut bvn = a
                * (-(bs / r#as + hk) / 2.0).exp()
                * (1.0 - c * (bs - r#as) * (1.0 - d * bs / 5.0) / 3.0 + c * d * r#as * r#as / 5.0);
            //             IF ( HK .GT. -160 ) THEN
            if hk > -160.0 {
                //                B = SQRT(BS)
                let b = bs.sqrt();
                //                BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*MVPHI(-B/A)*B
                //      +                    *( 1 - C*BS*( 1 - D*BS/5 )/3 )
                bvn -= (-hk / 2.0).exp()
                    * TWOPI.sqrt()
                    * mvphi(-b / a)
                    * b
                    * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0);
                //             ENDIF
            }
            //             A = A/2
            let a = a / 2.0;
            //             DO I = 1, LG
            //                XS = ( A*(X(I,NG)+1) )**2
            //                RS = SQRT( 1 - XS )
            //                BVN = BVN + A*W(I,NG)*
            //      +              ( EXP( -BS/(2*XS) - HK/(1+RS) )/RS
            //      +              - EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) )
            //                XS = AS*(-X(I,NG)+1)**2/4
            //                RS = SQRT( 1 - XS )
            //                BVN = BVN + A*W(I,NG)*EXP( -(BS/XS + HK)/2 )
            //      +                    *( EXP( -HK*(1-RS)/(2*(1+RS)) )/RS
            //      +                       - ( 1 + C*XS*( 1 + D*XS ) ) )
            //             END DO
            -(bvn
                + (if r.abs() < 0.3 {
                        X0.iter().zip(W0.iter())
                    } else if r.abs() < 0.75 {
                        X1.iter().zip(W1.iter())
                    } else {
                        X2.iter().zip(W2.iter())
                    }).map(|(x,w)| {
                        let xs = (a * (x + 1.0)) * (a * (x + 1.0));
                        let rs = (1.0 - xs).sqrt();
                        let bvn = a
                            * w
                            * ((-bs / (2.0 * xs) - hk / (1.0 + rs)).exp() / rs
                                - (-(bs / xs + hk) / 2.0).exp() * (1.0 + c * xs * (1.0 + d * xs)));
                        let xs =
                            r#as * (-x + 1.0) * (-x + 1.0) / 4.0;
                        let rs = (1.0 - xs).sqrt();
                        bvn + a
                            * w
                            * (-(bs / xs + hk) / 2.0).exp()
                            * ((-hk * (1.0 - rs) / (2.0 * (1.0 + rs))).exp() / rs
                                - (1.0 + c * xs * (1.0 + d * xs)))
                    })
                    .sum::<f64>())
                / TWOPI
        //             BVN = -BVN/TWOPI

        //          ENDIF
        } else {
            0.0
        };
        //          IF ( R .GT. 0 ) BVN =  BVN + MVPHI( -MAX( H, K ) )
        //          IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, MVPHI(-H) - MVPHI(-K) )
        if r > 0.0 {
            bvn + mvphi(-(h.max(k)))
        } else {
            -bvn + 0f64.max(mvphi(-h) - mvphi(-k))
        }
        //       ENDIF
    }
    //       MVBVU = BVN
    //       END
}


fn mvphi(z: f64) -> f64 {
    //       INTEGER I, IM
    //       DOUBLE PRECISION A(0:43), BM, B, BP, P, RTWO, T, XA, Z
    //       PARAMETER( RTWO = 1.414213562373095048801688724209D0, IM = 24 )
    const IM: i32 = 24;
    const RTWO: f64 = 1.414213562373095048801688724209;
    //       SAVE A
    //       DATA ( A(I), I = 0, 43 )/
    //      &    6.10143081923200417926465815756D-1,
    //      &   -4.34841272712577471828182820888D-1,
    //      &    1.76351193643605501125840298123D-1,
    //      &   -6.0710795609249414860051215825D-2,
    //      &    1.7712068995694114486147141191D-2,
    //      &   -4.321119385567293818599864968D-3,
    //      &    8.54216676887098678819832055D-4,
    //      &   -1.27155090609162742628893940D-4,
    //      &    1.1248167243671189468847072D-5, 3.13063885421820972630152D-7,
    //      &   -2.70988068537762022009086D-7, 3.0737622701407688440959D-8,
    //      &    2.515620384817622937314D-9, -1.028929921320319127590D-9,
    //      &    2.9944052119949939363D-11, 2.6051789687266936290D-11,
    //      &   -2.634839924171969386D-12, -6.43404509890636443D-13,
    //      &    1.12457401801663447D-13, 1.7281533389986098D-14,
    //      &   -4.264101694942375D-15, -5.45371977880191D-16,
    //      &    1.58697607761671D-16, 2.0899837844334D-17,
    //      &   -5.900526869409D-18, -9.41893387554D-19, 2.14977356470D-19,
    //      &    4.6660985008D-20, -7.243011862D-21, -2.387966824D-21,
    //      &    1.91177535D-22, 1.20482568D-22, -6.72377D-25, -5.747997D-24,
    //      &   -4.28493D-25, 2.44856D-25, 4.3793D-26, -8.151D-27, -3.089D-27,
    //      &    9.3D-29, 1.74D-28, 1.6D-29, -8.0D-30, -2.0D-30 /
    // *
    //       XA = ABS(Z)/RTWO
    let xa = z.abs() / RTWO;
    //       IF ( XA .GT. 100 ) THEN
    let p = if xa >= 100.0 {
        0.0
    //          P = 0
    //       ELSE
    } else {
        //          T = ( 8*XA - 30 ) / ( 4*XA + 15 )
        let t = (8.0 * xa - 30.0) / (4.0 * xa + 15.0);
        let mut bm = 0.0;
        let mut b = 0.0;
        let mut bp = 0.0;
        //          BM = 0
        //          B  = 0
        //          DO I = IM, 0, -1
        //             BP = B
        //             B  = BM
        //             BM = T*B - BP  + A(I)
        //          END DO
        for i in (0..=IM).rev() {
            bp = b;
            b = bm;
            bm = t * b - bp + A[i as usize]
        }
        //          P = EXP( -XA*XA )*( BM - BP )/4
        (-xa * xa).exp() * (bm - bp) / 4.0
        //       END IF
    };
    //       IF ( Z .GT. 0 ) P = 1 - P
    //       MVPHI = P
    if z > 0.0 {
        1.0 - p
    } else {
        p
    }
    //       END
}

src/param.rs

pub static SW0: [f64; 3] =
        [
            0.1713244923791705,
            0.3607615730481384,
            0.4679139345726904,
        ];
pub static SW1:[f64;6]=
        [
            0.04717533638651177,
            0.1069393259953183,
            0.1600783285433464,
            0.2031674267230659,
            0.2334925365383547,
            0.2491470458134029,
        ];
pub static SW2:[f64;10]=
        [
            0.01761400713915212,
            0.04060142980038694,
            0.06267204833410905,
            0.08327674157670475,
            0.1019301198172404,
            0.1181945319615184,
            0.1316886384491766,
            0.1420961093183821,
            0.1491729864726037,
            0.1527533871307259,
        ];
pub static SX0: [f64; 3] =
        [
            -0.9324695142031522,
            -0.6612093864662647,
            -0.238619186083197,
        ];
pub static SX1:[f64; 6] = [
            -0.9815606342467191,
            -0.904117256370475,
            -0.769902674194305,
            -0.5873179542866171,
            -0.3678314989981802,
            -0.1252334085114692,
        ];
pub static SX2:[f64;10]=[
            -0.9931285991850949,
            -0.9639719272779138,
            -0.912234428251326,
            -0.8391169718222188,
            -0.7463319064601508,
            -0.636053680726515,
            -0.5108670019508271,
            -0.3737060887154196,
            -0.2277858511416451,
            -0.07652652113349732,
        ];

pub static SA: [f64; 44] = [
    6.10143081923200417926465815756e-1,
    -4.34841272712577471828182820888e-1,
    1.76351193643605501125840298123e-1,
    -6.0710795609249414860051215825e-2,
    1.7712068995694114486147141191e-2,
    -4.321119385567293818599864968e-3,
    8.54216676887098678819832055e-4,
    -1.27155090609162742628893940e-4,
    1.1248167243671189468847072e-5,
    3.13063885421820972630152e-7,
    -2.70988068537762022009086e-7,
    3.0737622701407688440959e-8,
    2.515620384817622937314e-9,
    -1.028929921320319127590e-9,
    2.9944052119949939363e-11,
    2.6051789687266936290e-11,
    -2.634839924171969386e-12,
    -6.43404509890636443e-13,
    1.12457401801663447e-13,
    1.7281533389986098e-14,
    -4.264101694942375e-15,
    -5.45371977880191e-16,
    1.58697607761671e-16,
    2.0899837844334e-17,
    -5.900526869409e-18,
    -9.41893387554e-19,
    2.14977356470e-19,
    4.6660985008e-20,
    -7.243011862e-21,
    -2.387966824e-21,
    1.91177535e-22,
    1.20482568e-22,
    -6.72377e-25,
    -5.747997e-24,
    -4.28493e-25,
    2.44856e-25,
    4.3793e-26,
    -8.151e-27,
    -3.089e-27,
    9.3e-29,
    1.74e-28,
    1.6e-29,
    -8.0e-30,
    -2.0e-30,
];

pub fn mvbvn(mut li: f64, mut lj: f64, mut ui: f64, mut uj: f64, mut corr: f64, W0:&[f64;3],W1:&[f64;6],W2:&[f64;10],X0:&[f64;3],X1:&[f64;6],X2:&[f64;10],A:&[f64;44]) -> f64 {
    if li >= ui || lj >= uj {
        0.0
    } else {
        if li.abs() > ui.abs() {
            (li, ui, corr) = (-ui, -li, -corr)
        }
        if lj.abs() > uj.abs() {
            (lj, uj, corr) = (-uj, -lj, -corr)
        }
        //         if li==f64::NEG_INFINITY {MVPHI(uj)-MVPHI(lj)}
        mvbvu(li, lj, corr,W0,W1,W2,X0,X1,X2,A)
            - if uj != f64::INFINITY {
                mvbvu(li, uj, corr,W0,W1,W2,X0,X1,X2,A)
            } else {
                0.0
            }
            - if ui != f64::INFINITY {
                mvbvu(ui, lj, corr,W0,W1,W2,X0,X1,X2,A)
            } else {
                0.0
            }
            + if ui + uj != f64::INFINITY {
                mvbvu(ui, uj, corr,W0,W1,W2,X0,X1,X2,A)
            } else {
                0.0
            }
    }
}

fn mvbvu(sh: f64, sk: f64, r: f64, W0:&[f64;3],W1:&[f64;6],W2:&[f64;10],X0:&[f64;3],X1:&[f64;6],X2:&[f64;10],A:&[f64;44]) -> f64 {
    //      DOUBLE PRECISION BVN, SH, SK, R, ZERO, TWOPI
    //     let mut bvn=0f64;
    //      INTEGER I, LG, NG
    //      PARAMETER ( ZERO = 0, TWOPI = 6.283185307179586D0 )
    const TWOPI: f64 = 6.283185307179586;
    //      DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS
    //      DOUBLE PRECISION MVPHI, SN, ASR, H, K, BS, HS, HK
    //      SAVE X, W
    // *     Gauss Legendre Points and Weights, N =  6
    //       DATA ( W(I,1), X(I,1), I = 1, 3 ) /
    //      *  0.1713244923791705D+00,-0.9324695142031522D+00,
    //      *  0.3607615730481384D+00,-0.6612093864662647D+00,
    //      *  0.4679139345726904D+00,-0.2386191860831970D+00/
    // *     Gauss Legendre Points and Weights, N = 12
    //       DATA ( W(I,2), X(I,2), I = 1, 6 ) /
    //      *  0.4717533638651177D-01,-0.9815606342467191D+00,
    //      *  0.1069393259953183D+00,-0.9041172563704750D+00,
    //      *  0.1600783285433464D+00,-0.7699026741943050D+00,
    //      *  0.2031674267230659D+00,-0.5873179542866171D+00,
    //      *  0.2334925365383547D+00,-0.3678314989981802D+00,
    //      *  0.2491470458134029D+00,-0.1252334085114692D+00/
    // *     Gauss Legendre Points and Weights, N = 20
    //       DATA ( W(I,3), X(I,3), I = 1, 10 ) /
    //      *  0.1761400713915212D-01,-0.9931285991850949D+00,
    //      *  0.4060142980038694D-01,-0.9639719272779138D+00,
    //      *  0.6267204833410906D-01,-0.9122344282513259D+00,
    //      *  0.8327674157670475D-01,-0.8391169718222188D+00,
    //      *  0.1019301198172404D+00,-0.7463319064601508D+00,
    //      *  0.1181945319615184D+00,-0.6360536807265150D+00,
    //      *  0.1316886384491766D+00,-0.5108670019508271D+00,
    //      *  0.1420961093183821D+00,-0.3737060887154196D+00,
    //      *  0.1491729864726037D+00,-0.2277858511416451D+00,
    //      *  0.1527533871307259D+00,-0.7652652113349733D-01/

    //       IF ( ABS(R) .LT. 0.3 ) THEN
    //          NG = 1
    //          LG = 3
    //       ELSE IF ( ABS(R) .LT. 0.75 ) THEN
    //          NG = 2
    //          LG = 6
    //       ELSE
    //          NG = 3
    //          LG = 10
    //       ENDIF
    //       H = SH
    //       K = SK
    let h = sh;
    let mut k = sk;
    //       HK = H*K
    let mut hk = h * k;
    //       BVN = 0
    //       IF ( ABS(R) .LT. 0.925 ) THEN
    if r.abs() < 0.925 {
        //          HS = ( H*H + K*K )/2
        let hs = (h * h + k * k) / 2.0;
        //          ASR = ASIN(R)
        let asr = r.asin();
        //          DO I = 1, LG
        //         for i in 0..lg{
        // //             SN = SIN(ASR*( X(I,NG)+1 )/2)
        //             sn=(asr*(( X[ng][i]+1f64 )/2f64)).sin();
        // //             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
        //             bvn+=W[ng][i]*((sn*hk-hs)/(1-sn*sn)).exp()
        // //             SN = SIN(ASR*(-X(I,NG)+1 )/2)
        // //             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
        // //          END DO
        //         }
        //          BVN = BVN*ASR/(2*TWOPI) + MVPHI(-H)*MVPHI(-K)
        (if r.abs() < 0.3 {
            X0.iter().zip(W0.iter())
        } else if r.abs() < 0.75 {
            X1.iter().zip(W1.iter())
        } else {
            X2.iter().zip(W2.iter())
        }).map(|(x,w)| {
                let sn1 = (asr * ((x + 1.0) / 2.0)).sin();
                let sn2 = (asr * ((-x + 1.0) / 2.0)).sin();
                w
                    * (((sn1 * hk - hs) / (1.0 - sn1 * sn1)).exp()
                        + ((sn2 * hk - hs) / (1.0 - sn2 * sn2)).exp())
            })
            .sum::<f64>()
            * asr
            / (TWOPI * 2.0)
            + mvphi(-h,A) * mvphi(-k,A)
    } else {
        //       ELSE
        //          IF ( R .LT. 0 ) THEN
        //             K = -K
        //             HK = -HK
        //          ENDIF
        if r < 0.0 {
            k = -k;
            hk = -hk;
        }
        //          IF ( ABS(R) .LT. 1 ) THEN
        let bvn = if r.abs() < 1.0 {
            //             AS = ( 1 - R )*( 1 + R )
            let r#as = (1.0 - r) * (1.0 + r);
            //             A = SQRT(AS)
            let a = r#as.sqrt();
            //             BS = ( H - K )**2
            let bs = (h - k) * (h - k);
            //             C = ( 4 - HK )/8
            //             D = ( 12 - HK )/16
            let c = (4.0 - hk) / 8.0;
            let d = (12.0 - hk) / 16.0;
            //             BVN = A*EXP( -(BS/AS + HK)/2 )
            //      +             *( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 )
            let mut bvn = a
                * (-(bs / r#as + hk) / 2.0).exp()
                * (1.0 - c * (bs - r#as) * (1.0 - d * bs / 5.0) / 3.0 + c * d * r#as * r#as / 5.0);
            //             IF ( HK .GT. -160 ) THEN
            if hk > -160.0 {
                //                B = SQRT(BS)
                let b = bs.sqrt();
                //                BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*MVPHI(-B/A)*B
                //      +                    *( 1 - C*BS*( 1 - D*BS/5 )/3 )
                bvn -= (-hk / 2.0).exp()
                    * TWOPI.sqrt()
                    * mvphi(-b / a,A)
                    * b
                    * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0);
                //             ENDIF
            }
            //             A = A/2
            let a = a / 2.0;
            //             DO I = 1, LG
            //                XS = ( A*(X(I,NG)+1) )**2
            //                RS = SQRT( 1 - XS )
            //                BVN = BVN + A*W(I,NG)*
            //      +              ( EXP( -BS/(2*XS) - HK/(1+RS) )/RS
            //      +              - EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) )
            //                XS = AS*(-X(I,NG)+1)**2/4
            //                RS = SQRT( 1 - XS )
            //                BVN = BVN + A*W(I,NG)*EXP( -(BS/XS + HK)/2 )
            //      +                    *( EXP( -HK*(1-RS)/(2*(1+RS)) )/RS
            //      +                       - ( 1 + C*XS*( 1 + D*XS ) ) )
            //             END DO
            -(bvn
                + (if r.abs() < 0.3 {
                        X0.iter().zip(W0.iter())
                    } else if r.abs() < 0.75 {
                        X1.iter().zip(W1.iter())
                    } else {
                        X2.iter().zip(W2.iter())
                    }).map(|(x,w)| {
                        let xs = (a * (x + 1.0)) * (a * (x + 1.0));
                        let rs = (1.0 - xs).sqrt();
                        let bvn = a
                            * w
                            * ((-bs / (2.0 * xs) - hk / (1.0 + rs)).exp() / rs
                                - (-(bs / xs + hk) / 2.0).exp() * (1.0 + c * xs * (1.0 + d * xs)));
                        let xs =
                            r#as * (-x + 1.0) * (-x + 1.0) / 4.0;
                        let rs = (1.0 - xs).sqrt();
                        bvn + a
                            * w
                            * (-(bs / xs + hk) / 2.0).exp()
                            * ((-hk * (1.0 - rs) / (2.0 * (1.0 + rs))).exp() / rs
                                - (1.0 + c * xs * (1.0 + d * xs)))
                    })
                    .sum::<f64>())
                / TWOPI
        //             BVN = -BVN/TWOPI

        //          ENDIF
        } else {
            0.0
        };
        //          IF ( R .GT. 0 ) BVN =  BVN + MVPHI( -MAX( H, K ) )
        //          IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, MVPHI(-H) - MVPHI(-K) )
        if r > 0.0 {
            bvn + mvphi(-(h.max(k)),A)
        } else {
            -bvn + 0f64.max(mvphi(-h,A) - mvphi(-k,A))
        }
        //       ENDIF
    }
    //       MVBVU = BVN
    //       END
}


fn mvphi(z: f64, A:&[f64;44]) -> f64 {
    //       INTEGER I, IM
    //       DOUBLE PRECISION A(0:43), BM, B, BP, P, RTWO, T, XA, Z
    //       PARAMETER( RTWO = 1.414213562373095048801688724209D0, IM = 24 )
    const IM: i32 = 24;
    const RTWO: f64 = 1.414213562373095048801688724209;
    //       SAVE A
    //       DATA ( A(I), I = 0, 43 )/
    //      &    6.10143081923200417926465815756D-1,
    //      &   -4.34841272712577471828182820888D-1,
    //      &    1.76351193643605501125840298123D-1,
    //      &   -6.0710795609249414860051215825D-2,
    //      &    1.7712068995694114486147141191D-2,
    //      &   -4.321119385567293818599864968D-3,
    //      &    8.54216676887098678819832055D-4,
    //      &   -1.27155090609162742628893940D-4,
    //      &    1.1248167243671189468847072D-5, 3.13063885421820972630152D-7,
    //      &   -2.70988068537762022009086D-7, 3.0737622701407688440959D-8,
    //      &    2.515620384817622937314D-9, -1.028929921320319127590D-9,
    //      &    2.9944052119949939363D-11, 2.6051789687266936290D-11,
    //      &   -2.634839924171969386D-12, -6.43404509890636443D-13,
    //      &    1.12457401801663447D-13, 1.7281533389986098D-14,
    //      &   -4.264101694942375D-15, -5.45371977880191D-16,
    //      &    1.58697607761671D-16, 2.0899837844334D-17,
    //      &   -5.900526869409D-18, -9.41893387554D-19, 2.14977356470D-19,
    //      &    4.6660985008D-20, -7.243011862D-21, -2.387966824D-21,
    //      &    1.91177535D-22, 1.20482568D-22, -6.72377D-25, -5.747997D-24,
    //      &   -4.28493D-25, 2.44856D-25, 4.3793D-26, -8.151D-27, -3.089D-27,
    //      &    9.3D-29, 1.74D-28, 1.6D-29, -8.0D-30, -2.0D-30 /
    // *
    //       XA = ABS(Z)/RTWO
    let xa = z.abs() / RTWO;
    //       IF ( XA .GT. 100 ) THEN
    let p = if xa >= 100.0 {
        0.0
    //          P = 0
    //       ELSE
    } else {
        //          T = ( 8*XA - 30 ) / ( 4*XA + 15 )
        let t = (8.0 * xa - 30.0) / (4.0 * xa + 15.0);
        let mut bm = 0.0;
        let mut b = 0.0;
        let mut bp = 0.0;
        //          BM = 0
        //          B  = 0
        //          DO I = IM, 0, -1
        //             BP = B
        //             B  = BM
        //             BM = T*B - BP  + A(I)
        //          END DO
        for i in (0..=IM).rev() {
            bp = b;
            b = bm;
            bm = t * b - bp + A[i as usize]
        }
        //          P = EXP( -XA*XA )*( BM - BP )/4
        (-xa * xa).exp() * (bm - bp) / 4.0
        //       END IF
    };
    //       IF ( Z .GT. 0 ) P = 1 - P
    //       MVPHI = P
    if z > 0.0 {
        1.0 - p
    } else {
        p
    }
    //       END
}

benchmark

use criterion::{black_box, criterion_group, criterion_main, Criterion};


pub fn bench_orig(c: &mut Criterion) {
    use benchmark_static::orig::*;
    c.bench_function("orig", |b| b.iter(|| black_box([[-1f64,-1.,1.,1.];100].into_iter().zip((0..100).map(|x|(x*2-99) as f64/100.))).map(|x|mvbvn(x.0[0],x.0[1],x.0[2],x.0[3],x.1)).sum::<f64>()));
}
pub fn bench_global(c: &mut Criterion) {
    use benchmark_static::r#static::*;
    c.bench_function("global", |b| b.iter(|| black_box([[-1f64,-1.,1.,1.];100].into_iter().zip((0..100).map(|x|(x*2-99) as f64/100.))).map(|x|mvbvn(x.0[0],x.0[1],x.0[2],x.0[3],x.1)).sum::<f64>()));
}
pub fn bench_param(c: &mut Criterion) {
    use benchmark_static::param::*;
    let (W0,W1,W2,X0,X1,X2,A)=(&SW0,&SW1,&SW2,&SX0,&SX1,&SX2,&SA);
    c.bench_function("param", |b| b.iter(|| black_box([[-1f64,-1.,1.,1.];100].into_iter().zip((0..100).map(|x|(x*2-99) as f64/100.))).map(|x|mvbvn(x.0[0],x.0[1],x.0[2],x.0[3],x.1,W0,W1,W2,X0,X1,X2,A)).sum::<f64>()));
}


criterion_group!(orig, bench_orig);
criterion_group!(global, bench_global);
criterion_group!(param, bench_param);

criterion_main!(orig, global, param);

differences between static.rs and param.rs:

$ diff static.rs param.rs 
1c1
< static W0: [f64; 3] =
---
> pub static SW0: [f64; 3] =
7c7
< static W1:[f64;6]=
---
> pub static SW1:[f64;6]=
16c16
< static W2:[f64;10]=
---
> pub static SW2:[f64;10]=
29c29
< static X0: [f64; 3] =
---
> pub static SX0: [f64; 3] =
35c35
< static X1:[f64; 6] = [
---
> pub static SX1:[f64; 6] = [
43c43
< static X2:[f64;10]=[
---
> pub static SX2:[f64;10]=[
56c56
< static A: [f64; 44] = [
---
> pub static SA: [f64; 44] = [
103c103
< pub fn mvbvn(mut li: f64, mut lj: f64, mut ui: f64, mut uj: f64, mut corr: f64) -> f64 {
---
> pub fn mvbvn(mut li: f64, mut lj: f64, mut ui: f64, mut uj: f64, mut corr: f64, W0:&[f64;3],W1:&[f64;6],W2:&[f64;10],X0:&[f64;3],X1:&[f64;6],X2:&[f64;10],A:&[f64;44]) -> f64 {
114c114
<         mvbvu(li, lj, corr)
---
>         mvbvu(li, lj, corr,W0,W1,W2,X0,X1,X2,A)
116c116
<                 mvbvu(li, uj, corr)
---
>                 mvbvu(li, uj, corr,W0,W1,W2,X0,X1,X2,A)
121c121
<                 mvbvu(ui, lj, corr)
---
>                 mvbvu(ui, lj, corr,W0,W1,W2,X0,X1,X2,A)
126c126
<                 mvbvu(ui, uj, corr)
---
>                 mvbvu(ui, uj, corr,W0,W1,W2,X0,X1,X2,A)
133c133
< fn mvbvu(sh: f64, sk: f64, r: f64) -> f64 {
---
> fn mvbvu(sh: f64, sk: f64, r: f64, W0:&[f64;3],W1:&[f64;6],W2:&[f64;10],X0:&[f64;3],X1:&[f64;6],X2:&[f64;10],A:&[f64;44]) -> f64 {
218c218
<             + mvphi(-h) * mvphi(-k)
---
>             + mvphi(-h,A) * mvphi(-k,A)
254c254
<                     * mvphi(-b / a)
---
>                     * mvphi(-b / a,A)
307c307
<             bvn + mvphi(-(h.max(k)))
---
>             bvn + mvphi(-(h.max(k)),A)
309c309
<             -bvn + 0f64.max(mvphi(-h) - mvphi(-k))
---
>             -bvn + 0f64.max(mvphi(-h,A) - mvphi(-k,A))
318c318
< fn mvphi(z: f64) -> f64 {
---
> fn mvphi(z: f64, A:&[f64;44]) -> f64 {

Results:

    Finished bench [optimized] target(s) in 0.04s
     Running unittests src/lib.rs (/me/ben/target/release/deps/benchmark_static-56b725c82118e0da)

running 0 tests

test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s

     Running benches/bench.rs (/me/ben/target/release/deps/bench-7b5d4883affc7873)
orig                    time:   [187.64 µs 187.71 µs 187.83 µs]
                        change: [-0.3591% -0.2672% -0.1605%] (p = 0.00 < 0.05)
                        Change within noise threshold.
Found 13 outliers among 100 measurements (13.00%)
  2 (2.00%) high mild
  11 (11.00%) high severe

global                  time:   [188.64 µs 188.78 µs 189.09 µs]
                        change: [+0.0078% +0.1025% +0.2603%] (p = 0.14 > 0.05)
                        No change in performance detected.
Found 10 outliers among 100 measurements (10.00%)
  1 (1.00%) low severe
  1 (1.00%) low mild
  3 (3.00%) high mild
  5 (5.00%) high severe

param                   time:   [190.78 µs 190.80 µs 190.81 µs]
                        change: [-0.6868% -0.6719% -0.6578%] (p = 0.00 < 0.05)
                        Change within noise threshold.
Found 6 outliers among 100 measurements (6.00%)
  1 (1.00%) low mild
  4 (4.00%) high mild
  1 (1.00%) high severe

[farnz]

@Neutron3529

@farnz

Rust has been able to inline in both cases

Will, this is not true, since we could have REALLY COMPLEX function that could not be inlined.

If you are simply going to accuse me of saying something that's not true, on the basis that I could have written completely different code and then Rust wouldn't have behaved the same way, then there's no basis for communication here. That was a statement about what I was seeing with my example, not about all code - there is at least one case where Rust was able to do better with a data table passed as &[…; N] rather than available globally as static […; N].

You've then gone on a different direction - you've shown that where Rust can't exploit inlining, there's no difference between const and static. But my claim is different; mine is that there is no benefit to static DATA: DataType versus a parameter &DataType, and that in some cases (but not others), &DataType will enable better optimization.

@toc

But yes I have looked at init-token and considered the approach. In this case (coding a chess ai) there is a big public init which calls a number of private ones and it will all definitely blow up if you haven't called init. Most paths towards calling that API force you through some bottleneck which will call init (for instance, creating a new UCI interface). I considered it important to be able to call individual free functions fluently; simple_calc(x, y) looks the same as table_based_calc(x, y, z). And a part of me is more confident that this will optimize correctly:

I've been thinking about this - are the tables you create dependent on properties of the running system, or are you generating tables on startup because it's non-trivial to make the initializer const-evaluated?

If the latter (the tables are the same in every binary), then an easier feature to get right would be making it simpler to run code at compile time that outputs the tables you need. This works for your example of a sin lookup table; if the reason for the table to exist is to have precomputed values of f(x) for interesting values of x, there's no requirement to do the computation at startup instead of at compile time.

And in theory, you could write yourself a build.rs that used quote - Rust to write out the tables. But that's a lot harder than writing an init_table function and calling it at runtime - would making it easier to run code generators at compile time reduce the need for static initializers?

[kpreid]

@farnz

And in theory, you could write yourself a build.rs that used quote - Rust to write out the tables.

Tool to make this easier than quote: uneval - Rust

[toc]

@farnz

I've been thinking about this - are the tables you create dependent on properties of the running system, or are you generating tables on startup because it's non-trivial to make the initializer const-evaluated?

Some of these tables will be const at some point in the future. Some of the tables are big (enough to inflate the binary size), and potentially would inflate compile times as const eval is slow (though that should get better). These tables get allocated, computed, and Vec::leaked.

For this latter case, when/as const becomes more powerful, my purposes would be satisfied by allowing static initializers which can only call const functions. That would not work for me at present, but would be sufficient at some point.