From 00e513bd294e8a2fe1b8cb50784ff309c8db440f Mon Sep 17 00:00:00 2001 From: Mark Schultz-Wu Date: Sun, 27 Jul 2025 19:23:54 -0700 Subject: [PATCH 1/6] feat(bench): Add criterion benchmarking suite This introduces a new benchmark suite using the `criterion` crate to establish a performance baseline for the library's core functions. A `benches/performance.rs` file is added to benchmark a representative set of functions, including `cos`, `sin`, `atan`, `ln`, and `asin`, across various inputs. The `ln` function has been made public to allow it to be accessed by the benchmark runner. --- Cargo.toml | 7 ++++- benches/performance.rs | 64 ++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 +- 3 files changed, 71 insertions(+), 2 deletions(-) create mode 100644 benches/performance.rs diff --git a/Cargo.toml b/Cargo.toml index 89e5b5b..f175d27 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -12,4 +12,9 @@ homepage = "https://github.com/michaelciraci/trig-const" keywords = ["const", "trig", "no_std"] license = "MIT" -[dependencies] +[dev-dependencies] +criterion = "0.7.0" + +[[bench]] +name = "performance" +harness = false diff --git a/benches/performance.rs b/benches/performance.rs new file mode 100644 index 0000000..5979eda --- /dev/null +++ b/benches/performance.rs @@ -0,0 +1,64 @@ +use criterion::{criterion_group, criterion_main, Criterion}; +use std::hint::black_box; +use trig_const::{acos, asin, asinh, atan, cos, ln, sin, tan}; + +/// Benchmarks for the core trigonometric functions (sin, cos, tan). +fn bench_core_trig(c: &mut Criterion) { + let mut group = c.benchmark_group("Core Trig"); + + // Test a simple value. + group.bench_function("cos(1.5)", |b| b.iter(|| cos(black_box(1.5)))); + group.bench_function("sin(1.5)", |b| b.iter(|| sin(black_box(1.5)))); + + // Test a value that requires range reduction. + group.bench_function("cos(10.0)", |b| b.iter(|| cos(black_box(10.0)))); + + // Test tan, which involves two function calls and a division. + group.bench_function("tan(1.0)", |b| b.iter(|| tan(black_box(1.0)))); + + group.finish(); +} + +/// Benchmarks for the inverse trigonometric functions. +fn bench_inverse_trig(c: &mut Criterion) { + let mut group = c.benchmark_group("Inverse Trig"); + + // Test asin in its fast region (no range reduction needed). + group.bench_function("asin(0.4)", |b| b.iter(|| asin(black_box(0.4)))); + // Test asin in its slow region (requires range reduction). + group.bench_function("asin(0.9)", |b| b.iter(|| asin(black_box(0.9)))); + + // Test acos, which derives from asin. + group.bench_function("acos(0.5)", |b| b.iter(|| acos(black_box(0.5)))); + + // Test atan in its fast region. + group.bench_function("atan(0.5)", |b| b.iter(|| atan(black_box(0.5)))); + // Test atan at its slowest-converging input. + group.bench_function("atan(0.99)", |b| b.iter(|| atan(black_box(0.99)))); + + group.finish(); +} + +/// Benchmarks for logarithmic and hyperbolic functions. +fn bench_log_hyperbolic(c: &mut Criterion) { + let mut group = c.benchmark_group("Log & Hyperbolic"); + + // Test ln where its series converges relatively quickly. + group.bench_function("ln(1.1)", |b| b.iter(|| ln(black_box(1.1)))); + // Test ln at its slowest-converging input. + group.bench_function("ln(1.99)", |b| b.iter(|| ln(black_box(1.99)))); + + // Test asinh, which depends on ln and sqrt. + group.bench_function("asinh(2.0)", |b| b.iter(|| asinh(black_box(2.0)))); + + group.finish(); +} + +// Register all benchmark groups with criterion's main harness. +criterion_group!( + benches, + bench_core_trig, + bench_inverse_trig, + bench_log_hyperbolic +); +criterion_main!(benches); diff --git a/src/lib.rs b/src/lib.rs index 4d0d63f..68d7ac0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -452,7 +452,7 @@ const fn sqrt(x: f64) -> f64 { } /// Computes natural log using Taylor series approximation -const fn ln(x: f64) -> f64 { +pub const fn ln(x: f64) -> f64 { if x.is_nan() || x < 0.0 { return f64::NAN; } else if x == 0.0 { From 257fa46e1f43a28b2e6ed34aa9d45e2b18e9999c Mon Sep 17 00:00:00 2001 From: Mark Schultz-Wu Date: Sun, 27 Jul 2025 19:51:30 -0700 Subject: [PATCH 2/6] fix(atan): Improve precision and performance The previous `atan` implementation used a high-iteration Taylor series that was slow and only provided ~5-6 decimal places of precision. This refactors the function to use iterative argument reduction. The input is now mapped to a small value (`< 0.125`) before applying a much shorter, more efficient Taylor series. This new approach achieves full f64 machine precision while also providing a significant performance improvement. --- src/lib.rs | 25 ++++++++++++++----------- tests/brute.rs | 4 ++-- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 68d7ac0..08e147c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -56,9 +56,6 @@ use core::f64::{ const TAYLOR_SERIES_SUMS: usize = 16; /// Number of sum iterations for ln const LN_SUM_TERMS: f64 = 1001.0; -/// Number of sum iterations for atan. This series -/// takes a while to converge -const ATAN_SUMS: usize = 100_000; /// Cosine /// @@ -300,8 +297,7 @@ pub const fn atan(x: f64) -> f64 { let x_squared = x * x; let mut n = 0; - // This series takes a bit longer to converge - while n < ATAN_SUMS { + while n < TAYLOR_SERIES_SUMS { let denom = (2 * n + 1) as f64; s += sign * term / denom; term *= x_squared; @@ -315,13 +311,20 @@ pub const fn atan(x: f64) -> f64 { s } - - if x > 1.0 { - FRAC_PI_2 - atan_taylor_series(1.0 / x) - } else if x < -1.0 { - -FRAC_PI_2 - atan_taylor_series(1.0 / x) + let sign = x.signum(); + let mut val = x.abs(); + if val > 1.0 { + sign * (FRAC_PI_2 - atan(1.0 / val)) } else { - atan_taylor_series(x) + // Range reduction using the identity + // atan(x) = 2 atan(x / (1 + sqrt(1 + x^2))) + // until |x| \in [0, 0.125] + let mut multiplicand = 1.0; + while val > 0.125 { + val /= 1.0 + sqrt(1.0 + val * val); + multiplicand *= 2.0; + } + sign * multiplicand * atan_taylor_series(val) } } diff --git a/tests/brute.rs b/tests/brute.rs index 41bb59e..b62982b 100644 --- a/tests/brute.rs +++ b/tests/brute.rs @@ -58,7 +58,7 @@ fn test_acos() { #[test] fn test_atan() { for x in float_loop(-2.0 * PI, 2.0 * PI, 0.1) { - float_eq!(atan(x), x.atan(), 0.00001); + float_eq!(atan(x), x.atan()); } } @@ -66,7 +66,7 @@ fn test_atan() { fn test_atan2() { for x in float_loop(-2.0 * PI, 2.0 * PI, 0.1) { for y in float_loop(-2.0 * PI, 2.0 * PI, 0.1) { - float_eq!(atan2(x, y), x.atan2(y), 0.00001); + float_eq!(atan2(x, y), x.atan2(y)); } } } From 6672b8910075b843f1de40a1674cd5a42ae2a1a9 Mon Sep 17 00:00:00 2001 From: Mark Schultz-Wu Date: Sun, 27 Jul 2025 21:05:13 -0700 Subject: [PATCH 3/6] test(ln): add passing brute-force test for ln --- tests/brute.rs | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/tests/brute.rs b/tests/brute.rs index b62982b..14fada0 100644 --- a/tests/brute.rs +++ b/tests/brute.rs @@ -1,6 +1,6 @@ use core::f64::consts::PI; -use trig_const::{acos, acosh, asin, asinh, atan, atan2, cos, sin, tan}; +use trig_const::{acos, acosh, asin, asinh, atan, atan2, cos, ln, sin, tan}; fn float_loop(start: f64, stop: f64, step: f64) -> impl Iterator { core::iter::successors(Some(start), move |prev| { @@ -55,6 +55,14 @@ fn test_acos() { } } +#[test] +fn test_ln() { + for x in float_loop(0.01, 1000.0, 0.1) { + // Smallest n s.t. current impl passes with accuracy 10^{-n} + float_eq!(ln(x), x.ln(), 0.001); + } +} + #[test] fn test_atan() { for x in float_loop(-2.0 * PI, 2.0 * PI, 0.1) { From 50717ea5648f2c47bcedfed253aa768c67ee8223 Mon Sep 17 00:00:00 2001 From: Mark Schultz-Wu Date: Sun, 27 Jul 2025 21:20:38 -0700 Subject: [PATCH 4/6] feat(atanh): add atanh via taylor approximation Add a taylor approximation for atanh in preparation for a reimplementation of `ln`. --- benches/performance.rs | 18 ++++++++++++++++-- src/lib.rs | 36 ++++++++++++++++++++++++++++++++++++ tests/brute.rs | 9 ++++++++- 3 files changed, 60 insertions(+), 3 deletions(-) diff --git a/benches/performance.rs b/benches/performance.rs index 5979eda..55bd481 100644 --- a/benches/performance.rs +++ b/benches/performance.rs @@ -1,6 +1,6 @@ use criterion::{criterion_group, criterion_main, Criterion}; use std::hint::black_box; -use trig_const::{acos, asin, asinh, atan, cos, ln, sin, tan}; +use trig_const::{acos, asin, asinh, atan, atanh, cos, ln, sin, tan}; /// Benchmarks for the core trigonometric functions (sin, cos, tan). fn bench_core_trig(c: &mut Criterion) { @@ -54,11 +54,25 @@ fn bench_log_hyperbolic(c: &mut Criterion) { group.finish(); } +/// Benchmarks for the atanh function. +fn bench_atanh(c: &mut Criterion) { + let mut group = c.benchmark_group("atanh"); + + // Test a value where it converges quickly + group.bench_function("atanh(0.5)", |b| b.iter(|| atanh(black_box(0.5)))); + + // Test a value where it converges slowest (close to 1) + group.bench_function("atanh(0.99)", |b| b.iter(|| atanh(black_box(0.99)))); + + group.finish(); +} + // Register all benchmark groups with criterion's main harness. criterion_group!( benches, bench_core_trig, bench_inverse_trig, - bench_log_hyperbolic + bench_log_hyperbolic, + bench_atanh ); criterion_main!(benches); diff --git a/src/lib.rs b/src/lib.rs index 08e147c..4a89994 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -395,6 +395,42 @@ pub const fn acosh(x: f64) -> f64 { } } +/// Inverse hyperbolic tangent using range reduction. +pub const fn atanh(x: f64) -> f64 { + if x.is_nan() || x.abs() > 1.0 { + return f64::NAN; + } else if x == 1.0 { + return f64::INFINITY; + } else if x == -1.0 { + return f64::NEG_INFINITY; + } else if x == 0.0 { + return 0.0; + } + + let sign = x.signum(); + let mut y = x.abs(); + + // Iteratively reduce `|x|` until it's in `[0,0.125]` + // using the identity arctanh(x) = 2arctanh(x / (1 + sqrt(1 - x^2))) + let mut multiplicand = 1.0; + while y > 0.125 { + y /= 1.0 + sqrt(1.0 - y * y); + multiplicand *= 2.0; + } + + // atanh(y) = y + y^3/3 + y^5/5 + ... + let y_squared = y * y; + let mut term = y; + let mut s = y; + let mut n = 1; + while n < TAYLOR_SERIES_SUMS { + term *= y_squared; + s += term / ((2 * n + 1) as f64); + n += 1; + } + sign * multiplicand * s +} + /// e^x const fn exp(x: f64) -> f64 { let mut i = 1; diff --git a/tests/brute.rs b/tests/brute.rs index 14fada0..e70eac2 100644 --- a/tests/brute.rs +++ b/tests/brute.rs @@ -1,6 +1,6 @@ use core::f64::consts::PI; -use trig_const::{acos, acosh, asin, asinh, atan, atan2, cos, ln, sin, tan}; +use trig_const::{acos, acosh, asin, asinh, atan, atan2, atanh, cos, ln, sin, tan}; fn float_loop(start: f64, stop: f64, step: f64) -> impl Iterator { core::iter::successors(Some(start), move |prev| { @@ -92,3 +92,10 @@ fn test_acosh() { float_eq!(acosh(x), x.acosh()); } } + +#[test] +fn test_atanh() { + for x in float_loop(-0.99, 0.99, 0.01) { + float_eq!(atanh(x), x.atanh()); + } +} From 2e9a905b75818a6f85b9c7fde8596c58a87675cf Mon Sep 17 00:00:00 2001 From: Mark Schultz-Wu Date: Sun, 27 Jul 2025 21:22:54 -0700 Subject: [PATCH 5/6] feat(ln): higher precision, faster implementation of `ln(x)` --- src/lib.rs | 35 +++-------------------------------- tests/brute.rs | 3 +-- 2 files changed, 4 insertions(+), 34 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 4a89994..2bf33fe 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -54,8 +54,6 @@ use core::f64::{ /// Number of sum iterations for Taylor series const TAYLOR_SERIES_SUMS: usize = 16; -/// Number of sum iterations for ln -const LN_SUM_TERMS: f64 = 1001.0; /// Cosine /// @@ -490,7 +488,8 @@ const fn sqrt(x: f64) -> f64 { current_guess } -/// Computes natural log using Taylor series approximation +/// Computes natural log using the identity +/// `ln(x) = 2 arctanh((x - 1) / (x+1))` pub const fn ln(x: f64) -> f64 { if x.is_nan() || x < 0.0 { return f64::NAN; @@ -501,35 +500,7 @@ pub const fn ln(x: f64) -> f64 { } else if x.is_infinite() { return f64::INFINITY; } - - // Put into form ln(x) = ln(a * 2^k) = ln(a) + k * ln(2) - - let mut a = x; - let mut k = 0; - - // Normalize `a` to [1.0, 2.0) - while a >= 2.0 { - a /= 2.0; - k += 1; - } - while a < 1.0 { - a *= 2.0; - k -= 1; - } - - let x = a - 1.0; - - let mut s = 0.0; - let mut term = x; - let mut n = 1.0; - - while n < LN_SUM_TERMS { - s += term; - n += 1.0; - term = -term * x * (n - 1.0) / n; - } - - s + (k as f64) * f64::consts::LN_2 + 2.0 * atanh((x - 1.0) / (x + 1.0)) } #[cfg(test)] diff --git a/tests/brute.rs b/tests/brute.rs index e70eac2..bfd5ea2 100644 --- a/tests/brute.rs +++ b/tests/brute.rs @@ -58,8 +58,7 @@ fn test_acos() { #[test] fn test_ln() { for x in float_loop(0.01, 1000.0, 0.1) { - // Smallest n s.t. current impl passes with accuracy 10^{-n} - float_eq!(ln(x), x.ln(), 0.001); + float_eq!(ln(x), x.ln()); } } From 15a84ec99bd0b9053d164b242d2162df10ca4345 Mon Sep 17 00:00:00 2001 From: Mark Schultz-Wu Date: Sun, 27 Jul 2025 22:33:19 -0700 Subject: [PATCH 6/6] fix: improve precision and performance of all trig functions This fixes major precision and performance issues, allowing all functions to pass a full brute-force test suite at machine precision. - Rewrites `sin` and `cos` to use robust argument reduction. This fixes precision bugs for large inputs and resolves the persistent test failures in `cot` and `csc`. - Rewrites `exp` to use range reduction, fixing bugs that caused failures in `sinh` and `cosh`. - Adds brute-force tests for all remaining public functions and makes the `cot` and `csc` tests more robust to avoid false negatives near poles. --- src/lib.rs | 177 +++++++++++++++++++++++++++++++++++++++---------- tests/brute.rs | 52 ++++++++++++++- 2 files changed, 193 insertions(+), 36 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 2bf33fe..be0881f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -52,6 +52,19 @@ use core::f64::{ consts::{FRAC_PI_2, PI}, }; +// A const-compatible floor function. +const fn floor(x: f64) -> f64 { + let i = x as i64; + let f = i as f64; + if f == x || x >= 0.0 { + f + } else { + // For negative x, casting truncates towards zero, so we subtract 1 + // to get the correct floor. + f - 1.0 + } +} + /// Number of sum iterations for Taylor series const TAYLOR_SERIES_SUMS: usize = 16; @@ -64,36 +77,46 @@ const TAYLOR_SERIES_SUMS: usize = 16; /// const COS_PI: f64 = cos(PI); /// float_eq(COS_PI, -1.0); /// ``` -pub const fn cos(mut x: f64) -> f64 { - // If value is large, fold into smaller value - while x < -0.1 { - x += 2.0 * PI; - } - while x > 2.0 * PI + 0.1 { - x -= 2.0 * PI; - } - let div = (x / PI) as u32; - x -= div as f64 * PI; - let sign = if div % 2 != 0 { -1.0 } else { 1.0 }; - - let mut result = 1.0; - let mut inter = 1.0; - let num = x * x; - - let mut i = 1; - while i <= TAYLOR_SERIES_SUMS { - let comp = 2.0 * i as f64; - let den = comp * (comp - 1.0); - inter *= num / den; - if i % 2 == 0 { - result += inter; - } else { - result -= inter; +pub const fn cos(x: f64) -> f64 { + // Handle edge cases. cos(inf) is undefined. + if x.is_nan() || x.is_infinite() { + return f64::NAN; + } + + const fn taylor_cos(x: f64) -> f64 { + let x_squared = x * x; + let mut s = 1.0; + let mut term = 1.0; + let mut n = 1; + while n <= TAYLOR_SERIES_SUMS { + let comp = 2.0 * n as f64; + term *= -x_squared / (comp * (comp - 1.0)); + s += term; + n += 1; } - i += 1; + s } - sign * result + const FRAC_PI_2: f64 = core::f64::consts::FRAC_PI_2; + const PI: f64 = core::f64::consts::PI; + const TWO_PI: f64 = 2.0 * PI; + + let x = x - floor(x / TWO_PI) * TWO_PI; + + // Corrected Quadrant Reduction + if x < FRAC_PI_2 { + // Quadrant I: [0, PI/2) + taylor_cos(x) + } else if x < PI { + // Quadrant II: [PI/2, PI) + -taylor_cos(PI - x) + } else if x < 1.5 * PI { + // Quadrant III: [PI, 1.5*PI) + -taylor_cos(x - PI) + } else { + // Quadrant IV: [1.5*PI, 2*PI) + taylor_cos(TWO_PI - x) + } } /// Sine @@ -106,7 +129,57 @@ pub const fn cos(mut x: f64) -> f64 { /// float_eq(SIN_PI, 0.0); /// ``` pub const fn sin(x: f64) -> f64 { - cos(x - PI / 2.0) + if x.is_nan() || x.is_infinite() { + return f64::NAN; + } + + const fn floor(x: f64) -> f64 { + let i = x as i64; + let f = i as f64; + if f == x || x >= 0.0 { + f + } else { + f - 1.0 + } + } + + // Taylor series for sin(x) = x - x^3/3! + x^5/5! - ... + // This is only called with a small input in [0, PI/2]. + const fn taylor_sin(x: f64) -> f64 { + let x_squared = x * x; + let mut s = x; + let mut term = x; + let mut n = 1; + while n <= TAYLOR_SERIES_SUMS { + let comp = 2.0 * n as f64 + 1.0; + term *= -x_squared / (comp * (comp - 1.0)); + s += term; + n += 1; + } + s + } + + const PI: f64 = core::f64::consts::PI; + const TWO_PI: f64 = 2.0 * PI; + const FRAC_PI_2: f64 = core::f64::consts::FRAC_PI_2; + + // Perform a correct Euclidean-style remainder to map x to [0, 2*PI). + let x = x - floor(x / TWO_PI) * TWO_PI; + + // Reduce to [0, PI/2] and apply the correct sign based on quadrant. + if x < FRAC_PI_2 { + // Quadrant I: [0, PI/2) + taylor_sin(x) + } else if x < PI { + // Quadrant II: [PI/2, PI) -> sin(x) = sin(PI - x) + taylor_sin(PI - x) + } else if x < 1.5 * PI { + // Quadrant III: [PI, 1.5*PI) -> sin(x) = -sin(x - PI) + -taylor_sin(x - PI) + } else { + // Quadrant IV: [1.5*PI, 2*PI) -> sin(x) = -sin(2*PI - x) + -taylor_sin(TWO_PI - x) + } } /// Tangent @@ -429,17 +502,51 @@ pub const fn atanh(x: f64) -> f64 { sign * multiplicand * s } -/// e^x +/// e^x using range reduction for high precision. const fn exp(x: f64) -> f64 { - let mut i = 1; - let mut s = 1.0; + if x.is_nan() { + return f64::NAN; + } else if x == f64::INFINITY { + return f64::INFINITY; + } else if x == f64::NEG_INFINITY { + return 0.0; + } - while i < 16 { - s += expi(x, i) / factorial(i as f64); - i += 1; + const LN2: f64 = core::f64::consts::LN_2; + let k = (x / LN2 + 0.5) as i64; + + // Add checks to handle overflow before bit manipulation. + // The maximum exponent for a normal f64 corresponds to k=1023. + // k=1024 results in INFINITY. + if k >= 1024 { + return f64::INFINITY; + } + // The minimum exponent for a normal f64 corresponds to roughly k=-1022. + // Anything smaller will underflow to zero. -1075 is a safe threshold. + if k <= -1075 { + return 0.0; } - s + let r = x - k as f64 * LN2; + + // Taylor series for e^r + let mut s = 1.0; + let mut term = 1.0; + let mut n = 1; + while n < 20 { + term *= r / n as f64; + s += term; + n += 1; + } + + // Calculate 2^k by directly constructing the f64's bit representation. + // An f64's exponent is stored with a bias of 1023. We shift the + // biased exponent left by 52 bits to position it correctly, leaving + // the mantissa as zero, effectively creating the number 1.0 * 2^k. + let two_k_bits = ((1023 + k) as u64) << 52; + let scale = f64::from_bits(two_k_bits); + + scale * s } /// x^pow diff --git a/tests/brute.rs b/tests/brute.rs index bfd5ea2..5b91601 100644 --- a/tests/brute.rs +++ b/tests/brute.rs @@ -1,6 +1,8 @@ use core::f64::consts::PI; -use trig_const::{acos, acosh, asin, asinh, atan, atan2, atanh, cos, ln, sin, tan}; +use trig_const::{ + acos, acosh, asin, asinh, atan, atan2, atanh, cos, cosh, cot, csc, ln, sec, sin, sinh, tan, +}; fn float_loop(start: f64, stop: f64, step: f64) -> impl Iterator { core::iter::successors(Some(start), move |prev| { @@ -98,3 +100,51 @@ fn test_atanh() { float_eq!(atanh(x), x.atanh()); } } + +#[test] +fn test_sinh() { + // Test over a moderate range to avoid huge numbers. + for x in float_loop(-5.0, 5.0, 0.1) { + float_eq!(sinh(x), x.sinh()); + } +} + +#[test] +fn test_cosh() { + // Test over a moderate range to avoid huge numbers. + for x in float_loop(-5.0, 5.0, 0.1) { + float_eq!(cosh(x), x.cosh()); + } +} + +#[test] +fn test_sec() { + for x in float_loop(-2.0 * PI, 2.0 * PI, 0.1) { + float_eq!(sec(x), 1.0 / x.cos()); + } +} + +#[test] +fn test_cot() { + for x in float_loop(-2.0 * PI, 2.0 * PI, 0.1) { + let our_cot = cot(x); + // Test the stable identity: cot(x) * sin(x) = cos(x) + if our_cot.is_finite() { + let std_sin = x.sin(); + let std_cos = x.cos(); + float_eq!(our_cot * std_sin, std_cos); + } + } +} + +#[test] +fn test_csc() { + for x in float_loop(-2.0 * PI, 2.0 * PI, 0.1) { + let our_csc = csc(x); + // Test the stable identity: csc(x) * sin(x) = 1 + if our_csc.is_finite() { + let std_sin = x.sin(); + float_eq!(our_csc * std_sin, 1.0_f64); + } + } +}