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..55bd481
--- /dev/null
+++ b/benches/performance.rs
@@ -0,0 +1,78 @@
+use criterion::{criterion_group, criterion_main, Criterion};
+use std::hint::black_box;
+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) {
+    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();
+}
+
+/// 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_atanh
+);
+criterion_main!(benches);
diff --git a/src/lib.rs b/src/lib.rs
index 4d0d63f..be0881f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -52,13 +52,21 @@ 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;
-/// 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
 ///
@@ -69,36 +77,46 @@ const ATAN_SUMS: usize = 100_000;
 /// 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
@@ -111,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
@@ -300,8 +368,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 +382,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)
     }
 }
 
@@ -392,17 +466,87 @@ pub const fn acosh(x: f64) -> f64 {
     }
 }
 
-/// e^x
+/// 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 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
@@ -451,8 +595,9 @@ const fn sqrt(x: f64) -> f64 {
     current_guess
 }
 
-/// Computes natural log using Taylor series approximation
-const fn ln(x: f64) -> f64 {
+/// 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;
     } else if x == 0.0 {
@@ -462,35 +607,7 @@ 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 41bb59e..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, cos, 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<Item = f64> {
     core::iter::successors(Some(start), move |prev| {
@@ -55,10 +57,17 @@ fn test_acos() {
     }
 }
 
+#[test]
+fn test_ln() {
+    for x in float_loop(0.01, 1000.0, 0.1) {
+        float_eq!(ln(x), x.ln());
+    }
+}
+
 #[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 +75,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));
         }
     }
 }
@@ -84,3 +93,58 @@ 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());
+    }
+}
+
+#[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);
+        }
+    }
+}