ellip_rayon/
lib.rs

1/*
2 * Ellip is licensed under The 3-Clause BSD, see LICENSE.
3 * Copyright 2025 Sira Pornsiriprasert <code@psira.me>
4 */
5
6//! # Ellip Rayon
7//! Parallelized elliptic integral computation for Rust based on [Ellip](https://github.com/p-sira/ellip).
8//!
9//! ## Installation
10//! ```shell
11//! cargo add ellip-rayon
12//! ```
13//!
14//! ## Machine-specific Threshold
15//! Ellip Rayon employs parallelization if the length of the arguments exceesds certain thresholds. These thresholds depend on the core count, cache size, and architecture. For the most efficiency, these thresholds should be tuned on the target machine.
16//!
17//! 1. Generating benchmark data
18//!
19//! ```shell
20//! cargo bench
21//! ```
22//! The process should take about 30-40 minutes to complete, and the file `benches/par_threshold.md` will be created. This reports the threshold for each function.
23//!
24//! 2. Using the generated threshold
25//!
26//! ```shell
27//! cargo run --example generate_threshold_code
28//! ```
29//! This script automatically replaces the thresholds in the source code.
30//!
31//! 3. Adding locally compiled library
32//!
33//! From your working directory, run
34//! ```shell
35//! cargo add --path path/to/your/ellip-rayon
36//! ```
37
38use ellip::*;
39use itertools::izip;
40use rayon::iter::{IndexedParallelIterator, IntoParallelRefIterator, ParallelIterator};
41
42macro_rules! par_zip {
43    ($a:expr) => {
44        $a.par_iter()
45    };
46    ($a:expr, $b:expr) => {
47        $a.par_iter().zip($b.par_iter())
48    };
49    ($a:expr, $b:expr, $c:expr) => {
50        $a.par_iter()
51            .zip($b.par_iter())
52            .zip($c.par_iter())
53            .map(|((a, b), c)| (a, b, c))
54    };
55    ($a:expr, $b:expr, $c:expr, $d:expr) => {
56        $a.par_iter()
57            .zip($b.par_iter())
58            .zip($c.par_iter())
59            .zip($d.par_iter())
60            .map(|(((a, b), c), d)| (a, b, c, d))
61    };
62}
63
64macro_rules! impl_par {
65    (@inner, $fn:ident, 2) => {
66        |(&a, &b)| ellip::$fn(a, b)
67    };
68    (@inner, $fn:ident, 3) => {
69        |(&a, &b, &c)| ellip::$fn(a, b, c)
70    };
71    (@inner, $fn:ident, 4) => {
72        |(&a, &b, &c, &d)| ellip::$fn(a, b, c, d)
73    };
74    ($fn:ident, [$arg:ident], 1) => {
75        #[doc=concat!["Computes [", stringify!($fn), "](ellip::", stringify!($fn), ") in parallel."]]
76        pub fn $fn($arg: &[f64]) -> Result<Vec<f64>, StrErr> {
77            $arg.iter().map(|&a| ellip::$fn(a)).collect()
78        }
79    };
80    ($fn:ident, [$arg:ident], 1, $threshold:expr) => {
81        #[doc=concat!["Computes [", stringify!($fn), "](ellip::", stringify!($fn), ") in parallel."]]
82        pub fn $fn($arg: &[f64]) -> Result<Vec<f64>, StrErr> {
83            if $arg.len() < $threshold {
84                $arg.iter().map(|&a| ellip::$fn(a)).collect()
85            } else {
86                $arg.par_iter().map(|&a| ellip::$fn(a)).collect()
87            }
88        }
89    };
90    ($fn:ident, [$first:ident, $($args:ident),*], $n_arg:tt) => {
91        #[doc=concat!["Computes [", stringify!($fn), "](ellip::", stringify!($fn), ") in parallel."]]
92        pub fn $fn($first: &[f64], $($args: &[f64],)*) -> Result<Vec<f64>, StrErr> {
93            $(
94                if $first.len() != $args.len() {
95                    return Err(concat![stringify!($fn), ": All arguments must have the same length."]);
96                }
97            )*
98            izip!($first, $($args),*).map(impl_par!(@inner, $fn, $n_arg)).collect()
99        }
100    };
101    ($fn:ident, [$first:ident, $($args:ident),*], $n_arg:tt, $threshold:expr) => {
102        #[doc=concat!["Computes [", stringify!($fn), "](ellip::", stringify!($fn), ") in parallel."]]
103        pub fn $fn($first: &[f64], $($args: &[f64],)*) -> Result<Vec<f64>, StrErr> {
104            $(
105                if $first.len() != $args.len() {
106                    return Err(concat![stringify!($fn), ": All arguments must have the same length."]);
107                }
108            )*
109            if $first.len() < $threshold {
110                izip!($first, $($args),*).map(impl_par!(@inner, $fn, $n_arg)).collect()
111            } else {
112                par_zip!($first, $($args),*).map(impl_par!(@inner, $fn, $n_arg)).collect()
113            }
114        }
115    };
116}
117
118// {BEGIN_IMPL_PAR}
119// Generated threshold values from benchmark results
120impl_par!(ellipk, [m], 1, 1000);
121impl_par!(ellipe, [m], 1, 1500);
122impl_par!(ellipf, [phi, m], 2, 400);
123impl_par!(ellipeinc, [phi, m], 2, 300);
124impl_par!(ellippi, [n, m], 2, 200);
125impl_par!(ellippiinc, [phi, n, m], 3, 200);
126impl_par!(ellippiinc_bulirsch, [phi, n, m], 3, 300);
127impl_par!(ellipd, [m], 1, 600);
128impl_par!(ellipdinc, [phi, m], 2, 500);
129impl_par!(cel, [kc, p, a, b], 4, 600);
130impl_par!(cel1, [kc], 1, 1500);
131impl_par!(cel2, [kc, a, b], 3, 700);
132impl_par!(el1, [x, kc], 2, 600);
133impl_par!(el2, [x, kc, a, b], 4, 600);
134impl_par!(el3, [x, kc, p], 3, 500);
135impl_par!(elliprf, [x, y, z], 3, 600);
136impl_par!(elliprg, [x, y, z], 3, 500);
137impl_par!(elliprj, [x, y, z, p], 4, 300);
138impl_par!(elliprc, [x, y], 2, 800);
139impl_par!(elliprd, [x, y, z], 3, 500);
140impl_par!(jacobi_zeta, [phi, m], 2, 200);
141impl_par!(heuman_lambda, [phi, m], 2, 200);
142
143// {END_IMPL_PAR}
ellip_rayon/lib.rs

ellip_rayon/
lib.rs
