Skip to content

Commit

Permalink
rust: Add cmplx!
Browse files Browse the repository at this point in the history
  • Loading branch information
felixhekhorn committed Oct 23, 2023
1 parent 1fa54f3 commit 408b5be
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 197 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -54,14 +54,15 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {

#[cfg(test)]
mod tests {
use crate::cmplx;
use crate::{anomalous_dimensions::unpolarized::spacelike::as1::*, harmonics::cache::Cache};
use float_cmp::assert_approx_eq;
use num::complex::Complex;
const NF: u8 = 5;

#[test]
fn number_conservation() {
const N: Complex<f64> = Complex::<f64> { re: 1., im: 0. };
const N: Complex<f64> = cmplx![1., 0.];
let mut c = Cache::new(N);
let me = gamma_ns(&mut c, NF);
assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
Expand All @@ -70,7 +71,7 @@ mod tests {

#[test]
fn quark_momentum_conservation() {
const N: Complex<f64> = Complex::<f64> { re: 2., im: 0. };
const N: Complex<f64> = cmplx![2., 0.];
let mut c = Cache::new(N);
let me = gamma_ns(&mut c, NF) + gamma_gq(&mut c, NF);
assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
Expand All @@ -79,7 +80,7 @@ mod tests {

#[test]
fn gluon_momentum_conservation() {
const N: Complex<f64> = Complex::<f64> { re: 2., im: 0. };
const N: Complex<f64> = cmplx![2., 0.];
let mut c = Cache::new(N);
let me = gamma_qg(&mut c, NF) + gamma_gg(&mut c, NF);
assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
Expand All @@ -88,7 +89,7 @@ mod tests {

#[test]
fn gamma_qg_() {
const N: Complex<f64> = Complex::<f64> { re: 1., im: 0. };
const N: Complex<f64> = cmplx![1., 0.];
let mut c = Cache::new(N);
let me = gamma_qg(&mut c, NF);
assert_approx_eq!(f64, me.re, -20. / 3.0, ulps = 32);
Expand All @@ -97,7 +98,7 @@ mod tests {

#[test]
fn gamma_gq_() {
const N: Complex<f64> = Complex::<f64> { re: 0., im: 1. };
const N: Complex<f64> = cmplx![0., 1.];
let mut c = Cache::new(N);
let me = gamma_gq(&mut c, NF);
assert_approx_eq!(f64, me.re, 4. / 3.0, ulps = 32);
Expand All @@ -106,7 +107,7 @@ mod tests {

#[test]
fn gamma_gg_() {
const N: Complex<f64> = Complex::<f64> { re: 0., im: 1. };
const N: Complex<f64> = cmplx![0., 1.];
let mut c = Cache::new(N);
let me = gamma_gg(&mut c, NF);
assert_approx_eq!(f64, me.re, 5.195725159621, ulps = 32, epsilon = 1e-11);
Expand Down
248 changes: 57 additions & 191 deletions crates/ekore/src/harmonics/polygamma.rs
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
//! Tools to compute [polygamma functions](https://en.wikipedia.org/wiki/Polygamma_function).
use num::complex::Complex;
use num::{complex::Complex, Zero};
use std::f64::consts::PI;

#[cfg_attr(doc, katexit::katexit)]
Expand Down Expand Up @@ -72,7 +72,7 @@ pub fn cern_polygamma(Z: Complex<f64>, K: usize) -> Complex<f64> {
U = -U;
}
let mut V = U;
let mut H = Complex::<f64> { re: 0., im: 0. };
let mut H = Complex::zero();
if A < 15. {
H = 1. / V.powu(K1 as u32);
for _ in 1..(14 - A_as_int + 1) {
Expand Down Expand Up @@ -121,215 +121,81 @@ pub fn cern_polygamma(Z: Complex<f64>, K: usize) -> Complex<f64> {

#[cfg(test)]
mod tests {
use crate::cmplx;
use crate::harmonics::polygamma::cern_polygamma;
use float_cmp::assert_approx_eq;
use num::complex::Complex;

#[test]
fn fortran() {
const ZS: [Complex<f64>; 9] = [
Complex::<f64> { re: 1.0, im: 0. },
Complex::<f64> { re: 2.0, im: 0. },
Complex::<f64> { re: 3.0, im: 0. },
Complex::<f64> { re: 0., im: 1. },
Complex::<f64> { re: -1., im: 1. },
Complex::<f64> { re: -2., im: 1. },
Complex::<f64> { re: -1., im: 2. },
Complex::<f64> { re: -2., im: 2. },
Complex::<f64> { re: -3., im: 2. },
cmplx![1.0, 0.],
cmplx![2.0, 0.],
cmplx![3.0, 0.],
cmplx![0., 1.],
cmplx![-1., 1.],
cmplx![-2., 1.],
cmplx![-1., 2.],
cmplx![-2., 2.],
cmplx![-3., 2.],
];
const KS: [usize; 5] = [0, 1, 2, 3, 4];

const FORTRAN_REF: [[Complex<f64>; 9]; 5] = [
[
Complex::<f64> {
re: -0.5772156649015332,
im: 0.,
},
Complex::<f64> {
re: 0.42278433509846636,
im: 0.,
},
Complex::<f64> {
re: 0.9227843350984672,
im: 0.,
},
Complex::<f64> {
re: 0.09465032062247669,
im: 2.0766740474685816,
},
Complex::<f64> {
re: 0.5946503206224767,
im: 2.5766740474685816,
},
Complex::<f64> {
re: 0.9946503206224772,
im: 2.7766740474685814,
},
Complex::<f64> {
re: 0.9145915153739776,
im: 2.2208072826422303,
},
Complex::<f64> {
re: 1.1645915153739772,
im: 2.47080728264223,
},
Complex::<f64> {
re: 1.395360746143208,
im: 2.624653436488384,
},
cmplx![-0.5772156649015332, 0.],
cmplx![0.42278433509846636, 0.],
cmplx![0.9227843350984672, 0.],
cmplx![0.09465032062247669, 2.0766740474685816],
cmplx![0.5946503206224767, 2.5766740474685816],
cmplx![0.9946503206224772, 2.7766740474685814],
cmplx![0.9145915153739776, 2.2208072826422303],
cmplx![1.1645915153739772, 2.47080728264223],
cmplx![1.395360746143208, 2.624653436488384],
],
[
Complex::<f64> {
re: 1.6449340668482264,
im: 0.,
},
Complex::<f64> {
re: 0.6449340668482264,
im: 0.,
},
Complex::<f64> {
re: 0.3949340668482264,
im: 0.,
},
Complex::<f64> {
re: -0.5369999033772361,
im: -0.7942335427593189,
},
Complex::<f64> {
re: -0.5369999033772362,
im: -0.2942335427593189,
},
Complex::<f64> {
re: -0.4169999033772362,
im: -0.13423354275931887,
},
Complex::<f64> {
re: -0.24506883785905695,
im: -0.3178255501472297,
},
Complex::<f64> {
re: -0.24506883785905695,
im: -0.19282555014722969,
},
Complex::<f64> {
re: -0.21548303904248892,
im: -0.12181963298746643,
},
cmplx![1.6449340668482264, 0.],
cmplx![0.6449340668482264, 0.],
cmplx![0.3949340668482264, 0.],
cmplx![-0.5369999033772361, -0.7942335427593189],
cmplx![-0.5369999033772362, -0.2942335427593189],
cmplx![-0.4169999033772362, -0.13423354275931887],
cmplx![-0.24506883785905695, -0.3178255501472297],
cmplx![-0.24506883785905695, -0.19282555014722969],
cmplx![-0.21548303904248892, -0.12181963298746643],
],
[
Complex::<f64> {
re: -2.404113806319188,
im: 0.,
},
Complex::<f64> {
re: -0.40411380631918853,
im: 0.,
},
Complex::<f64> {
re: -0.15411380631918858,
im: 0.,
},
Complex::<f64> {
re: 0.3685529315879351,
im: -1.233347149654934,
},
Complex::<f64> {
re: -0.13144706841206477,
im: -0.7333471496549337,
},
Complex::<f64> {
re: -0.09944706841206462,
im: -0.5573471496549336,
},
Complex::<f64> {
re: 0.03902435405364951,
im: -0.15743252404131272,
},
Complex::<f64> {
re: -0.02347564594635048,
im: -0.09493252404131272,
},
Complex::<f64> {
re: -0.031668636387861625,
im: -0.053057239562477945,
},
cmplx![-2.404113806319188, 0.],
cmplx![-0.40411380631918853, 0.],
cmplx![-0.15411380631918858, 0.],
cmplx![0.3685529315879351, -1.233347149654934],
cmplx![-0.13144706841206477, -0.7333471496549337],
cmplx![-0.09944706841206462, -0.5573471496549336],
cmplx![0.03902435405364951, -0.15743252404131272],
cmplx![-0.02347564594635048, -0.09493252404131272],
cmplx![-0.031668636387861625, -0.053057239562477945],
],
[
Complex::<f64> {
re: 6.49393940226683,
im: 0.,
},
Complex::<f64> {
re: 0.49393940226682925,
im: 0.,
},
Complex::<f64> {
re: 0.11893940226682913,
im: 0.,
},
Complex::<f64> {
re: 4.4771255510465044,
im: -0.31728657866196064,
},
Complex::<f64> {
re: 2.9771255510464925,
im: -0.3172865786619599,
},
Complex::<f64> {
re: 2.909925551046492,
im: -0.08688657866195917,
},
Complex::<f64> {
re: 0.12301766661068443,
im: -0.05523068481179527,
},
Complex::<f64> {
re: 0.02926766661068438,
im: -0.055230684811795216,
},
Complex::<f64> {
re: 0.004268541930176011,
im: -0.03002148345329936,
},
cmplx![6.49393940226683, 0.],
cmplx![0.49393940226682925, 0.],
cmplx![0.11893940226682913, 0.],
cmplx![4.4771255510465044, -0.31728657866196064],
cmplx![2.9771255510464925, -0.3172865786619599],
cmplx![2.909925551046492, -0.08688657866195917],
cmplx![0.12301766661068443, -0.05523068481179527],
cmplx![0.02926766661068438, -0.055230684811795216],
cmplx![0.004268541930176011, -0.03002148345329936],
],
[
Complex::<f64> {
re: -24.88626612344089,
im: 0.,
},
Complex::<f64> {
re: -0.8862661234408784,
im: 0.,
},
Complex::<f64> {
re: -0.13626612344087824,
im: 0.,
},
Complex::<f64> {
re: 3.2795081690440493,
im: 21.41938794863803,
},
Complex::<f64> {
re: 0.2795081690440445,
im: 18.419387948637894,
},
Complex::<f64> {
re: -0.012331830955960252,
im: 18.734267948637896,
},
Complex::<f64> {
re: 0.14223316576854003,
im: 0.10023607930398608,
},
Complex::<f64> {
re: 0.04848316576854002,
im: 0.006486079303986134,
},
Complex::<f64> {
re: 0.009893695996688708,
im: 0.014372034600746361,
},
cmplx![-24.88626612344089, 0.],
cmplx![-0.8862661234408784, 0.],
cmplx![-0.13626612344087824, 0.],
cmplx![3.2795081690440493, 21.41938794863803],
cmplx![0.2795081690440445, 18.419387948637894],
cmplx![-0.012331830955960252, 18.734267948637896],
cmplx![0.14223316576854003, 0.10023607930398608],
cmplx![0.04848316576854002, 0.006486079303986134],
cmplx![0.009893695996688708, 0.014372034600746361],
],
];
for kit in KS.iter().enumerate() {
Expand Down
1 change: 1 addition & 0 deletions crates/ekore/src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
pub mod anomalous_dimensions;
mod constants;
pub mod harmonics;
pub mod util;

/// References
pub mod bib {
Expand Down
6 changes: 6 additions & 0 deletions crates/ekore/src/util.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
#[macro_export]
macro_rules! cmplx {
($re:expr, $im:expr) => {
Complex::new($re, $im)
};
}

0 comments on commit 408b5be

Please sign in to comment.