diff --git a/Cargo.toml b/Cargo.toml
index 3e6817851..acfea457a 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,36 +1,8 @@
-[package]
-authors = ["Jorge Aparicio <jorge@japaric.io>"]
-categories = ["no-std"]
-description = "libm in pure Rust"
-documentation = "https://docs.rs/libm"
-keywords = ["libm", "math"]
-license = "MIT OR Apache-2.0"
-name = "libm"
-repository = "https://github.com/rust-lang-nursery/libm"
-version = "0.2.0"
-edition = "2018"
-
-[features]
-# only used to run our test suite
-default = ['stable']
-stable = []
-
-# Generate tests which are random inputs and the outputs are calculated with
-# musl libc.
-musl-reference-tests = ['rand']
-
-# Used checked array indexing instead of unchecked array indexing in this
-# library.
-checked = []
-
 [workspace]
 members = [
+  "crates/libm",
+  "crates/libm-analyze",
+  "crates/libm-test",
   "crates/compiler-builtins-smoke-test",
   "crates/libm-bench",
-]
-
-[dev-dependencies]
-no-panic = "0.1.8"
-
-[build-dependencies]
-rand = { version = "0.6.5", optional = true }
+]
\ No newline at end of file
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index c89346c73..257ece7d1 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -47,14 +47,14 @@ jobs:
       - template: ci/azure-install-rust.yml
       - script: rustup target add wasm32-unknown-unknown
         displayName: "Install rust wasm target"
-      - script: cargo build --target wasm32-unknown-unknown
+      - script: cargo build -p libm --target wasm32-unknown-unknown
         displayName: "Build for wasm"
-      - script: cargo build --target wasm32-unknown-unknown --no-default-features
+      - script: cargo build -p libm --target wasm32-unknown-unknown --no-default-features
         displayName: "Build for wasm (no default features)"
     variables:
       TOOLCHAIN: nightly
 
-  - job: rustfmt
+  - job: tools
     pool:
       vmImage: ubuntu-16.04
     steps:
@@ -63,6 +63,10 @@ jobs:
         displayName: "install rustfmt"
       - bash: cargo fmt --all -- --check
         displayName: "check formatting"
+      - bash: rustup component add clippy
+        displayName: "install clippy"
+      - bash: cargo clippy --all -- -D clippy::pedantic
+        displayName: "check clippy"
 
   - job: compiler_builtins_works
     pool:
diff --git a/build.rs b/build.rs
deleted file mode 100644
index 9af6dec93..000000000
--- a/build.rs
+++ /dev/null
@@ -1,444 +0,0 @@
-use std::env;
-
-fn main() {
-    println!("cargo:rerun-if-changed=build.rs");
-
-    #[cfg(feature = "musl-reference-tests")]
-    musl_reference_tests::generate();
-
-    if !cfg!(feature = "checked") {
-        let lvl = env::var("OPT_LEVEL").unwrap();
-        if lvl != "0" {
-            println!("cargo:rustc-cfg=assert_no_panic");
-        }
-    }
-}
-
-#[cfg(feature = "musl-reference-tests")]
-mod musl_reference_tests {
-    use rand::seq::SliceRandom;
-    use rand::Rng;
-    use std::fs;
-    use std::process::Command;
-
-    // Number of tests to generate for each function
-    const NTESTS: usize = 500;
-
-    // These files are all internal functions or otherwise miscellaneous, not
-    // defining a function we want to test.
-    const IGNORED_FILES: &[&str] = &["fenv.rs"];
-
-    struct Function {
-        name: String,
-        args: Vec<Ty>,
-        ret: Vec<Ty>,
-        tests: Vec<Test>,
-    }
-
-    enum Ty {
-        F32,
-        F64,
-        I32,
-        Bool,
-    }
-
-    struct Test {
-        inputs: Vec<i64>,
-        outputs: Vec<i64>,
-    }
-
-    pub fn generate() {
-        let files = fs::read_dir("src/math")
-            .unwrap()
-            .map(|f| f.unwrap().path())
-            .collect::<Vec<_>>();
-
-        let mut math = Vec::new();
-        for file in files {
-            if IGNORED_FILES.iter().any(|f| file.ends_with(f)) {
-                continue;
-            }
-
-            println!("generating musl reference tests in {:?}", file);
-
-            let contents = fs::read_to_string(file).unwrap();
-            let mut functions = contents.lines().filter(|f| f.starts_with("pub fn"));
-            while let Some(function_to_test) = functions.next() {
-                math.push(parse(function_to_test));
-            }
-        }
-
-        // Generate a bunch of random inputs for each function. This will
-        // attempt to generate a good set of uniform test cases for exercising
-        // all the various functionality.
-        generate_random_tests(&mut math, &mut rand::thread_rng());
-
-        // After we have all our inputs, use the x86_64-unknown-linux-musl
-        // target to generate the expected output.
-        generate_test_outputs(&mut math);
-        //panic!("Boo");
-        // ... and now that we have both inputs and expected outputs, do a bunch
-        // of codegen to create the unit tests which we'll actually execute.
-        generate_unit_tests(&math);
-    }
-
-    /// A "poor man's" parser for the signature of a function
-    fn parse(s: &str) -> Function {
-        let s = eat(s, "pub fn ");
-        let pos = s.find('(').unwrap();
-        let name = &s[..pos];
-        let s = &s[pos + 1..];
-        let end = s.find(')').unwrap();
-        let args = s[..end]
-            .split(',')
-            .map(|arg| {
-                let colon = arg.find(':').unwrap();
-                parse_ty(arg[colon + 1..].trim())
-            })
-            .collect::<Vec<_>>();
-        let tail = &s[end + 1..];
-        let tail = eat(tail, " -> ");
-        let ret = parse_retty(tail.replace("{", "").trim());
-
-        return Function {
-            name: name.to_string(),
-            args,
-            ret,
-            tests: Vec::new(),
-        };
-
-        fn parse_ty(s: &str) -> Ty {
-            match s {
-                "f32" => Ty::F32,
-                "f64" => Ty::F64,
-                "i32" => Ty::I32,
-                "bool" => Ty::Bool,
-                other => panic!("unknown type `{}`", other),
-            }
-        }
-
-        fn parse_retty(s: &str) -> Vec<Ty> {
-            match s {
-                "(f32, f32)" => vec![Ty::F32, Ty::F32],
-                "(f32, i32)" => vec![Ty::F32, Ty::I32],
-                "(f64, f64)" => vec![Ty::F64, Ty::F64],
-                "(f64, i32)" => vec![Ty::F64, Ty::I32],
-                other => vec![parse_ty(other)],
-            }
-        }
-
-        fn eat<'a>(s: &'a str, prefix: &str) -> &'a str {
-            if s.starts_with(prefix) {
-                &s[prefix.len()..]
-            } else {
-                panic!("{:?} didn't start with {:?}", s, prefix)
-            }
-        }
-    }
-
-    fn generate_random_tests<R: Rng>(functions: &mut [Function], rng: &mut R) {
-        for function in functions {
-            for _ in 0..NTESTS {
-                function.tests.push(generate_test(function, rng));
-            }
-        }
-
-        fn generate_test<R: Rng>(function: &Function, rng: &mut R) -> Test {
-            let mut inputs = function
-                .args
-                .iter()
-                .map(|ty| ty.gen_i64(rng))
-                .collect::<Vec<_>>();
-
-            // First argument to this function appears to be a number of
-            // iterations, so passing in massive random numbers causes it to
-            // take forever to execute, so make sure we're not running random
-            // math code until the heat death of the universe.
-            if function.name == "jn" || function.name == "jnf" {
-                inputs[0] &= 0xffff;
-            }
-
-            Test {
-                inputs,
-                // zero output for now since we'll generate it later
-                outputs: vec![],
-            }
-        }
-    }
-
-    impl Ty {
-        fn gen_i64<R: Rng>(&self, r: &mut R) -> i64 {
-            use std::f32;
-            use std::f64;
-
-            return match self {
-                Ty::F32 => {
-                    if r.gen_range(0, 20) < 1 {
-                        let i = *[f32::NAN, f32::INFINITY, f32::NEG_INFINITY]
-                            .choose(r)
-                            .unwrap();
-                        i.to_bits().into()
-                    } else {
-                        r.gen::<f32>().to_bits().into()
-                    }
-                }
-                Ty::F64 => {
-                    if r.gen_range(0, 20) < 1 {
-                        let i = *[f64::NAN, f64::INFINITY, f64::NEG_INFINITY]
-                            .choose(r)
-                            .unwrap();
-                        i.to_bits() as i64
-                    } else {
-                        r.gen::<f64>().to_bits() as i64
-                    }
-                }
-                Ty::I32 => {
-                    if r.gen_range(0, 10) < 1 {
-                        let i = *[i32::max_value(), 0, i32::min_value()].choose(r).unwrap();
-                        i.into()
-                    } else {
-                        r.gen::<i32>().into()
-                    }
-                }
-                Ty::Bool => r.gen::<bool>() as i64,
-            };
-        }
-
-        fn libc_ty(&self) -> &'static str {
-            match self {
-                Ty::F32 => "f32",
-                Ty::F64 => "f64",
-                Ty::I32 => "i32",
-                Ty::Bool => "i32",
-            }
-        }
-
-        fn libc_pty(&self) -> &'static str {
-            match self {
-                Ty::F32 => "*mut f32",
-                Ty::F64 => "*mut f64",
-                Ty::I32 => "*mut i32",
-                Ty::Bool => "*mut i32",
-            }
-        }
-
-        fn default(&self) -> &'static str {
-            match self {
-                Ty::F32 => "0_f32",
-                Ty::F64 => "0_f64",
-                Ty::I32 => "0_i32",
-                Ty::Bool => "false",
-            }
-        }
-
-        fn to_i64(&self) -> &'static str {
-            match self {
-                Ty::F32 => ".to_bits() as i64",
-                Ty::F64 => ".to_bits() as i64",
-                Ty::I32 => " as i64",
-                Ty::Bool => " as i64",
-            }
-        }
-    }
-
-    fn generate_test_outputs(functions: &mut [Function]) {
-        let mut src = String::new();
-        let dst = std::env::var("OUT_DIR").unwrap();
-
-        // Generate a program which will run all tests with all inputs in
-        // `functions`. This program will write all outputs to stdout (in a
-        // binary format).
-        src.push_str("use std::io::Write;");
-        src.push_str("fn main() {");
-        src.push_str("let mut result = Vec::new();");
-        for function in functions.iter_mut() {
-            src.push_str("unsafe {");
-            src.push_str("extern { fn ");
-            src.push_str(&function.name);
-            src.push_str("(");
-
-            let (ret, retptr) = match function.name.as_str() {
-                "sincos" | "sincosf" => (None, &function.ret[..]),
-                _ => (Some(&function.ret[0]), &function.ret[1..]),
-            };
-            for (i, arg) in function.args.iter().enumerate() {
-                src.push_str(&format!("arg{}: {},", i, arg.libc_ty()));
-            }
-            for (i, ret) in retptr.iter().enumerate() {
-                src.push_str(&format!("argret{}: {},", i, ret.libc_pty()));
-            }
-            src.push_str(")");
-            if let Some(ty) = ret {
-                src.push_str(" -> ");
-                src.push_str(ty.libc_ty());
-            }
-            src.push_str("; }");
-
-            src.push_str(&format!("static TESTS: &[[i64; {}]]", function.args.len()));
-            src.push_str(" = &[");
-            for test in function.tests.iter() {
-                src.push_str("[");
-                for val in test.inputs.iter() {
-                    src.push_str(&val.to_string());
-                    src.push_str(",");
-                }
-                src.push_str("],");
-            }
-            src.push_str("];");
-
-            src.push_str("for test in TESTS {");
-            for (i, arg) in retptr.iter().enumerate() {
-                src.push_str(&format!("let mut argret{} = {};", i, arg.default()));
-            }
-            src.push_str("let output = ");
-            src.push_str(&function.name);
-            src.push_str("(");
-            for (i, arg) in function.args.iter().enumerate() {
-                src.push_str(&match arg {
-                    Ty::F32 => format!("f32::from_bits(test[{}] as u32)", i),
-                    Ty::F64 => format!("f64::from_bits(test[{}] as u64)", i),
-                    Ty::I32 => format!("test[{}] as i32", i),
-                    Ty::Bool => format!("test[{}] as i32", i),
-                });
-                src.push_str(",");
-            }
-            for (i, _) in retptr.iter().enumerate() {
-                src.push_str(&format!("&mut argret{},", i));
-            }
-            src.push_str(");");
-            if let Some(ty) = &ret {
-                src.push_str(&format!("let output = output{};", ty.to_i64()));
-                src.push_str("result.extend_from_slice(&output.to_le_bytes());");
-            }
-
-            for (i, ret) in retptr.iter().enumerate() {
-                src.push_str(&format!(
-                    "result.extend_from_slice(&(argret{}{}).to_le_bytes());",
-                    i,
-                    ret.to_i64(),
-                ));
-            }
-            src.push_str("}");
-
-            src.push_str("}");
-        }
-
-        src.push_str("std::io::stdout().write_all(&result).unwrap();");
-
-        src.push_str("}");
-
-        let path = format!("{}/gen.rs", dst);
-        fs::write(&path, src).unwrap();
-
-        // Make it somewhat pretty if something goes wrong
-        drop(Command::new("rustfmt").arg(&path).status());
-
-        // Compile and execute this tests for the musl target, assuming we're an
-        // x86_64 host effectively.
-        let status = Command::new("rustc")
-            .current_dir(&dst)
-            .arg(&path)
-            .arg("--target=x86_64-unknown-linux-musl")
-            .status()
-            .unwrap();
-        assert!(status.success());
-        let output = Command::new("./gen").current_dir(&dst).output().unwrap();
-        assert!(output.status.success());
-        assert!(output.stderr.is_empty());
-
-        // Map all the output bytes back to an `i64` and then shove it all into
-        // the expected results.
-        let mut results = output.stdout.chunks_exact(8).map(|buf| {
-            let mut exact = [0; 8];
-            exact.copy_from_slice(buf);
-            i64::from_le_bytes(exact)
-        });
-
-        for f in functions.iter_mut() {
-            for test in f.tests.iter_mut() {
-                test.outputs = (0..f.ret.len()).map(|_| results.next().unwrap()).collect();
-            }
-        }
-        assert!(results.next().is_none());
-    }
-
-    /// Codegens a file which has a ton of `#[test]` annotations for all the
-    /// tests that we generated above.
-    fn generate_unit_tests(functions: &[Function]) {
-        let mut src = String::new();
-        let dst = std::env::var("OUT_DIR").unwrap();
-
-        for function in functions {
-            src.push_str("#[test]");
-            src.push_str("fn ");
-            src.push_str(&function.name);
-            src.push_str("_matches_musl() {");
-            src.push_str(&format!(
-                "static TESTS: &[([i64; {}], [i64; {}])]",
-                function.args.len(),
-                function.ret.len(),
-            ));
-            src.push_str(" = &[");
-            for test in function.tests.iter() {
-                src.push_str("([");
-                for val in test.inputs.iter() {
-                    src.push_str(&val.to_string());
-                    src.push_str(",");
-                }
-                src.push_str("],");
-                src.push_str("[");
-                for val in test.outputs.iter() {
-                    src.push_str(&val.to_string());
-                    src.push_str(",");
-                }
-                src.push_str("],");
-                src.push_str("),");
-            }
-            src.push_str("];");
-
-            src.push_str("for (test, expected) in TESTS {");
-            src.push_str("let output = ");
-            src.push_str(&function.name);
-            src.push_str("(");
-            for (i, arg) in function.args.iter().enumerate() {
-                src.push_str(&match arg {
-                    Ty::F32 => format!("f32::from_bits(test[{}] as u32)", i),
-                    Ty::F64 => format!("f64::from_bits(test[{}] as u64)", i),
-                    Ty::I32 => format!("test[{}] as i32", i),
-                    Ty::Bool => format!("test[{}] as i32", i),
-                });
-                src.push_str(",");
-            }
-            src.push_str(");");
-
-            for (i, ret) in function.ret.iter().enumerate() {
-                let get = if function.ret.len() == 1 {
-                    String::new()
-                } else {
-                    format!(".{}", i)
-                };
-                src.push_str(&(match ret {
-                    Ty::F32 => format!("if _eqf(output{}, f32::from_bits(expected[{}] as u32)).is_ok() {{ continue }}", get, i),
-                    Ty::F64 => format!("if _eq(output{}, f64::from_bits(expected[{}] as u64)).is_ok() {{ continue }}", get, i),
-                    Ty::I32 => format!("if output{} as i64 == expected[{}] {{ continue }}", get, i),
-                    Ty::Bool => unreachable!(),
-                }));
-            }
-
-            src.push_str(
-                r#"
-                panic!("INPUT: {:?} EXPECTED: {:?} ACTUAL {:?}", test, expected, output);
-            "#,
-            );
-            src.push_str("}");
-
-            src.push_str("}");
-        }
-
-        let path = format!("{}/musl-tests.rs", dst);
-        fs::write(&path, src).unwrap();
-
-        // Try to make it somewhat pretty
-        drop(Command::new("rustfmt").arg(&path).status());
-    }
-}
diff --git a/ci/run.sh b/ci/run.sh
index 37ffb8793..a9c2792c7 100755
--- a/ci/run.sh
+++ b/ci/run.sh
@@ -1,9 +1,19 @@
 #!/usr/bin/env sh
 
+if [ -z "$1" ]; then
+    echo "This script takes the $TARGET triple as its argument"
+    exit 1
+fi
+
 set -ex
+
 TARGET=$1
 
-CMD="cargo test --all --no-default-features --target $TARGET"
+CMD_="cargo test \
+  --manifest-path=crates/libm-test/Cargo.toml --all \
+  --no-default-features "
+
+CMD="${CMD_} --target ${TARGET}"
 
 $CMD
 $CMD --release
@@ -11,5 +21,11 @@ $CMD --release
 $CMD --features 'stable'
 $CMD --release --features 'stable'
 
-$CMD --features 'stable checked musl-reference-tests'
-$CMD --release --features  'stable checked musl-reference-tests'
+cargo build --release -p libm-test --features 'exhaustive'
+
+if [ "$TARGET" = "x86_64-unknown-linux-gnu" ]; then
+    export TARGET=x86_64-unknown-linux-musl
+
+    $CMD_ --target $TARGET --features 'stable checked system_libm'
+    $CMD_ --target $TARGET --release --features  'stable checked system_libm'
+fi
diff --git a/crates/compiler-builtins-smoke-test/src/lib.rs b/crates/compiler-builtins-smoke-test/src/lib.rs
index 7fad301b9..9407aa765 100644
--- a/crates/compiler-builtins-smoke-test/src/lib.rs
+++ b/crates/compiler-builtins-smoke-test/src/lib.rs
@@ -5,5 +5,5 @@
 #![allow(dead_code)]
 #![no_std]
 
-#[path = "../../../src/math/mod.rs"]
+#[path = "../../libm/src/math/mod.rs"]
 mod libm;
diff --git a/crates/libm-analyze/Cargo.toml b/crates/libm-analyze/Cargo.toml
new file mode 100644
index 000000000..7debe8948
--- /dev/null
+++ b/crates/libm-analyze/Cargo.toml
@@ -0,0 +1,19 @@
+[package]
+name = "libm-analyze"
+version = "0.1.0"
+authors = ["Gonzalo Brito Gadeschi <gonzalobg88@gmail.com>"]
+edition = "2018"
+
+[lib]
+proc-macro = true
+test = false
+
+[dependencies]
+syn = { version = "0.15", features = ["full", "extra-traits"] }
+quote = "0.6"
+walkdir = "2.2.8"
+proc-macro2 = "0.4"
+
+[features]
+default = ["analyze"]
+analyze = []
\ No newline at end of file
diff --git a/crates/libm-analyze/README.md b/crates/libm-analyze/README.md
new file mode 100644
index 000000000..742fbc852
--- /dev/null
+++ b/crates/libm-analyze/README.md
@@ -0,0 +1,22 @@
+# libm-analyze
+
+This crate provides a single macro, `for_each_api`:
+
+```rust
+macro_rules! nop {
+    (
+        id: $id:ident;
+        arg_tys: $($arg_tys:ty),*;
+        arg_ids: $($arg_ids:ident),*;
+        ret_ty: $ty:ty;
+    ) => {};
+}
+
+libm_analyze::for_each_api!(nop);
+```
+
+This macro takes a user-provided macro, and expands it for all libm APIs.
+
+For example, see how the `libm-test` crate `tests/system.rs` test uses it to
+test all `libm` APIs against random inputs, and verify the results against the
+system's libm library.
diff --git a/crates/libm-analyze/src/lib.rs b/crates/libm-analyze/src/lib.rs
new file mode 100644
index 000000000..be0aefd82
--- /dev/null
+++ b/crates/libm-analyze/src/lib.rs
@@ -0,0 +1,360 @@
+//! A proc macro to analyze the libm APIs We want to exhaustively match all
+// fields here to create a compilation error if new fields are added.
+#![allow(clippy::unneeded_field_pattern)]
+
+extern crate proc_macro;
+use self::proc_macro::TokenStream;
+use quote::quote;
+use std::collections::HashSet;
+use syn::parse_macro_input;
+
+/// `input` contains a single identifier, corresponding to a user-defined macro.
+/// This identifier is expanded for each libm public API.
+///
+/// See tests/analyze or below for the API.
+#[proc_macro]
+pub fn for_each_api(input: TokenStream) -> TokenStream {
+    let input = parse_macro_input!(input as Input);
+    let files = get_libm_files();
+    let functions = get_functions(&files, &input.ignored);
+    let mut tokens = proc_macro2::TokenStream::new();
+    let input_macro = input.macro_id;
+    for function in functions {
+        let id = function.ident;
+        let api_kind = function.api_kind;
+        let ret_ty = function.ret_ty;
+        let arg_tys = function.arg_tys;
+        let arg_ids = get_arg_ids(arg_tys.len());
+        let t = quote! {
+            #input_macro! {
+                id: #id;
+                api_kind: #api_kind;
+                arg_tys: #(#arg_tys),*;
+                arg_ids: #(#arg_ids),*;
+                ret_ty: #ret_ty;
+            }
+
+        };
+        tokens.extend(t);
+    }
+    tokens.into()
+}
+
+/// Traverses the libm crate directory, parsing all .rs files
+fn get_libm_files() -> Vec<syn::File> {
+    // Find the directory of the libm crate:
+    let root_dir = std::path::Path::new(env!("CARGO_MANIFEST_DIR"));
+    let libm_dir = root_dir
+        .parent()
+        .expect("couldn't access crates/ dir")
+        .join("libm");
+    let libm_src_dir = libm_dir.join("src");
+
+    // Traverse all Rust files, parsing them as `syn::File`
+    let mut files = Vec::new();
+    for entry in walkdir::WalkDir::new(libm_src_dir)
+        .into_iter()
+        .filter_map(Result::ok)
+    {
+        use std::io::Read;
+        let file_path = entry.path();
+        if file_path.is_dir()
+            || !file_path
+                .to_str()
+                .expect("can't format file path")
+                .ends_with(".rs")
+        {
+            // If the path is a directory or not a ".rs" file => skip it.
+            continue;
+        }
+
+        // Read the file into a string, and parse it into an AST using syn.
+        let mut file_string = String::new();
+        std::fs::File::open(&file_path)
+            .unwrap_or_else(|_| panic!("can't open file at path: {}", file_path.display()))
+            .read_to_string(&mut file_string)
+            .expect("failed to read file to string");
+        let file = syn::parse_file(&file_string).expect("failed to parse");
+        files.push(file);
+    }
+    files
+}
+
+/// Function signature that will be expanded for the user macro.
+struct FnSig {
+    ident: syn::Ident,
+    api_kind: syn::Ident,
+    c_abi: bool,
+    ret_ty: Option<syn::Type>,
+    arg_tys: Vec<syn::Type>,
+}
+
+impl FnSig {
+    fn name(&self) -> String {
+        self.ident.to_string()
+    }
+}
+
+macro_rules! syn_to_str {
+    ($e:expr) => {{
+        let t = $e;
+        let tokens = quote! {
+            #t
+        };
+        format!("{}", tokens)
+    }};
+}
+
+/// Extracts all public functions from the libm files while
+/// doing some sanity checks on the function signatures.
+fn get_functions(files: &[syn::File], ignored: &Option<HashSet<String>>) -> Vec<FnSig> {
+    let mut error = false;
+    let mut functions = Vec::new();
+    // Traverse all files matching function items
+    for item in files.iter().flat_map(|f| f.items.iter()) {
+        let mut e = false;
+        if let syn::Item::Fn(syn::ItemFn {
+            vis: syn::Visibility::Public(_),
+            ident,
+            constness,
+            asyncness,
+            unsafety,
+            attrs,
+            abi,
+            decl,
+            block: _,
+        }) = item
+        {
+            // Build a function signature while doing some sanity checks
+            let mut fn_sig = FnSig {
+                ident: ident.clone(),
+                api_kind: to_api_kind(ident.clone()),
+                c_abi: false,
+                arg_tys: Vec::new(),
+                ret_ty: None,
+            };
+            // Skip ignored functions:
+            if let Some(ignored) = ignored {
+                if ignored.contains(&fn_sig.name()) {
+                    continue;
+                }
+            }
+            macro_rules! err {
+                ($msg:expr) => {{
+                    #[cfg(feature = "analyze")]
+                    {
+                        eprintln!("[error]: Function \"{}\" {}", fn_sig.name(), $msg);
+                    }
+                    #[allow(unused_assignments)]
+                    {
+                        e = true;
+                    }
+                    ()
+                }};
+            }
+            if let Some(syn::Abi {
+                name: Some(l),
+                extern_token: _,
+            }) = abi
+            {
+                if l.value() == "C" {
+                    fn_sig.c_abi = true;
+                }
+            }
+            // If the function signature isn't extern "C", we aren't ABI compatible
+            // with libm.
+            if !fn_sig.c_abi {
+                err!("not `extern \"C\"`");
+            }
+            // Right now there are no const fn functions. We might add them
+            // in the future, and at that point, we should tune this here.
+            // In the mean time, error if somebody tries.
+            if constness.is_some() {
+                err!("is const");
+            }
+            // No functions should be async fn
+            if asyncness.is_some() {
+                err!("is async");
+            }
+            // FIXME: Math functions are not unsafe. Some functions in the
+            // libm C API take pointers, but in our API take repr(Rust)
+            // tuples (for some reason). Once we fix those to have the same
+            // API as C libm, we should use references on their signature
+            // instead, and make them safe.
+            if unsafety.is_some() {
+                let e2 = e;
+                err!("is unsafe");
+                e = e2;
+            }
+            let syn::FnDecl {
+                fn_token: _,
+                generics,
+                paren_token: _,
+                inputs,
+                variadic,
+                output,
+            } = (**decl).clone();
+
+            // Forbid generic parameters, lifetimes, and consts in public APIs:
+            if variadic.is_some() {
+                err!(format!(
+                    "contains variadic arguments \"{}\"",
+                    syn_to_str!(variadic.unwrap())
+                ));
+            }
+            if generics.type_params().count() != 0 {
+                err!(format!(
+                    "contains generic parameters \"{}\"",
+                    syn_to_str!(generics.clone())
+                ));
+            }
+            if generics.lifetimes().count() != 0 {
+                err!(format!(
+                    "contains lifetime parameters \"{}\"",
+                    syn_to_str!(generics.clone())
+                ));
+            }
+            if generics.const_params().count() != 0 {
+                err!(format!(
+                    "contains const parameters \"{}\"",
+                    syn_to_str!(generics.clone())
+                ));
+            }
+            if attrs.is_empty() {
+                err!("missing `#[inline]` and `#[no_panic]` attributes");
+            } else {
+                let attrs = attrs
+                    .iter()
+                    .map(|a| syn_to_str!(a))
+                    .collect::<Vec<_>>()
+                    .join(",");
+                if !attrs.contains("inline") {
+                    err!("missing `#[inline]` attribute");
+                }
+                if !attrs.contains("no_panic") {
+                    err!("missing `#[no_panic]` attributes");
+                }
+            }
+            // Validate and parse output parameters and function arguments:
+            match output {
+                syn::ReturnType::Default => (),
+                syn::ReturnType::Type(_, ref b) if valid_ty(&b) => fn_sig.ret_ty = Some(*b.clone()),
+                other => err!(format!("returns unsupported type {}", syn_to_str!(other))),
+            }
+            for input in inputs {
+                match input {
+                    syn::FnArg::Captured(ref c) if valid_ty(&c.ty) => {
+                        fn_sig.arg_tys.push(c.ty.clone())
+                    }
+                    other => err!(format!(
+                        "takes unsupported argument type {}",
+                        syn_to_str!(other)
+                    )),
+                }
+            }
+            // If there was an error, we skip the function.
+            // Otherwise, the user macro is expanded with
+            // the function:
+            if e {
+                error = true;
+            } else {
+                functions.push(fn_sig);
+            }
+        }
+    }
+    if error {
+        // too many errors:
+        //        panic!("errors found");
+    }
+    functions
+}
+
+/// Parses a type into a String - arg is true if the type is an argument, and
+/// false if its a return value.
+fn valid_ty(t: &syn::Type) -> bool {
+    match t {
+        syn::Type::Ptr(p) => {
+            let c = p.const_token.is_some();
+            let m = p.mutability.is_some();
+            assert!(!(c && m));
+            match &*p.elem {
+                syn::Type::Path(_) => valid_ty(&p.elem),
+                // Only one layer of pointers allowed:
+                _ => false,
+            }
+        }
+        syn::Type::Path(p) => {
+            assert!(p.qself.is_none());
+            assert_eq!(p.path.segments.len(), 1);
+            let s = p
+                .path
+                .segments
+                .first()
+                .unwrap()
+                .into_value()
+                .ident
+                .to_string();
+            match s.as_str() {
+                "i8" | "i16" | "i32" | "i64" | "isize" | "u8" | "u16" | "u32" | "u64" | "usize"
+                | "f32" | "f64" => true,
+                _ => false,
+            }
+        }
+        _ => false,
+    }
+}
+
+/// Returns a vector containing `len` identifiers.
+fn get_arg_ids(len: usize) -> Vec<syn::Ident> {
+    let mut ids = Vec::new();
+    for i in 0..len {
+        let x = format!("x{}", i);
+        ids.push(syn::Ident::new(&x, proc_macro2::Span::call_site()));
+    }
+    ids
+}
+
+/// Returns the `ApiKind` enum variant for this function
+fn to_api_kind(id: syn::Ident) -> syn::Ident {
+    let name = syn_to_str!(id);
+    let first = name.chars().nth(0).unwrap();
+    let first_upper = first.to_uppercase().nth(0).unwrap();
+    let name = name.replacen(first, &first_upper.to_string(), 1);
+    syn::Ident::new(&name, proc_macro2::Span::call_site())
+}
+
+#[derive(Debug)]
+struct Input {
+    macro_id: syn::Ident,
+    ignored: Option<HashSet<String>>,
+}
+
+impl syn::parse::Parse for Input {
+    fn parse(input: syn::parse::ParseStream) -> syn::Result<Self> {
+        let content;
+        let macro_id: syn::Ident = input.parse()?;
+        let lookahead = input.lookahead1();
+        if lookahead.peek(syn::token::Paren) {
+            let _paren_token = syn::parenthesized!(content in input);
+            let ignored: syn::Lit = content.parse::<syn::Lit>()?;
+            if let syn::Lit::Str(c) = ignored {
+                let s = c.value();
+                let mut hash_set = HashSet::<String>::new();
+                for i in s.split(',') {
+                    hash_set.insert(i.to_string());
+                }
+                Ok(Self {
+                    macro_id,
+                    ignored: Some(hash_set),
+                })
+            } else {
+                Err(lookahead.error())
+            }
+        } else {
+            Ok(Self {
+                macro_id,
+                ignored: None,
+            })
+        }
+    }
+}
diff --git a/crates/libm-analyze/tests/analyze.rs b/crates/libm-analyze/tests/analyze.rs
new file mode 100644
index 000000000..0b1ba83b3
--- /dev/null
+++ b/crates/libm-analyze/tests/analyze.rs
@@ -0,0 +1,14 @@
+//! Tests that the proc-macro accepts macros with
+//! the following pattern:
+
+macro_rules! nop {
+    (
+        id: $id:ident;
+        api_kind: $api_kind:ident;
+        arg_tys: $($arg_tys:ty),*;
+        arg_ids: $($arg_ids:ident),*;
+        ret_ty: $ret_ty:ty;
+    ) => {};
+}
+
+libm_analyze::for_each_api!(nop);
diff --git a/crates/libm-analyze/tests/analyze_opt.rs b/crates/libm-analyze/tests/analyze_opt.rs
new file mode 100644
index 000000000..c5a2598eb
--- /dev/null
+++ b/crates/libm-analyze/tests/analyze_opt.rs
@@ -0,0 +1,14 @@
+//! Tests that the proc-macro accepts macros with
+//! the following pattern:
+
+macro_rules! nop {
+    (
+        id: $id:ident;
+        api_kind: $api_kind:ident;
+        arg_tys: $($arg_tys:ty),*;
+        arg_ids: $($arg_ids:ident),*;
+        ret_ty: $ret_ty:ty;
+    ) => {};
+}
+
+libm_analyze::for_each_api!(nop("j1f,jn"));
diff --git a/crates/libm-bench/Cargo.toml b/crates/libm-bench/Cargo.toml
index ba65dbd5f..24211d1bc 100644
--- a/crates/libm-bench/Cargo.toml
+++ b/crates/libm-bench/Cargo.toml
@@ -6,9 +6,10 @@ edition = "2018"
 license = "MIT OR Apache-2.0"
 
 [dependencies]
-libm = { path = "../..", default-features = false }
+libm = { path = "../libm", default-features = false }
+libm-analyze = { path = "../libm-analyze", default-features = false }
+libm-test = { path = "../libm-test", default-features = false }
 rand = "0.6.5"
-paste = "0.1.5"
 
 [features]
 default = []
diff --git a/crates/libm-bench/benches/bench.rs b/crates/libm-bench/benches/bench.rs
index b6d874153..17c2bb33a 100644
--- a/crates/libm-bench/benches/bench.rs
+++ b/crates/libm-bench/benches/bench.rs
@@ -1,115 +1,36 @@
 #![feature(test)]
 extern crate test;
 
+use libm_test::{ApiKind, CallFn};
 use rand::Rng;
 use test::Bencher;
 
-macro_rules! unary {
-  ($($func:ident),*) => ($(
-      paste::item! {
-        #[bench]
-        pub fn [<$func>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f64>();
-            bh.iter(|| test::black_box(libm::[<$func>](x)))
-        }
-        #[bench]
-        pub fn [<$func f>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f32>();
-            bh.iter(|| test::black_box(libm::[<$func f>](x)))
-        }
-    }
-  )*);
-}
-macro_rules! binary {
-  ($($func:ident),*) => ($(
-      paste::item! {
-        #[bench]
-        pub fn [<$func>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f64>();
-            let y = rng.gen::<f64>();
-            bh.iter(|| test::black_box(libm::[<$func>](x, y)))
-        }
-        #[bench]
-        pub fn [<$func f>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f32>();
-            let y = rng.gen::<f32>();
-            bh.iter(|| test::black_box(libm::[<$func f>](x, y)))
-        }
-    }
-  )*);
-  ($($func:ident);*) => ($(
-      paste::item! {
-        #[bench]
-        pub fn [<$func>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f64>();
-            let n = rng.gen::<i32>();
-            bh.iter(|| test::black_box(libm::[<$func>](x, n)))
-        }
-        #[bench]
-        pub fn [<$func f>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f32>();
-            let n = rng.gen::<i32>();
-            bh.iter(|| test::black_box(libm::[<$func f>](x, n)))
-        }
-    }
-  )*);
-}
-macro_rules! trinary {
-  ($($func:ident),*) => ($(
-      paste::item! {
-        #[bench]
-        pub fn [<$func>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f64>();
-            let y = rng.gen::<f64>();
-            let z = rng.gen::<f64>();
-            bh.iter(|| test::black_box(libm::[<$func>](x, y, z)))
-        }
-        #[bench]
-        pub fn [<$func f>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let x = rng.gen::<f32>();
-            let y = rng.gen::<f32>();
-            let z = rng.gen::<f32>();
-            bh.iter(|| test::black_box(libm::[<$func f>](x, y, z)))
-        }
-    }
-  )*);
-}
-macro_rules! bessel {
-  ($($func:ident),*) => ($(
-      paste::item! {
-        #[bench]
-        pub fn [<$func>](bh: &mut Bencher) {
-            let mut rng = rand::thread_rng();
-            let mut n = rng.gen::<i32>();
-            n &= 0xffff;
-            let x = rng.gen::<f64>();
-            bh.iter(|| test::black_box(libm::[<$func>](n, x)))
-        }
-        #[bench]
-        pub fn [<$func f>](bh: &mut Bencher) {
+macro_rules! bench_fn {
+    (
+        id: $id:ident;
+        api_kind: $api_kind:ident;
+        arg_tys: $($arg_tys:ty),*;
+        arg_ids: $($arg_ids:ident),*;
+        ret_ty: $ret_ty:ty;
+    ) => {
+        #[bench]
+        #[allow(unused_mut)]
+        pub fn $id(bh: &mut Bencher) {
+            type FnTy
+                = unsafe extern "C" fn ($($arg_ids: $arg_tys),*) -> $ret_ty;
+
+            // Generate a tuple of arguments containing random values:
             let mut rng = rand::thread_rng();
-            let mut n = rng.gen::<i32>();
-            n &= 0xffff;
-            let x = rng.gen::<f32>();
-            bh.iter(|| test::black_box(libm::[<$func f>](n, x)))
+            let mut x: ( $($arg_tys,)+ ) = ( $(rng.gen::<$arg_tys>(),)+ );
+
+            if let ApiKind::Jn = ApiKind::$api_kind {
+                let ptr = &mut x as *mut _ as *mut i32;
+                unsafe { ptr.write(ptr.read() & 0xffff) };
+            }
+
+            bh.iter(|| test::black_box(x).call(libm::$id as FnTy))
         }
-    }
-  )*);
+    };
 }
 
-unary!(
-    acos, acosh, asin, atan, cbrt, ceil, cos, cosh, erf, exp, exp2, exp10, expm1, fabs, floor, j0,
-    j1, lgamma, log, log1p, log2, log10, round, sin, sinh, sqrt, tan, tanh, tgamma, trunc, y0, y1
-);
-binary!(atan2, copysign, fdim, fmax, fmin, fmod, hypot, pow);
-trinary!(fma);
-bessel!(jn, yn);
-binary!(ldexp; scalbn);
+libm_analyze::for_each_api!(bench_fn(/*ignore:*/ "jnf"));
diff --git a/crates/libm-test/Cargo.toml b/crates/libm-test/Cargo.toml
new file mode 100644
index 000000000..46691669a
--- /dev/null
+++ b/crates/libm-test/Cargo.toml
@@ -0,0 +1,19 @@
+[package]
+name = "libm-test"
+version = "0.1.0"
+authors = ["Gonzalo Brito Gadeschi <gonzalobg88@gmail.com>"]
+edition = "2018"
+
+[dependencies]
+rand = "0.7"
+libm = { path = "../libm", default-features = false }
+
+[dev-dependencies]
+libm-analyze = { path = "../libm-analyze", default-features = true }
+
+[features]
+default = []
+checked = ["libm/checked"]
+stable = ["libm/stable"]
+system_libm = []
+exhaustive = []
diff --git a/crates/libm-test/README.md b/crates/libm-test/README.md
new file mode 100644
index 000000000..3e420e9fb
--- /dev/null
+++ b/crates/libm-test/README.md
@@ -0,0 +1,9 @@
+# libm-test
+
+This crate contains different types of test for the `libm` crate.
+
+* `tests/system.rs`: generate random inputs, and tests that the results of the
+  `libm` crate are within the tolerance required by the IEEE from those of the
+  system's libm library (e.g. musl, glibc's libm, libSystem_m, etc.).
+  
+* `tests/unit.rs`: contains some small unit tests. 
diff --git a/crates/libm-test/build.rs b/crates/libm-test/build.rs
new file mode 100644
index 000000000..345e60378
--- /dev/null
+++ b/crates/libm-test/build.rs
@@ -0,0 +1,22 @@
+use std::env;
+
+fn main() {
+    println!("cargo:rerun-if-changed=build.rs");
+    let target = env::var("TARGET").expect("TARGET was not set");
+    let profile = env::var("PROFILE").unwrap_or(String::new());
+    let opt_level: i32 = env::var("OPT_LEVEL").unwrap().parse().unwrap();
+    if !cfg!(feature = "checked") {
+        if opt_level != 0 {
+            println!("cargo:rustc-cfg=assert_no_panic");
+        }
+    }
+
+    if profile == "release" || opt_level > 0 {
+        match target.as_str() {
+            "x86_64-unknown-linux-gnu" | "x86_64-apple-darwin" | "x86_64-pc-windows-msvc" => {
+                println!("cargo:rustc-cfg=exhaustive32");
+            }
+            _ => (),
+        }
+    }
+}
diff --git a/crates/libm-test/src/api_kind.rs b/crates/libm-test/src/api_kind.rs
new file mode 100644
index 000000000..d1ca5a1c6
--- /dev/null
+++ b/crates/libm-test/src/api_kind.rs
@@ -0,0 +1,120 @@
+//! An enum with one variant per API
+
+// FIXME: this should be auto-generated
+
+/// Kind of libm API - used to handle generating tests
+/// for some functions slightly differently.
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub enum ApiKind {
+    Acos,
+    Acosf,
+    Acosh,
+    Acoshf,
+    Asin,
+    Asinf,
+    Asinh,
+    Asinhf,
+    Atan,
+    Atan2,
+    Atan2f,
+    Atanf,
+    Atanh,
+    Atanhf,
+    Cbrt,
+    Cbrtf,
+    Ceil,
+    Ceilf,
+    Copysign,
+    Copysignf,
+    Cos,
+    Cosf,
+    Cosh,
+    Coshf,
+    Erf,
+    Erfc,
+    Erfcf,
+    Erff,
+    Exp,
+    Exp10,
+    Exp10f,
+    Exp2,
+    Exp2f,
+    Expf,
+    Expm1,
+    Expm1f,
+    Fabs,
+    Fabsf,
+    Fdim,
+    Fdimf,
+    Floor,
+    Floorf,
+    Fma,
+    Fmaf,
+    Fmax,
+    Fmaxf,
+    Fmin,
+    Fminf,
+    Fmod,
+    Fmodf,
+    Hypot,
+    Hypotf,
+    Ilogb,
+    Ilogbf,
+    J0,
+    J0f,
+    J1,
+    J1f,
+    Jn,
+    Ldexp,
+    Ldexpf,
+    Lgamma,
+    Lgammaf,
+    Log,
+    Log10,
+    Log10f,
+    Log1p,
+    Log1pf,
+    Log2,
+    Log2f,
+    Logf,
+    Modf,
+    Modff,
+    Nextafter,
+    Nextafterf,
+    Nexttoward,
+    Nexttowardf,
+    Pow,
+    Powf,
+    Remainder,
+    Remainderf,
+    Rint,
+    Rintf,
+    Round,
+    Roundf,
+    Scalbn,
+    Scalbnf,
+    Sin,
+    SinCos,
+    SinCosPi,
+    SinCosPif,
+    SinCosf,
+    Sinf,
+    Sinh,
+    Sinhf,
+    Sqrt,
+    Sqrtf,
+    Tan,
+    Tanf,
+    Tanh,
+    Tanhf,
+    Tgamma,
+    Tgammaf,
+    Trunc,
+    Truncf,
+    Y0,
+    Y0f,
+    Y1,
+    Y1f,
+    Yn,
+    Ynf,
+}
diff --git a/crates/libm-test/src/call_fn.rs b/crates/libm-test/src/call_fn.rs
new file mode 100644
index 000000000..0d879ac30
--- /dev/null
+++ b/crates/libm-test/src/call_fn.rs
@@ -0,0 +1,38 @@
+//! Trait to do function dispatch based on tuples:
+
+/// This implements function dispatch for tuples of arguments used in the tests
+/// above, so that we can: (f32, 32).call(fn(f32, f32) -> f32) generically.
+///
+/// We need the input parameter F to support dispatching, e.g., (f32,f32) with
+/// functions that return both f32 or i32. Those are two different types, so we
+/// need to be parametric over them.
+pub trait CallFn<F> {
+    type Ret;
+    fn call(self, f: F) -> Self::Ret;
+}
+
+macro_rules! impl_call {
+    (($($arg_tys:ty),*) -> $ret_ty:ty: $self_:ident: $($xs:expr),*)  => {
+        impl CallFn<unsafe extern"C" fn($($arg_tys),*) -> $ret_ty> for ($($arg_tys,)+) {
+            type Ret = $ret_ty;
+            fn call(self, f: unsafe extern "C" fn($($arg_tys),*) -> $ret_ty) -> Self::Ret {
+                let $self_ = self;
+                unsafe { f($($xs),*) }
+            }
+        }
+    };
+}
+
+impl_call!((f32) -> f32: x: x.0);
+impl_call!((f64) -> f64: x: x.0);
+impl_call!((f64) -> i32: x: x.0);
+impl_call!((f32) -> i32: x: x.0);
+impl_call!((f32, f32) -> f32: x: x.0, x.1);
+impl_call!((f32, f64) -> f32: x: x.0, x.1);
+impl_call!((f64, f64) -> f64: x: x.0, x.1);
+impl_call!((f64, i32) -> f64: x: x.0, x.1);
+impl_call!((f32, i32) -> f32: x: x.0, x.1);
+impl_call!((i32, f64) -> f64: x: x.0, x.1);
+impl_call!((i32, f32) -> f32: x: x.0, x.1);
+impl_call!((f32, f32, f32) -> f32: x: x.0, x.1, x.2);
+impl_call!((f64, f64, f64) -> f64: x: x.0, x.1, x.2);
diff --git a/crates/libm-test/src/float_rng.rs b/crates/libm-test/src/float_rng.rs
new file mode 100644
index 000000000..07ea8ecdf
--- /dev/null
+++ b/crates/libm-test/src/float_rng.rs
@@ -0,0 +1,49 @@
+//! Traits for generating ranges of floating-point values
+
+pub trait Toward: Sized {
+    /// Generates `len` minimal incremental steps from self to other.
+    ///
+    /// Other is often never reached.
+    fn toward(self, other: Self, len: usize) -> Vec<Self>;
+}
+
+pub trait Distribute: Sized {
+    /// Distributes `len` values in range `[self, other]`.
+    fn distribute(self, other: Self, len: usize) -> Vec<Self>;
+}
+
+macro_rules! impl_f {
+    ($float_ty:ident, $toward_fn:path) => {
+        impl Toward for $float_ty {
+            fn toward(self, other: Self, len: usize) -> Vec<Self> {
+                let mut vec = Vec::with_capacity(len);
+                let mut current = self;
+                vec.push(self);
+                for _ in 0..=len {
+                    current = $toward_fn(current, other as _);
+                    vec.push(self);
+                    if current.to_bits() == other.to_bits() {
+                        break;
+                    }
+                }
+                vec
+            }
+        }
+        impl Distribute for $float_ty {
+            fn distribute(self, other: Self, mut len: usize) -> Vec<Self> {
+                let mut vec = Vec::with_capacity(len + 1);
+                vec.push(self);
+                // Bresenham's alg:
+                let mut x = self;
+                while len > 0 {
+                    x += (other - self) / (len as Self);
+                    len -= 1;
+                    vec.push(x);
+                }
+                vec
+            }
+        }
+    };
+}
+impl_f!(f32, libm::nextafterf);
+impl_f!(f64, libm::nextafter);
diff --git a/crates/libm-test/src/lib.rs b/crates/libm-test/src/lib.rs
new file mode 100644
index 000000000..692579bb2
--- /dev/null
+++ b/crates/libm-test/src/lib.rs
@@ -0,0 +1,42 @@
+//! Testing utilities required by most tests.
+
+mod within_ulps;
+pub use self::within_ulps::WithinUlps;
+
+mod call_fn;
+pub use self::call_fn::CallFn;
+
+mod tuple_vec;
+pub use self::tuple_vec::TupleVec;
+
+mod api_kind;
+pub use self::api_kind::ApiKind;
+
+mod float_rng;
+pub use self::float_rng::{Distribute, Toward};
+
+mod rand_seq;
+pub use self::rand_seq::RandSeq;
+
+/// Asserts that two values are approximately equal up-to ULP tolerance
+#[macro_export]
+macro_rules! assert_approx_eq {
+    ($result:ident == $expected:ident,
+     id: $id:ident, arg: $arg:ident, ulp: $ulps:expr) => {
+        if !$crate::WithinUlps::within_ulps($result, $expected, $ulps) {
+            let f = format!(
+                "{}{:?} returns = {:?} != {:?} (expected)",
+                stringify!($id),
+                $arg,
+                $result,
+                $expected
+            );
+            panic!(f);
+        }
+    };
+    ($result:expr, $expected:expr, ulp: $ulps:expr) => {
+        if !$crate::WithinUlps::within_ulps($result, $expected, $ulps) {
+            panic!("{:?} != {:?}", $result, $expected);
+        }
+    };
+}
diff --git a/crates/libm-test/src/rand_seq.rs b/crates/libm-test/src/rand_seq.rs
new file mode 100644
index 000000000..ab89b298f
--- /dev/null
+++ b/crates/libm-test/src/rand_seq.rs
@@ -0,0 +1,191 @@
+//! Trait for generating rand sequence testing relevant values for each API
+
+use crate::ApiKind;
+
+pub trait RandSeq: Sized {
+    /// Generates a sequence containing exactly `len` values computed from the RNG `rng`
+    /// according to the `api_kind`
+    fn rand_seq<R: rand::Rng>(rng: &mut R, api_kind: ApiKind, len: usize) -> Vec<Self>;
+}
+
+macro_rules! impl_rand_seq_f {
+    ($float_ty:ident) => {
+        #[allow(clippy::use_self)]
+        impl RandSeq for $float_ty {
+            fn rand_seq<R: rand::Rng>(rng: &mut R, api_kind: ApiKind, len: usize) -> Vec<Self> {
+                use crate::{Distribute, Toward};
+                use rand::seq::SliceRandom;
+                use std::$float_ty::*;
+
+                let mut vec = Vec::with_capacity(len);
+                let double = std::mem::size_of::<Self>() == 8;
+
+                // These inputs are always tested
+                const BOUNDS: [$float_ty; 9] = [
+                    NAN,
+                    INFINITY,
+                    NEG_INFINITY,
+                    EPSILON,
+                    -EPSILON,
+                    MAX,
+                    MIN,
+                    MIN_POSITIVE,
+                    -MIN_POSITIVE,
+                ];
+                vec.extend(&BOUNDS);
+
+                // The domain close to these inputs is also tested:
+                const LIMIT_NSTEPS: usize = 10_001;
+                vec.extend(INFINITY.toward(0., LIMIT_NSTEPS));
+                vec.extend(NEG_INFINITY.toward(0., LIMIT_NSTEPS));
+                vec.extend((0. as Self).toward(MIN_POSITIVE, LIMIT_NSTEPS));
+                vec.extend((0. as Self).toward(-MIN_POSITIVE, LIMIT_NSTEPS));
+                vec.extend((-1. as Self).distribute(1. as Self, 200_001));
+
+                // These domains are extended with tests specific to each
+                // type of API:
+                {
+                    macro_rules! dist {
+                        ($from:expr, $to:expr, $steps:expr) => {
+                            vec.extend(($from as Self).distribute($to as Self, $steps));
+                        };
+                    }
+                    use crate::ApiKind::*;
+                    match api_kind {
+                        Sin | Cos | Tan => {
+                            if double {
+                                dist!(1e-300, 1e+8, 200_001);
+                            } else {
+                                dist!(1e-30, 1e+8, 200_001);
+                            }
+                        }
+                        SinCos | SinCosPi => {
+                            dist!(-1e-14, 1e+14, 200_001);
+                        }
+                        Log | Log2 | Log10 | Log1p => {
+                            if double {
+                                dist!(-1e-300, 1e+14, 200_001);
+                            } else {
+                                dist!(-1e-30, 1e+14, 200_001);
+                            }
+                        }
+
+                        Exp | Exp2 | Exp10 | Expm1 => {
+                            dist!(-1000., 1000., 200_001);
+                        }
+                        Pow => {
+                            dist!(-100., 100., 200_001);
+                        }
+                        Cbrt => {
+                            dist!(1e-14, 1e+14, 100_001);
+                            dist!(-1e-14, -1e+14, 100_001);
+                        }
+                        Hypot => {
+                            dist!(-1e-7, 1e7, 200_001);
+                        }
+                        Atan => {
+                            dist!(1e-3, 1e+7, 100_001);
+                            dist!(1e-2, 1e+8, 100_001);
+                            dist!(1e-3, 1e+7, 100_001);
+                            dist!(-1e-2, 1e+8, 100_001);
+                        }
+                        Atan2 => {
+                            dist!(-10, 10, 200_001);
+                        }
+                        Sinh | Cosh | Tanh | Asinh | Atanh => {
+                            if double {
+                                dist!(-700, 700, 200_001);
+                            } else {
+                                dist!(-88, 88, 200_001)
+                            }
+                        }
+
+                        Acosh => {
+                            if double {
+                                dist!(1, 700, 200_001);
+                            } else {
+                                dist!(1, 88, 200_001);
+                            }
+                        }
+                        Lgamma => {
+                            dist!(-5000, 5000, 200_001);
+                        }
+                        Tgamma => {
+                            dist!(-10, 10, 200_001);
+                        }
+                        Erf => {
+                            dist!(-100, 100, 200_001);
+                        }
+                        Erfc => {
+                            dist!(-1, 100, 200_001);
+                        }
+                        Fabs => {
+                            dist!(-100.5, 100.5, 200_001);
+                        }
+                        Copysign | Fmax | Fmin | Fmod | Nextafter => {
+                            dist!(-1e+10, 1e+10, 200_001);
+                        }
+                        Modf => {
+                            dist!(-1e+14, 1e+14, 200_001);
+                        }
+                        Trunc | Floor | Ceil | Round | Rint => {
+                            dist!(-100, 100, 800);
+                        }
+                        _ => (),
+                    }
+                }
+
+                // ~NSTEPS * 4
+                let current_len = vec.len();
+                assert!(len > current_len);
+                let remaining_len = len.checked_sub(current_len).unwrap();
+
+                for _ in 0..remaining_len {
+                    let n = rng.gen::<Self>();
+                    vec.push(n);
+                }
+                assert_eq!(vec.len(), len);
+
+                // Duplicate the vector, randomly shuffle it, and
+                // concatenate it. Otherwise for n-ary functions
+                // all vectors might have the same values. But
+                // testing with the same values is also worth doing.
+                let mut vec2 = vec.clone();
+                vec2.shuffle(rng);
+                vec.extend(vec2);
+                vec
+            }
+        }
+    };
+}
+
+impl_rand_seq_f!(f32);
+impl_rand_seq_f!(f64);
+
+impl RandSeq for i32 {
+    fn rand_seq<R: rand::Rng>(rng: &mut R, api_kind: ApiKind, len: usize) -> Vec<Self> {
+        use rand::seq::SliceRandom;
+
+        let mut v = Vec::with_capacity(len);
+        for _ in 0..len {
+            let mut r = rng.gen::<Self>();
+            if let ApiKind::Jn = api_kind {
+                // The integer argument of these APIs is a number of iterations.
+                // Computational cost blows up if we pass huge values, so zero
+                // their lower bits.
+                r &= 0xffff;
+            }
+            v.push(r);
+        }
+        assert_eq!(v.len(), len);
+
+        // Duplicate the vector, randomly shuffle it, and
+        // concatenate it. Otherwise for n-ary functions
+        // all vectors might have the same values. But
+        // testing with the same values is also worth doing.
+        let mut v2 = v.clone();
+        v2.shuffle(rng);
+        v.extend(v2);
+        v
+    }
+}
diff --git a/crates/libm-test/src/tuple_vec.rs b/crates/libm-test/src/tuple_vec.rs
new file mode 100644
index 000000000..dbfaeecce
--- /dev/null
+++ b/crates/libm-test/src/tuple_vec.rs
@@ -0,0 +1,31 @@
+//! Trait for tuples of vectors
+
+/// Access to tuple elements of a tuple of vectors
+pub trait TupleVec {
+    type Output;
+    fn get(&self, i: usize) -> Self::Output;
+}
+
+macro_rules! impl_tuple_vec {
+    (($($arg_tys:ty),*): $self_:ident: $($xs:expr),*)  => {
+        impl TupleVec for ($(Vec<$arg_tys>,)+) {
+            type Output = ($($arg_tys,)+);
+            fn get(&self, i: usize) -> Self::Output {
+                let $self_ = self;
+                ($($xs[i],)*)
+            }
+        }
+    };
+}
+
+impl_tuple_vec!((f32): x: x.0);
+impl_tuple_vec!((f64): x: x.0);
+impl_tuple_vec!((f32, f32): x: x.0, x.1);
+impl_tuple_vec!((f32, f64): x: x.0, x.1);
+impl_tuple_vec!((f64, f64): x: x.0, x.1);
+impl_tuple_vec!((f64, i32): x: x.0, x.1);
+impl_tuple_vec!((f32, i32): x: x.0, x.1);
+impl_tuple_vec!((i32, f64): x: x.0, x.1);
+impl_tuple_vec!((i32, f32): x: x.0, x.1);
+impl_tuple_vec!((f32, f32, f32): x: x.0, x.1, x.2);
+impl_tuple_vec!((f64, f64, f64): x: x.0, x.1, x.2);
diff --git a/crates/libm-test/src/within_ulps.rs b/crates/libm-test/src/within_ulps.rs
new file mode 100644
index 000000000..88761852e
--- /dev/null
+++ b/crates/libm-test/src/within_ulps.rs
@@ -0,0 +1,50 @@
+//! Floating-point comparison with ULP tolerance
+
+pub trait WithinUlps {
+    /// Returns true if two numbers are closer than `ulp_tol` to each other.
+    fn within_ulps(self, other: Self, ulp_tol: usize) -> bool;
+}
+
+// Stamp the impls for floats:
+macro_rules! impl_within_ulps_f {
+    ($f_ty:ty, $i_ty:ty) => {
+        impl WithinUlps for $f_ty {
+            #[allow(clippy::cast_possible_truncation)]
+            fn within_ulps(self, y: Self, ulp_tol: usize) -> bool {
+                let x = self;
+                if x.is_nan() != y.is_nan() {
+                    // one is nan but the other is not
+                    return false;
+                }
+                if x.is_nan() && y.is_nan() {
+                    return true;
+                }
+                if x.is_infinite() != y.is_infinite() {
+                    // one is inf but the other is not
+                    return false;
+                }
+
+                let xi: $i_ty = unsafe { core::intrinsics::transmute(x) };
+                let yi: $i_ty = unsafe { core::intrinsics::transmute(y) };
+                if (xi < 0) != (yi < 0) {
+                    // different sign, e.g., -0.0 != +0.0:
+                    return false;
+                }
+                let ulps = (xi - yi).abs();
+                ulps <= ulp_tol as _
+            }
+        }
+    };
+}
+
+impl_within_ulps_f!(f32, i32);
+impl_within_ulps_f!(f64, i64);
+
+impl WithinUlps for i32 {
+    #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
+    fn within_ulps(self, y: Self, ulp_tol: usize) -> bool {
+        let x = self;
+        let ulps = (x - y).abs();
+        ulps <= ulp_tol as _
+    }
+}
diff --git a/crates/libm-test/tests/exhaustive32.rs b/crates/libm-test/tests/exhaustive32.rs
new file mode 100644
index 000000000..28bb7b43a
--- /dev/null
+++ b/crates/libm-test/tests/exhaustive32.rs
@@ -0,0 +1,46 @@
+//! Exhaustively test unary APIs taking 32-bit wide arguments.
+#![cfg(test)]
+#![cfg(all(exhaustive32, feature = "exhaustive"))]
+use libm_test::assert_approx_eq;
+
+macro_rules! exhaustive32 {
+    (
+        id: $id:ident;
+        api_kind: $api_kind:ident;
+        arg_tys: f32;
+        arg_ids: $arg_id:ident;
+        ret_ty: $ret_ty:ty;
+    ) => {
+        #[test]
+        #[allow(unused)]
+        fn $id() {
+            extern "C" {
+                // The system's libm function:
+                fn $id($arg_id: f32) -> $ret_ty;
+            }
+
+            for i in 0..=u32::max_value() {
+                let arg: f32 = unsafe { std::mem::transmute(i) };
+                let result = libm::$id(arg);
+                let expected = unsafe { $id(arg) };
+                assert_approx_eq!(
+                    result == expected,
+                    id: $id, arg: arg, ulp: 4
+                );
+            }
+        }
+    };
+    (
+        id: $id:ident;
+        api_kind: $api_kind:ident;
+        arg_tys: $($arg_tys:ty),*;
+        arg_ids: $($arg_ids:ident),*;
+        ret_ty: $ret_ty:ty;
+    ) => {};
+}
+
+libm_analyze::for_each_api!(exhaustive32(
+    // jn is to expensive - the others have weird API:
+    /*ignored:*/
+    "j0f,j1f,y0f,y1f,exp10f,jn"
+));
diff --git a/crates/libm-test/tests/system_libm.rs b/crates/libm-test/tests/system_libm.rs
new file mode 100644
index 000000000..e98d577b7
--- /dev/null
+++ b/crates/libm-test/tests/system_libm.rs
@@ -0,0 +1,61 @@
+//! Compare the results of the `libm` implementation against the system's libm.
+#![cfg(test)]
+#![cfg(feature = "system_libm")]
+
+use libm_test::{assert_approx_eq, CallFn, RandSeq, TupleVec};
+
+// Number of tests to generate for each function
+const NTESTS: usize = 1_000_000;
+
+const ULP_TOL: usize = 4;
+
+macro_rules! system_libm {
+    // Generate random tests for all others:
+    (
+        id: $id:ident;
+        api_kind: $api_kind:ident;
+        arg_tys: $($arg_tys:ty),*;
+        arg_ids: $($arg_ids:ident),*;
+        ret_ty: $ret_ty:ty;
+    ) => {
+        #[test]
+        #[allow(unused)]
+        fn $id() {
+            // Type of the system libm fn:
+            type FnTy
+                = unsafe extern "C" fn ($($arg_ids: $arg_tys),*) -> $ret_ty;
+            extern "C" {
+                // The system's libm function:
+                fn $id($($arg_ids: $arg_tys),*) -> $ret_ty;
+            }
+
+            let mut rng = rand::thread_rng();
+
+            // Generate a tuple of arguments containing random values:
+            let mut args: ( $(Vec<$arg_tys>,)+ )
+                = ( $(
+                        <$arg_tys as RandSeq>::rand_seq(
+                            &mut rng, libm_test::ApiKind::$api_kind, NTESTS
+                         ),
+                     )+
+                );
+
+            for i in 0..NTESTS {
+                let args: ( $($arg_tys,)+ ) = args.get(i);
+                let result = args.call(libm::$id as FnTy);
+                let expected = args.call($id as FnTy);
+                assert_approx_eq!(
+                    result == expected,
+                    id: $id, arg: args, ulp: ULP_TOL
+                );
+            }
+        }
+    }
+}
+
+libm_analyze::for_each_api!(system_libm(
+    // FIXME: Some are not exposed by the system's musl,
+    // others are incorrect. FMA is broken.
+    /*ignored:*/
+    "j0f,j1f,jnf,y0f,y1f,ynf,exp10,exp10f,fma"
+));
diff --git a/crates/libm-test/tests/unit.rs b/crates/libm-test/tests/unit.rs
new file mode 100644
index 000000000..2915b23da
--- /dev/null
+++ b/crates/libm-test/tests/unit.rs
@@ -0,0 +1,369 @@
+#![cfg(test)]
+
+use libm::*;
+use libm_test::assert_approx_eq;
+
+#[test]
+fn remquo_q_overflow() {
+    // 0xc000000000000001, 0x04c0000000000004
+    let _ = remquo(-2.0000000000000004, 8.406091369059082e-286);
+}
+
+#[test]
+fn ceil_sanity_check() {
+    assert_eq!(ceil(1.1), 2.0);
+    assert_eq!(ceil(2.9), 3.0);
+}
+
+#[test]
+fn atan2_sanity_check() {
+    use std::f64::consts::PI;
+    assert_eq!(atan2(0.0, 1.0), 0.0);
+    assert_eq!(atan2(0.0, -1.0), PI);
+    assert_eq!(atan2(-0.0, -1.0), -PI);
+    assert_eq!(atan2(3.0, 2.0), atan(3.0 / 2.0));
+    assert_eq!(atan2(2.0, -1.0), atan(2.0 / -1.0) + PI);
+    assert_eq!(atan2(-2.0, -1.0), atan(-2.0 / -1.0) - PI);
+}
+
+#[test]
+fn floor_overflow() {
+    assert_eq!(floorf(0.5), 0.0);
+}
+
+#[test]
+fn trunc_sanity_check() {
+    assert_eq!(trunc(1.1), 1.0);
+}
+
+mod pow {
+    use libm::*;
+    use std::f64::consts::{E, PI};
+    use std::f64::{EPSILON, INFINITY, MAX, MIN, MIN_POSITIVE, NAN, NEG_INFINITY};
+
+    const POS_ZERO: &[f64] = &[0.0];
+    const NEG_ZERO: &[f64] = &[-0.0];
+    const POS_ONE: &[f64] = &[1.0];
+    const NEG_ONE: &[f64] = &[-1.0];
+    const POS_FLOATS: &[f64] = &[99.0 / 70.0, E, PI];
+    const NEG_FLOATS: &[f64] = &[-99.0 / 70.0, -E, -PI];
+    const POS_SMALL_FLOATS: &[f64] = &[(1.0 / 2.0), MIN_POSITIVE, EPSILON];
+    const NEG_SMALL_FLOATS: &[f64] = &[-(1.0 / 2.0), -MIN_POSITIVE, -EPSILON];
+    const POS_EVENS: &[f64] = &[2.0, 6.0, 8.0, 10.0, 22.0, 100.0, MAX];
+    const NEG_EVENS: &[f64] = &[MIN, -100.0, -22.0, -10.0, -8.0, -6.0, -2.0];
+    const POS_ODDS: &[f64] = &[3.0, 7.0];
+    const NEG_ODDS: &[f64] = &[-7.0, -3.0];
+    const NANS: &[f64] = &[NAN];
+    const POS_INF: &[f64] = &[INFINITY];
+    const NEG_INF: &[f64] = &[NEG_INFINITY];
+
+    const ALL: &[&[f64]] = &[
+        POS_ZERO,
+        NEG_ZERO,
+        NANS,
+        NEG_SMALL_FLOATS,
+        POS_SMALL_FLOATS,
+        NEG_FLOATS,
+        POS_FLOATS,
+        NEG_EVENS,
+        POS_EVENS,
+        NEG_ODDS,
+        POS_ODDS,
+        NEG_INF,
+        POS_INF,
+        NEG_ONE,
+        POS_ONE,
+    ];
+    const POS: &[&[f64]] = &[POS_ZERO, POS_ODDS, POS_ONE, POS_FLOATS, POS_EVENS, POS_INF];
+    const NEG: &[&[f64]] = &[NEG_ZERO, NEG_ODDS, NEG_ONE, NEG_FLOATS, NEG_EVENS, NEG_INF];
+
+    fn pow_test(base: f64, exponent: f64, expected: f64) {
+        let res = pow(base, exponent);
+        assert!(
+            if expected.is_nan() {
+                res.is_nan()
+            } else {
+                pow(base, exponent) == expected
+            },
+            "{} ** {} was {} instead of {}",
+            base,
+            exponent,
+            res,
+            expected
+        );
+    }
+
+    fn test_sets_as_base(sets: &[&[f64]], exponent: f64, expected: f64) {
+        sets.iter()
+            .for_each(|s| s.iter().for_each(|val| pow_test(*val, exponent, expected)));
+    }
+
+    fn test_sets_as_exponent(base: f64, sets: &[&[f64]], expected: f64) {
+        sets.iter()
+            .for_each(|s| s.iter().for_each(|val| pow_test(base, *val, expected)));
+    }
+
+    fn test_sets(sets: &[&[f64]], computed: &dyn Fn(f64) -> f64, expected: &dyn Fn(f64) -> f64) {
+        sets.iter().for_each(|s| {
+            s.iter().for_each(|val| {
+                let exp = expected(*val);
+                let res = computed(*val);
+
+                assert!(
+                    if exp.is_nan() {
+                        res.is_nan()
+                    } else {
+                        exp == res
+                    },
+                    "test for {} was {} instead of {}",
+                    val,
+                    res,
+                    exp
+                );
+            })
+        });
+    }
+
+    #[test]
+    fn zero_as_exponent() {
+        test_sets_as_base(ALL, 0.0, 1.0);
+        test_sets_as_base(ALL, -0.0, 1.0);
+    }
+
+    #[test]
+    fn one_as_base() {
+        test_sets_as_exponent(1.0, ALL, 1.0);
+    }
+
+    #[test]
+    fn nan_inputs() {
+        // NAN as the base:
+        // (NAN ^ anything *but 0* should be NAN)
+        test_sets_as_exponent(NAN, &ALL[2..], NAN);
+
+        // NAN as the exponent:
+        // (anything *but 1* ^ NAN should be NAN)
+        test_sets_as_base(&ALL[..(ALL.len() - 2)], NAN, NAN);
+    }
+
+    #[test]
+    fn infinity_as_base() {
+        // Positive Infinity as the base:
+        // (+Infinity ^ positive anything but 0 and NAN should be +Infinity)
+        test_sets_as_exponent(INFINITY, &POS[1..], INFINITY);
+
+        // (+Infinity ^ negative anything except 0 and NAN should be 0.0)
+        test_sets_as_exponent(INFINITY, &NEG[1..], 0.0);
+
+        // Negative Infinity as the base:
+        // (-Infinity ^ positive odd ints should be -Infinity)
+        test_sets_as_exponent(NEG_INFINITY, &[POS_ODDS], NEG_INFINITY);
+
+        // (-Infinity ^ anything but odd ints should be == -0 ^ (-anything))
+        // We can lump in pos/neg odd ints here because they don't seem to
+        // cause panics (div by zero) in release mode (I think).
+        test_sets(ALL, &|v: f64| pow(NEG_INFINITY, v), &|v: f64| pow(-0.0, -v));
+    }
+
+    #[test]
+    fn infinity_as_exponent() {
+        // Positive/Negative base greater than 1:
+        // (pos/neg > 1 ^ Infinity should be Infinity - note this excludes NAN as the base)
+        test_sets_as_base(&ALL[5..(ALL.len() - 2)], INFINITY, INFINITY);
+
+        // (pos/neg > 1 ^ -Infinity should be 0.0)
+        test_sets_as_base(&ALL[5..ALL.len() - 2], NEG_INFINITY, 0.0);
+
+        // Positive/Negative base less than 1:
+        let base_below_one = &[POS_ZERO, NEG_ZERO, NEG_SMALL_FLOATS, POS_SMALL_FLOATS];
+
+        // (pos/neg < 1 ^ Infinity should be 0.0 - this also excludes NAN as the base)
+        test_sets_as_base(base_below_one, INFINITY, 0.0);
+
+        // (pos/neg < 1 ^ -Infinity should be Infinity)
+        test_sets_as_base(base_below_one, NEG_INFINITY, INFINITY);
+
+        // Positive/Negative 1 as the base:
+        // (pos/neg 1 ^ Infinity should be 1)
+        test_sets_as_base(&[NEG_ONE, POS_ONE], INFINITY, 1.0);
+
+        // (pos/neg 1 ^ -Infinity should be 1)
+        test_sets_as_base(&[NEG_ONE, POS_ONE], NEG_INFINITY, 1.0);
+    }
+
+    #[test]
+    fn zero_as_base() {
+        // Positive Zero as the base:
+        // (+0 ^ anything positive but 0 and NAN should be +0)
+        test_sets_as_exponent(0.0, &POS[1..], 0.0);
+
+        // (+0 ^ anything negative but 0 and NAN should be Infinity)
+        // (this should panic because we're dividing by zero)
+        test_sets_as_exponent(0.0, &NEG[1..], INFINITY);
+
+        // Negative Zero as the base:
+        // (-0 ^ anything positive but 0, NAN, and odd ints should be +0)
+        test_sets_as_exponent(-0.0, &POS[3..], 0.0);
+
+        // (-0 ^ anything negative but 0, NAN, and odd ints should be Infinity)
+        // (should panic because of divide by zero)
+        test_sets_as_exponent(-0.0, &NEG[3..], INFINITY);
+
+        // (-0 ^ positive odd ints should be -0)
+        test_sets_as_exponent(-0.0, &[POS_ODDS], -0.0);
+
+        // (-0 ^ negative odd ints should be -Infinity)
+        // (should panic because of divide by zero)
+        test_sets_as_exponent(-0.0, &[NEG_ODDS], NEG_INFINITY);
+    }
+
+    #[test]
+    fn special_cases() {
+        // One as the exponent:
+        // (anything ^ 1 should be anything - i.e. the base)
+        test_sets(ALL, &|v: f64| pow(v, 1.0), &|v: f64| v);
+
+        // Negative One as the exponent:
+        // (anything ^ -1 should be 1/anything)
+        test_sets(ALL, &|v: f64| pow(v, -1.0), &|v: f64| 1.0 / v);
+
+        // Factoring -1 out:
+        // (negative anything ^ integer should be (-1 ^ integer) * (positive anything ^ integer))
+        &[POS_ZERO, NEG_ZERO, POS_ONE, NEG_ONE, POS_EVENS, NEG_EVENS]
+            .iter()
+            .for_each(|int_set| {
+                int_set.iter().for_each(|int| {
+                    test_sets(ALL, &|v: f64| pow(-v, *int), &|v: f64| {
+                        pow(-1.0, *int) * pow(v, *int)
+                    });
+                })
+            });
+
+        // Negative base (imaginary results):
+        // (-anything except 0 and Infinity ^ non-integer should be NAN)
+        &NEG[1..(NEG.len() - 1)].iter().for_each(|set| {
+            set.iter().for_each(|val| {
+                test_sets(&ALL[3..7], &|v: f64| pow(*val, v), &|_| NAN);
+            })
+        });
+    }
+
+    #[test]
+    fn normal_cases() {
+        assert_eq!(pow(2.0, 20.0), (1 << 20) as f64);
+        assert_eq!(pow(-1.0, 9.0), -1.0);
+        assert!(pow(-1.0, 2.2).is_nan());
+        assert!(pow(-1.0, -1.14).is_nan());
+    }
+}
+
+mod atan {
+    use libm::atan;
+    use std::f64;
+
+    #[test]
+    fn sanity_check() {
+        for (input, answer) in [
+            (3.0_f64.sqrt() / 3.0, f64::consts::FRAC_PI_6),
+            (1.0, f64::consts::FRAC_PI_4),
+            (3.0_f64.sqrt(), f64::consts::FRAC_PI_3),
+            (-3.0_f64.sqrt() / 3.0, -f64::consts::FRAC_PI_6),
+            (-1.0, -f64::consts::FRAC_PI_4),
+            (-3.0_f64.sqrt(), -f64::consts::FRAC_PI_3),
+        ]
+        .iter()
+        {
+            assert!(
+                (atan(*input) - answer) / answer < 1e-5,
+                "\natan({:.4}/16) = {:.4}, actual: {}",
+                input * 16.0,
+                answer,
+                atan(*input)
+            );
+        }
+    }
+
+    #[test]
+    fn zero() {
+        assert_eq!(atan(0.0), 0.0);
+    }
+
+    #[test]
+    fn infinity() {
+        assert_eq!(atan(f64::INFINITY), f64::consts::FRAC_PI_2);
+    }
+
+    #[test]
+    fn minus_infinity() {
+        assert_eq!(atan(f64::NEG_INFINITY), -f64::consts::FRAC_PI_2);
+    }
+
+    #[test]
+    fn nan() {
+        assert!(atan(f64::NAN).is_nan());
+    }
+}
+
+#[test]
+fn sin_near_pi() {
+    let x = f64::from_bits(0x400921fb000FD5DD); // 3.141592026217707
+    let sx = f64::from_bits(0x3ea50d15ced1a4a2); // 6.273720864039205e-7
+    assert_eq!(sin(x), sx);
+}
+
+#[test]
+fn truncf_sanity_check() {
+    assert_eq!(truncf(1.1), 1.0);
+}
+
+#[test]
+fn expm1_sanity_check() {
+    assert_eq!(expm1(1.1), 2.0041660239464334);
+}
+
+#[test]
+fn roundf_negative_zero() {
+    assert_eq!(roundf(-0.0_f32).to_bits(), (-0.0_f32).to_bits());
+}
+
+#[test]
+fn exp2_i0_wrap_test() {
+    let x = -3.0 / 256.0;
+    assert_eq!(exp2(x), f64::from_bits(0x3fefbdba3692d514));
+}
+
+#[test]
+fn round_negative_zero() {
+    assert_eq!(round(-0.0_f64).to_bits(), (-0.0_f64).to_bits());
+}
+
+#[test]
+fn j1f_2488() {
+    // 0x401F3E49
+    assert_eq!(j1f(2.4881766_f32), 0.49999475_f32);
+}
+#[test]
+fn y1f_2002() {
+    assert_eq!(y1f(2.0000002_f32), -0.10703229_f32);
+}
+
+#[test]
+fn fma_segfault_bug() {
+    // An attempt to substract with overflow was causing a segfault
+    // on FMA for these inputs:
+    assert_approx_eq!(
+        fma(
+            -0.0000000000000002220446049250313,
+            -0.0000000000000002220446049250313,
+            -0.0000000000000002220446049250313
+        ),
+        -0.00000000000000022204460492503126,
+        ulp: 0
+    );
+
+    assert_approx_eq!(
+        fma(-0.992, -0.992, -0.992),
+        -0.00793599999988632,
+        ulp: 0
+    );
+}
diff --git a/crates/libm/Cargo.toml b/crates/libm/Cargo.toml
new file mode 100644
index 000000000..40cbd2b33
--- /dev/null
+++ b/crates/libm/Cargo.toml
@@ -0,0 +1,26 @@
+[package]
+authors = ["Jorge Aparicio <jorge@japaric.io>"]
+categories = ["no-std"]
+description = "libm in pure Rust"
+documentation = "https://docs.rs/libm"
+keywords = ["libm", "math"]
+license = "MIT OR Apache-2.0"
+name = "libm"
+repository = "https://github.com/rust-lang-nursery/libm"
+version = "0.2.0"
+edition = "2018"
+
+[features]
+# only used to run our test suite
+default = ['stable']
+stable = []
+
+# Used checked array indexing instead of unchecked array indexing in this
+# library.
+checked = []
+
+[dev-dependencies]
+no-panic = "0.1.8"
+
+[build-dependencies]
+rand = { version = "0.6.5", optional = true }
diff --git a/crates/libm/src/lib.rs b/crates/libm/src/lib.rs
new file mode 100644
index 000000000..527fcb4f9
--- /dev/null
+++ b/crates/libm/src/lib.rs
@@ -0,0 +1,10 @@
+//! libm in pure Rust
+#![deny(warnings)]
+#![no_std]
+#![cfg_attr(
+    all(target_arch = "wasm32", not(feature = "stable")),
+    feature(core_intrinsics)
+)]
+
+mod math;
+pub use self::math::*;
diff --git a/src/math/acos.rs b/crates/libm/src/math/acos.rs
similarity index 98%
rename from src/math/acos.rs
rename to crates/libm/src/math/acos.rs
index d5e1f6865..2d98a9948 100644
--- a/src/math/acos.rs
+++ b/crates/libm/src/math/acos.rs
@@ -62,7 +62,7 @@ fn r(z: f64) -> f64 {
 /// Returns values in radians, in the range of 0 to pi.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn acos(x: f64) -> f64 {
+pub extern "C" fn acos(x: f64) -> f64 {
     let x1p_120f = f64::from_bits(0x3870000000000000); // 0x1p-120 === 2 ^ -120
     let z: f64;
     let w: f64;
diff --git a/src/math/acosf.rs b/crates/libm/src/math/acosf.rs
similarity index 98%
rename from src/math/acosf.rs
rename to crates/libm/src/math/acosf.rs
index d0598e811..2aab1673c 100644
--- a/src/math/acosf.rs
+++ b/crates/libm/src/math/acosf.rs
@@ -36,7 +36,7 @@ fn r(z: f32) -> f32 {
 /// Returns values in radians, in the range of 0 to pi.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn acosf(x: f32) -> f32 {
+pub extern "C" fn acosf(x: f32) -> f32 {
     let x1p_120 = f32::from_bits(0x03800000); // 0x1p-120 === 2 ^ (-120)
 
     let z: f32;
diff --git a/src/math/acosh.rs b/crates/libm/src/math/acosh.rs
similarity index 87%
rename from src/math/acosh.rs
rename to crates/libm/src/math/acosh.rs
index ac7a5f1c6..344dd8700 100644
--- a/src/math/acosh.rs
+++ b/crates/libm/src/math/acosh.rs
@@ -7,7 +7,9 @@ const LN2: f64 = 0.693147180559945309417232121458176568; /* 0x3fe62e42,  0xfefa3
 /// Calculates the inverse hyperbolic cosine of `x`.
 /// Is defined as `log(x + sqrt(x*x-1))`.
 /// `x` must be a number greater than or equal to 1.
-pub fn acosh(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn acosh(x: f64) -> f64 {
     let u = x.to_bits();
     let e = ((u >> 52) as usize) & 0x7ff;
 
diff --git a/src/math/acoshf.rs b/crates/libm/src/math/acoshf.rs
similarity index 86%
rename from src/math/acoshf.rs
rename to crates/libm/src/math/acoshf.rs
index 0879e1edb..bc6af68f8 100644
--- a/src/math/acoshf.rs
+++ b/crates/libm/src/math/acoshf.rs
@@ -7,7 +7,9 @@ const LN2: f32 = 0.693147180559945309417232121458176568;
 /// Calculates the inverse hyperbolic cosine of `x`.
 /// Is defined as `log(x + sqrt(x*x-1))`.
 /// `x` must be a number greater than or equal to 1.
-pub fn acoshf(x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn acoshf(x: f32) -> f32 {
     let u = x.to_bits();
     let a = u & 0x7fffffff;
 
diff --git a/src/math/asin.rs b/crates/libm/src/math/asin.rs
similarity index 98%
rename from src/math/asin.rs
rename to crates/libm/src/math/asin.rs
index 774475e51..31bbcde46 100644
--- a/src/math/asin.rs
+++ b/crates/libm/src/math/asin.rs
@@ -69,7 +69,7 @@ fn comp_r(z: f64) -> f64 {
 /// Returns values in radians, in the range of -pi/2 to pi/2.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn asin(mut x: f64) -> f64 {
+pub extern "C" fn asin(mut x: f64) -> f64 {
     let z: f64;
     let r: f64;
     let s: f64;
diff --git a/src/math/asinf.rs b/crates/libm/src/math/asinf.rs
similarity index 97%
rename from src/math/asinf.rs
rename to crates/libm/src/math/asinf.rs
index ce0f4a997..b0bb0b4c4 100644
--- a/src/math/asinf.rs
+++ b/crates/libm/src/math/asinf.rs
@@ -38,7 +38,7 @@ fn r(z: f32) -> f32 {
 /// Returns values in radians, in the range of -pi/2 to pi/2.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn asinf(mut x: f32) -> f32 {
+pub extern "C" fn asinf(mut x: f32) -> f32 {
     let x1p_120 = f64::from_bits(0x3870000000000000); // 0x1p-120 === 2 ^ (-120)
 
     let hx = x.to_bits();
diff --git a/src/math/asinh.rs b/crates/libm/src/math/asinh.rs
similarity index 90%
rename from src/math/asinh.rs
rename to crates/libm/src/math/asinh.rs
index 14295357a..7c0f84f04 100644
--- a/src/math/asinh.rs
+++ b/crates/libm/src/math/asinh.rs
@@ -7,7 +7,9 @@ const LN2: f64 = 0.693147180559945309417232121458176568; /* 0x3fe62e42,  0xfefa3
 ///
 /// Calculates the inverse hyperbolic sine of `x`.
 /// Is defined as `sgn(x)*log(|x|+sqrt(x*x+1))`.
-pub fn asinh(mut x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn asinh(mut x: f64) -> f64 {
     let mut u = x.to_bits();
     let e = ((u >> 52) as usize) & 0x7ff;
     let sign = (u >> 63) != 0;
diff --git a/src/math/asinhf.rs b/crates/libm/src/math/asinhf.rs
similarity index 89%
rename from src/math/asinhf.rs
rename to crates/libm/src/math/asinhf.rs
index e22a29132..6072d3751 100644
--- a/src/math/asinhf.rs
+++ b/crates/libm/src/math/asinhf.rs
@@ -7,7 +7,9 @@ const LN2: f32 = 0.693147180559945309417232121458176568;
 ///
 /// Calculates the inverse hyperbolic sine of `x`.
 /// Is defined as `sgn(x)*log(|x|+sqrt(x*x+1))`.
-pub fn asinhf(mut x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn asinhf(mut x: f32) -> f32 {
     let u = x.to_bits();
     let i = u & 0x7fffffff;
     let sign = (u >> 31) != 0;
diff --git a/src/math/atan.rs b/crates/libm/src/math/atan.rs
similarity index 80%
rename from src/math/atan.rs
rename to crates/libm/src/math/atan.rs
index d2684ece8..57244146a 100644
--- a/src/math/atan.rs
+++ b/crates/libm/src/math/atan.rs
@@ -66,7 +66,7 @@ const AT: [f64; 11] = [
 /// Returns a value in radians, in the range of -pi/2 to pi/2.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn atan(x: f64) -> f64 {
+pub extern "C" fn atan(x: f64) -> f64 {
     let mut x = x;
     let mut ix = (x.to_bits() >> 32) as u32;
     let sign = ix >> 31;
@@ -135,51 +135,3 @@ pub fn atan(x: f64) -> f64 {
         z
     }
 }
-
-#[cfg(test)]
-mod tests {
-    use super::atan;
-    use core::f64;
-
-    #[test]
-    fn sanity_check() {
-        for (input, answer) in [
-            (3.0_f64.sqrt() / 3.0, f64::consts::FRAC_PI_6),
-            (1.0, f64::consts::FRAC_PI_4),
-            (3.0_f64.sqrt(), f64::consts::FRAC_PI_3),
-            (-3.0_f64.sqrt() / 3.0, -f64::consts::FRAC_PI_6),
-            (-1.0, -f64::consts::FRAC_PI_4),
-            (-3.0_f64.sqrt(), -f64::consts::FRAC_PI_3),
-        ]
-        .iter()
-        {
-            assert!(
-                (atan(*input) - answer) / answer < 1e-5,
-                "\natan({:.4}/16) = {:.4}, actual: {}",
-                input * 16.0,
-                answer,
-                atan(*input)
-            );
-        }
-    }
-
-    #[test]
-    fn zero() {
-        assert_eq!(atan(0.0), 0.0);
-    }
-
-    #[test]
-    fn infinity() {
-        assert_eq!(atan(f64::INFINITY), f64::consts::FRAC_PI_2);
-    }
-
-    #[test]
-    fn minus_infinity() {
-        assert_eq!(atan(f64::NEG_INFINITY), -f64::consts::FRAC_PI_2);
-    }
-
-    #[test]
-    fn nan() {
-        assert!(atan(f64::NAN).is_nan());
-    }
-}
diff --git a/src/math/atan2.rs b/crates/libm/src/math/atan2.rs
similarity index 91%
rename from src/math/atan2.rs
rename to crates/libm/src/math/atan2.rs
index 08385cd10..feef73ea0 100644
--- a/src/math/atan2.rs
+++ b/crates/libm/src/math/atan2.rs
@@ -50,7 +50,7 @@ const PI_LO: f64 = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
 /// Returns a value in radians, in the range of -pi to pi.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn atan2(y: f64, x: f64) -> f64 {
+pub extern "C" fn atan2(y: f64, x: f64) -> f64 {
     if x.is_nan() || y.is_nan() {
         return x + y;
     }
@@ -115,13 +115,3 @@ pub fn atan2(y: f64, x: f64) -> f64 {
         _ => (z - PI_LO) - PI, /* atan(-,-) */
     }
 }
-
-#[test]
-fn sanity_check() {
-    assert_eq!(atan2(0.0, 1.0), 0.0);
-    assert_eq!(atan2(0.0, -1.0), PI);
-    assert_eq!(atan2(-0.0, -1.0), -PI);
-    assert_eq!(atan2(3.0, 2.0), atan(3.0 / 2.0));
-    assert_eq!(atan2(2.0, -1.0), atan(2.0 / -1.0) + PI);
-    assert_eq!(atan2(-2.0, -1.0), atan(-2.0 / -1.0) - PI);
-}
diff --git a/src/math/atan2f.rs b/crates/libm/src/math/atan2f.rs
similarity index 98%
rename from src/math/atan2f.rs
rename to crates/libm/src/math/atan2f.rs
index 7bbe5f1d4..e5b169d39 100644
--- a/src/math/atan2f.rs
+++ b/crates/libm/src/math/atan2f.rs
@@ -26,7 +26,7 @@ const PI_LO: f32 = -8.7422776573e-08; /* 0xb3bbbd2e */
 /// Returns a value in radians, in the range of -pi to pi.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn atan2f(y: f32, x: f32) -> f32 {
+pub extern "C" fn atan2f(y: f32, x: f32) -> f32 {
     if x.is_nan() || y.is_nan() {
         return x + y;
     }
diff --git a/src/math/atanf.rs b/crates/libm/src/math/atanf.rs
similarity index 98%
rename from src/math/atanf.rs
rename to crates/libm/src/math/atanf.rs
index 363e11d64..e9250b468 100644
--- a/src/math/atanf.rs
+++ b/crates/libm/src/math/atanf.rs
@@ -43,7 +43,7 @@ const A_T: [f32; 5] = [
 /// Returns a value in radians, in the range of -pi/2 to pi/2.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn atanf(mut x: f32) -> f32 {
+pub extern "C" fn atanf(mut x: f32) -> f32 {
     let x1p_120 = f32::from_bits(0x03800000); // 0x1p-120 === 2 ^ (-120)
 
     let z: f32;
diff --git a/src/math/atanh.rs b/crates/libm/src/math/atanh.rs
similarity index 88%
rename from src/math/atanh.rs
rename to crates/libm/src/math/atanh.rs
index 79a989c42..58d3c1eb2 100644
--- a/src/math/atanh.rs
+++ b/crates/libm/src/math/atanh.rs
@@ -5,7 +5,9 @@ use super::log1p;
 ///
 /// Calculates the inverse hyperbolic tangent of `x`.
 /// Is defined as `log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2`.
-pub fn atanh(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn atanh(x: f64) -> f64 {
     let u = x.to_bits();
     let e = ((u >> 52) as usize) & 0x7ff;
     let sign = (u >> 63) != 0;
diff --git a/src/math/atanhf.rs b/crates/libm/src/math/atanhf.rs
similarity index 88%
rename from src/math/atanhf.rs
rename to crates/libm/src/math/atanhf.rs
index 7b2f34d97..5e66e4196 100644
--- a/src/math/atanhf.rs
+++ b/crates/libm/src/math/atanhf.rs
@@ -5,7 +5,9 @@ use super::log1pf;
 ///
 /// Calculates the inverse hyperbolic tangent of `x`.
 /// Is defined as `log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2`.
-pub fn atanhf(mut x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn atanhf(mut x: f32) -> f32 {
     let mut u = x.to_bits();
     let sign = (u >> 31) != 0;
 
diff --git a/src/math/cbrt.rs b/crates/libm/src/math/cbrt.rs
similarity index 99%
rename from src/math/cbrt.rs
rename to crates/libm/src/math/cbrt.rs
index 04469b159..3f5253175 100644
--- a/src/math/cbrt.rs
+++ b/crates/libm/src/math/cbrt.rs
@@ -32,7 +32,7 @@ const P4: f64 = 0.145996192886612446982; /* 0x3fc2b000, 0xd4e4edd7 */
 /// Computes the cube root of the argument.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn cbrt(x: f64) -> f64 {
+pub extern "C" fn cbrt(x: f64) -> f64 {
     let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
 
     let mut ui: u64 = x.to_bits();
diff --git a/src/math/cbrtf.rs b/crates/libm/src/math/cbrtf.rs
similarity index 98%
rename from src/math/cbrtf.rs
rename to crates/libm/src/math/cbrtf.rs
index 6e589c099..a129ef9bc 100644
--- a/src/math/cbrtf.rs
+++ b/crates/libm/src/math/cbrtf.rs
@@ -27,7 +27,7 @@ const B2: u32 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
 /// Computes the cube root of the argument.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn cbrtf(x: f32) -> f32 {
+pub extern "C" fn cbrtf(x: f32) -> f32 {
     let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
 
     let mut r: f64;
diff --git a/src/math/ceil.rs b/crates/libm/src/math/ceil.rs
similarity index 82%
rename from src/math/ceil.rs
rename to crates/libm/src/math/ceil.rs
index 59883a8a7..a263968c8 100644
--- a/src/math/ceil.rs
+++ b/crates/libm/src/math/ceil.rs
@@ -7,12 +7,12 @@ const TOINT: f64 = 1. / f64::EPSILON;
 /// Finds the nearest integer greater than or equal to `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn ceil(x: f64) -> f64 {
+pub extern "C" fn ceil(x: f64) -> f64 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f64.ceil` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::ceilf64(x) }
         }
     }
@@ -40,12 +40,3 @@ pub fn ceil(x: f64) -> f64 {
         x + y
     }
 }
-
-#[cfg(test)]
-mod tests {
-    #[test]
-    fn sanity_check() {
-        assert_eq!(super::ceil(1.1), 2.0);
-        assert_eq!(super::ceil(2.9), 3.0);
-    }
-}
diff --git a/src/math/ceilf.rs b/crates/libm/src/math/ceilf.rs
similarity index 90%
rename from src/math/ceilf.rs
rename to crates/libm/src/math/ceilf.rs
index 151a4f210..82924a0fe 100644
--- a/src/math/ceilf.rs
+++ b/crates/libm/src/math/ceilf.rs
@@ -5,12 +5,12 @@ use core::f32;
 /// Finds the nearest integer greater than or equal to `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn ceilf(x: f32) -> f32 {
+pub extern "C" fn ceilf(x: f32) -> f32 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f32.ceil` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::ceilf32(x) }
         }
     }
diff --git a/src/math/copysign.rs b/crates/libm/src/math/copysign.rs
similarity index 71%
rename from src/math/copysign.rs
rename to crates/libm/src/math/copysign.rs
index 1527fb6ea..fe0a020b6 100644
--- a/src/math/copysign.rs
+++ b/crates/libm/src/math/copysign.rs
@@ -2,7 +2,9 @@
 ///
 /// Constructs a number with the magnitude (absolute value) of its
 /// first argument, `x`, and the sign of its second argument, `y`.
-pub fn copysign(x: f64, y: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn copysign(x: f64, y: f64) -> f64 {
     let mut ux = x.to_bits();
     let uy = y.to_bits();
     ux &= (!0) >> 1;
diff --git a/src/math/copysignf.rs b/crates/libm/src/math/copysignf.rs
similarity index 71%
rename from src/math/copysignf.rs
rename to crates/libm/src/math/copysignf.rs
index 35148561a..1ecafc204 100644
--- a/src/math/copysignf.rs
+++ b/crates/libm/src/math/copysignf.rs
@@ -2,7 +2,9 @@
 ///
 /// Constructs a number with the magnitude (absolute value) of its
 /// first argument, `x`, and the sign of its second argument, `y`.
-pub fn copysignf(x: f32, y: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn copysignf(x: f32, y: f32) -> f32 {
     let mut ux = x.to_bits();
     let uy = y.to_bits();
     ux &= 0x7fffffff;
diff --git a/src/math/cos.rs b/crates/libm/src/math/cos.rs
similarity index 98%
rename from src/math/cos.rs
rename to crates/libm/src/math/cos.rs
index fe5a89919..4eeb271ec 100644
--- a/src/math/cos.rs
+++ b/crates/libm/src/math/cos.rs
@@ -43,7 +43,7 @@ use super::{k_cos, k_sin, rem_pio2};
 //
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn cos(x: f64) -> f64 {
+pub extern "C" fn cos(x: f64) -> f64 {
     let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
 
     /* |x| ~< pi/4 */
diff --git a/src/math/cosf.rs b/crates/libm/src/math/cosf.rs
similarity index 98%
rename from src/math/cosf.rs
rename to crates/libm/src/math/cosf.rs
index 48d76c8ee..fc912948b 100644
--- a/src/math/cosf.rs
+++ b/crates/libm/src/math/cosf.rs
@@ -26,7 +26,7 @@ const C4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn cosf(x: f32) -> f32 {
+pub extern "C" fn cosf(x: f32) -> f32 {
     let x64 = x as f64;
 
     let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
diff --git a/src/math/cosh.rs b/crates/libm/src/math/cosh.rs
similarity index 95%
rename from src/math/cosh.rs
rename to crates/libm/src/math/cosh.rs
index bac875566..fb7d28198 100644
--- a/src/math/cosh.rs
+++ b/crates/libm/src/math/cosh.rs
@@ -9,7 +9,7 @@ use super::k_expo2;
 /// Angles are specified in radians.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn cosh(mut x: f64) -> f64 {
+pub extern "C" fn cosh(mut x: f64) -> f64 {
     /* |x| */
     let mut ix = x.to_bits();
     ix &= 0x7fffffffffffffff;
diff --git a/src/math/coshf.rs b/crates/libm/src/math/coshf.rs
similarity index 95%
rename from src/math/coshf.rs
rename to crates/libm/src/math/coshf.rs
index bf99e42f0..f27fbd46c 100644
--- a/src/math/coshf.rs
+++ b/crates/libm/src/math/coshf.rs
@@ -9,7 +9,7 @@ use super::k_expo2f;
 /// Angles are specified in radians.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn coshf(mut x: f32) -> f32 {
+pub extern "C" fn coshf(mut x: f32) -> f32 {
     let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
 
     /* |x| */
diff --git a/src/math/erf.rs b/crates/libm/src/math/erf.rs
similarity index 98%
rename from src/math/erf.rs
rename to crates/libm/src/math/erf.rs
index a2c617d34..a3bf71096 100644
--- a/src/math/erf.rs
+++ b/crates/libm/src/math/erf.rs
@@ -219,7 +219,9 @@ fn erfc2(ix: u32, mut x: f64) -> f64 {
 /// Calculates an approximation to the “error function”, which estimates
 /// the probability that an observation will fall within x standard
 /// deviations of the mean (assuming a normal distribution).
-pub fn erf(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn erf(x: f64) -> f64 {
     let r: f64;
     let s: f64;
     let z: f64;
@@ -268,7 +270,9 @@ pub fn erf(x: f64) -> f64 {
 /// Is `1 - erf(x)`. Is computed directly, so that you can use it to avoid
 /// the loss of precision that would result from subtracting
 /// large probabilities (on large `x`) from 1.
-pub fn erfc(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn erfc(x: f64) -> f64 {
     let r: f64;
     let s: f64;
     let z: f64;
diff --git a/src/math/erff.rs b/crates/libm/src/math/erff.rs
similarity index 97%
rename from src/math/erff.rs
rename to crates/libm/src/math/erff.rs
index 384052293..cfb53b4e0 100644
--- a/src/math/erff.rs
+++ b/crates/libm/src/math/erff.rs
@@ -130,7 +130,9 @@ fn erfc2(mut ix: u32, mut x: f32) -> f32 {
 /// Calculates an approximation to the “error function”, which estimates
 /// the probability that an observation will fall within x standard
 /// deviations of the mean (assuming a normal distribution).
-pub fn erff(x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn erff(x: f32) -> f32 {
     let r: f32;
     let s: f32;
     let z: f32;
@@ -179,7 +181,9 @@ pub fn erff(x: f32) -> f32 {
 /// Is `1 - erf(x)`. Is computed directly, so that you can use it to avoid
 /// the loss of precision that would result from subtracting
 /// large probabilities (on large `x`) from 1.
-pub fn erfcf(x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn erfcf(x: f32) -> f32 {
     let r: f32;
     let s: f32;
     let z: f32;
diff --git a/src/math/exp.rs b/crates/libm/src/math/exp.rs
similarity index 99%
rename from src/math/exp.rs
rename to crates/libm/src/math/exp.rs
index 5465b5693..633d4bda6 100644
--- a/src/math/exp.rs
+++ b/crates/libm/src/math/exp.rs
@@ -83,7 +83,7 @@ const P5: f64 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 /// (where *e* is the base of the natural system of logarithms, approximately 2.71828).
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn exp(mut x: f64) -> f64 {
+pub extern "C" fn exp(mut x: f64) -> f64 {
     let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
     let x1p_149 = f64::from_bits(0x36a0000000000000); // 0x1p-149 === 2 ^ -149
 
diff --git a/src/math/exp10.rs b/crates/libm/src/math/exp10.rs
similarity index 85%
rename from src/math/exp10.rs
rename to crates/libm/src/math/exp10.rs
index 9537f76f1..99f65bf6f 100644
--- a/src/math/exp10.rs
+++ b/crates/libm/src/math/exp10.rs
@@ -6,7 +6,9 @@ const P10: &[f64] = &[
     1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15,
 ];
 
-pub fn exp10(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn exp10(x: f64) -> f64 {
     let (mut y, n) = modf(x);
     let u: u64 = n.to_bits();
     /* fabs(n) < 16 without raising invalid on nan */
diff --git a/src/math/exp10f.rs b/crates/libm/src/math/exp10f.rs
similarity index 85%
rename from src/math/exp10f.rs
rename to crates/libm/src/math/exp10f.rs
index d45fff36e..f73cc2e22 100644
--- a/src/math/exp10f.rs
+++ b/crates/libm/src/math/exp10f.rs
@@ -6,7 +6,9 @@ const P10: &[f32] = &[
     1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7,
 ];
 
-pub fn exp10f(x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn exp10f(x: f32) -> f32 {
     let (mut y, n) = modff(x);
     let u = n.to_bits();
     /* fabsf(n) < 8 without raising invalid on nan */
diff --git a/src/math/exp2.rs b/crates/libm/src/math/exp2.rs
similarity index 98%
rename from src/math/exp2.rs
rename to crates/libm/src/math/exp2.rs
index c2192fde5..7c7436dac 100644
--- a/src/math/exp2.rs
+++ b/crates/libm/src/math/exp2.rs
@@ -28,7 +28,7 @@ use super::scalbn;
 
 const TBLSIZE: usize = 256;
 
-#[cfg_attr(rustfmt, rustfmt_skip)]
+#[rustfmt::skip]
 static TBL: [u64; TBLSIZE * 2] = [
     //  exp2(z + eps)          eps
     0x3fe6a09e667f3d5d, 0x3d39880000000000,
@@ -324,7 +324,7 @@ static TBL: [u64; TBLSIZE * 2] = [
 /// Calculate `2^x`, that is, 2 raised to the power `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn exp2(mut x: f64) -> f64 {
+pub extern "C" fn exp2(mut x: f64) -> f64 {
     let redux = f64::from_bits(0x4338000000000000) / TBLSIZE as f64;
     let p1 = f64::from_bits(0x3fe62e42fefa39ef);
     let p2 = f64::from_bits(0x3fcebfbdff82c575);
@@ -387,9 +387,3 @@ pub fn exp2(mut x: f64) -> f64 {
 
     scalbn(r, ki)
 }
-
-#[test]
-fn i0_wrap_test() {
-    let x = -3.0 / 256.0;
-    assert_eq!(exp2(x), f64::from_bits(0x3fefbdba3692d514));
-}
diff --git a/src/math/exp2f.rs b/crates/libm/src/math/exp2f.rs
similarity index 99%
rename from src/math/exp2f.rs
rename to crates/libm/src/math/exp2f.rs
index 12c9e76a4..2249bb6a9 100644
--- a/src/math/exp2f.rs
+++ b/crates/libm/src/math/exp2f.rs
@@ -75,7 +75,7 @@ static EXP2FT: [u64; TBLSIZE] = [
 /// Calculate `2^x`, that is, 2 raised to the power `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn exp2f(mut x: f32) -> f32 {
+pub extern "C" fn exp2f(mut x: f32) -> f32 {
     let redux = f32::from_bits(0x4b400000) / TBLSIZE as f32;
     let p1 = f32::from_bits(0x3f317218);
     let p2 = f32::from_bits(0x3e75fdf0);
diff --git a/src/math/expf.rs b/crates/libm/src/math/expf.rs
similarity index 98%
rename from src/math/expf.rs
rename to crates/libm/src/math/expf.rs
index 09323ec8d..b9d185be7 100644
--- a/src/math/expf.rs
+++ b/crates/libm/src/math/expf.rs
@@ -32,7 +32,7 @@ const P2: f32 = -2.7667332906e-3; /* -0xb55215.0p-32 */
 /// (where *e* is the base of the natural system of logarithms, approximately 2.71828).
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn expf(mut x: f32) -> f32 {
+pub extern "C" fn expf(mut x: f32) -> f32 {
     let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
     let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126  /*original 0x1p-149f    ??????????? */
     let mut hx = x.to_bits();
diff --git a/src/math/expm1.rs b/crates/libm/src/math/expm1.rs
similarity index 96%
rename from src/math/expm1.rs
rename to crates/libm/src/math/expm1.rs
index 0d43b4e10..9e7dd0e15 100644
--- a/src/math/expm1.rs
+++ b/crates/libm/src/math/expm1.rs
@@ -32,7 +32,7 @@ const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
 /// where using `exp(x)-1` would lose many significant digits.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn expm1(mut x: f64) -> f64 {
+pub extern "C" fn expm1(mut x: f64) -> f64 {
     let hi: f64;
     let lo: f64;
     let k: i32;
@@ -135,11 +135,3 @@ pub fn expm1(mut x: f64) -> f64 {
     }
     y
 }
-
-#[cfg(test)]
-mod tests {
-    #[test]
-    fn sanity_check() {
-        assert_eq!(super::expm1(1.1), 2.0041660239464334);
-    }
-}
diff --git a/src/math/expm1f.rs b/crates/libm/src/math/expm1f.rs
similarity index 98%
rename from src/math/expm1f.rs
rename to crates/libm/src/math/expm1f.rs
index 9bb223448..6117b6a06 100644
--- a/src/math/expm1f.rs
+++ b/crates/libm/src/math/expm1f.rs
@@ -34,7 +34,7 @@ const Q2: f32 = 1.5807170421e-3; /*  0xcf3010.0p-33 */
 /// where using `exp(x)-1` would lose many significant digits.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn expm1f(mut x: f32) -> f32 {
+pub extern "C" fn expm1f(mut x: f32) -> f32 {
     let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
 
     let mut hx = x.to_bits();
diff --git a/src/math/expo2.rs b/crates/libm/src/math/expo2.rs
similarity index 100%
rename from src/math/expo2.rs
rename to crates/libm/src/math/expo2.rs
diff --git a/src/math/fabs.rs b/crates/libm/src/math/fabs.rs
similarity index 84%
rename from src/math/fabs.rs
rename to crates/libm/src/math/fabs.rs
index 52a9adcbf..e7bc34dd9 100644
--- a/src/math/fabs.rs
+++ b/crates/libm/src/math/fabs.rs
@@ -5,12 +5,12 @@ use core::u64;
 /// by direct manipulation of the bit representation of `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fabs(x: f64) -> f64 {
+pub extern "C" fn fabs(x: f64) -> f64 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f64.abs` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::fabsf64(x) }
         }
     }
diff --git a/src/math/fabsf.rs b/crates/libm/src/math/fabsf.rs
similarity index 83%
rename from src/math/fabsf.rs
rename to crates/libm/src/math/fabsf.rs
index 5942d983a..35b109d46 100644
--- a/src/math/fabsf.rs
+++ b/crates/libm/src/math/fabsf.rs
@@ -3,12 +3,12 @@
 /// by direct manipulation of the bit representation of `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fabsf(x: f32) -> f32 {
+pub extern "C" fn fabsf(x: f32) -> f32 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f32.abs` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::fabsf32(x) }
         }
     }
diff --git a/src/math/fdim.rs b/crates/libm/src/math/fdim.rs
similarity index 90%
rename from src/math/fdim.rs
rename to crates/libm/src/math/fdim.rs
index 06edc9960..555fb93c4 100644
--- a/src/math/fdim.rs
+++ b/crates/libm/src/math/fdim.rs
@@ -10,7 +10,7 @@ use core::f64;
 /// A range error may occur.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fdim(x: f64, y: f64) -> f64 {
+pub extern "C" fn fdim(x: f64, y: f64) -> f64 {
     if x.is_nan() {
         x
     } else if y.is_nan() {
diff --git a/src/math/fdimf.rs b/crates/libm/src/math/fdimf.rs
similarity index 90%
rename from src/math/fdimf.rs
rename to crates/libm/src/math/fdimf.rs
index f1ad5896b..1e001cbab 100644
--- a/src/math/fdimf.rs
+++ b/crates/libm/src/math/fdimf.rs
@@ -10,7 +10,7 @@ use core::f32;
 /// A range error may occur.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fdimf(x: f32, y: f32) -> f32 {
+pub extern "C" fn fdimf(x: f32, y: f32) -> f32 {
     if x.is_nan() {
         x
     } else if y.is_nan() {
diff --git a/crates/libm/src/math/fenv.rs b/crates/libm/src/math/fenv.rs
new file mode 100644
index 000000000..652e60324
--- /dev/null
+++ b/crates/libm/src/math/fenv.rs
@@ -0,0 +1,33 @@
+// src: musl/src/fenv/fenv.c
+/* Dummy functions for archs lacking fenv implementation */
+
+pub(crate) const FE_UNDERFLOW: i32 = 0;
+pub(crate) const FE_INEXACT: i32 = 0;
+
+pub(crate) const FE_TONEAREST: i32 = 0;
+pub(crate) const FE_TOWARDZERO: i32 = 0;
+
+#[inline]
+pub(crate) fn feclearexcept(_mask: i32) -> i32 {
+    0
+}
+
+#[inline]
+pub(crate) fn feraiseexcept(_mask: i32) -> i32 {
+    0
+}
+
+#[inline]
+pub(crate) fn fetestexcept(_mask: i32) -> i32 {
+    0
+}
+
+#[inline]
+pub(crate) fn fegetround() -> i32 {
+    FE_TONEAREST
+}
+
+#[inline]
+pub(crate) fn fesetround(_r: i32) -> i32 {
+    0
+}
diff --git a/src/math/floor.rs b/crates/libm/src/math/floor.rs
similarity index 90%
rename from src/math/floor.rs
rename to crates/libm/src/math/floor.rs
index f6068c697..96125cc2a 100644
--- a/src/math/floor.rs
+++ b/crates/libm/src/math/floor.rs
@@ -7,12 +7,12 @@ const TOINT: f64 = 1. / f64::EPSILON;
 /// Finds the nearest integer less than or equal to `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn floor(x: f64) -> f64 {
+pub extern "C" fn floor(x: f64) -> f64 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f64.floor` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::floorf64(x) }
         }
     }
diff --git a/src/math/floorf.rs b/crates/libm/src/math/floorf.rs
similarity index 84%
rename from src/math/floorf.rs
rename to crates/libm/src/math/floorf.rs
index c04f18aee..b22c760df 100644
--- a/src/math/floorf.rs
+++ b/crates/libm/src/math/floorf.rs
@@ -5,12 +5,12 @@ use core::f32;
 /// Finds the nearest integer less than or equal to `x`.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn floorf(x: f32) -> f32 {
+pub extern "C" fn floorf(x: f32) -> f32 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f32.floor` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::floorf32(x) }
         }
     }
@@ -40,11 +40,3 @@ pub fn floorf(x: f32) -> f32 {
     }
     f32::from_bits(ui)
 }
-
-#[cfg(test)]
-mod tests {
-    #[test]
-    fn no_overflow() {
-        assert_eq!(super::floorf(0.5), 0.0);
-    }
-}
diff --git a/src/math/fma.rs b/crates/libm/src/math/fma.rs
similarity index 97%
rename from src/math/fma.rs
rename to crates/libm/src/math/fma.rs
index 07d90f8b7..01d4faaa6 100644
--- a/src/math/fma.rs
+++ b/crates/libm/src/math/fma.rs
@@ -55,7 +55,7 @@ fn mul(x: u64, y: u64) -> (u64, u64) {
 /// according to the rounding mode characterized by the value of FLT_ROUNDS.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fma(x: f64, y: f64, z: f64) -> f64 {
+pub extern "C" fn fma(x: f64, y: f64, z: f64) -> f64 {
     let x1p63: f64 = f64::from_bits(0x43e0000000000000); // 0x1p63 === 2 ^ 63
     let x0_ffffff8p_63 = f64::from_bits(0x3bfffffff0000000); // 0x0.ffffff8p-63
 
@@ -126,8 +126,8 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
     } else {
         /* r -= z */
         let t = rlo;
-        rlo -= zlo;
-        rhi = rhi - zhi - (t < rlo) as u64;
+        rlo = rlo.wrapping_sub(zlo);
+        rhi = rhi.wrapping_sub(zhi.wrapping_sub((t < rlo) as u64));
         if (rhi >> 63) != 0 {
             rlo = (-(rlo as i64)) as u64;
             rhi = (-(rhi as i64)) as u64 - (rlo != 0) as u64;
diff --git a/src/math/fmaf.rs b/crates/libm/src/math/fmaf.rs
similarity index 98%
rename from src/math/fmaf.rs
rename to crates/libm/src/math/fmaf.rs
index e77e0fa4a..a6609d918 100644
--- a/src/math/fmaf.rs
+++ b/crates/libm/src/math/fmaf.rs
@@ -48,7 +48,7 @@ use super::fenv::{
 /// according to the rounding mode characterized by the value of FLT_ROUNDS.
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 {
+pub extern "C" fn fmaf(x: f32, y: f32, mut z: f32) -> f32 {
     let xy: f64;
     let mut result: f64;
     let mut ui: u64;
diff --git a/src/math/fmax.rs b/crates/libm/src/math/fmax.rs
similarity index 93%
rename from src/math/fmax.rs
rename to crates/libm/src/math/fmax.rs
index 22016d11c..9a7de6aa1 100644
--- a/src/math/fmax.rs
+++ b/crates/libm/src/math/fmax.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fmax(x: f64, y: f64) -> f64 {
+pub extern "C" fn fmax(x: f64, y: f64) -> f64 {
     // IEEE754 says: maxNum(x, y) is the canonicalized number y if x < y, x if y < x, the
     // canonicalized number if one operand is a number and the other a quiet NaN. Otherwise it
     // is either x or y, canonicalized (this means results might differ among implementations).
diff --git a/src/math/fmaxf.rs b/crates/libm/src/math/fmaxf.rs
similarity index 93%
rename from src/math/fmaxf.rs
rename to crates/libm/src/math/fmaxf.rs
index a883fdaef..84aded4e5 100644
--- a/src/math/fmaxf.rs
+++ b/crates/libm/src/math/fmaxf.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fmaxf(x: f32, y: f32) -> f32 {
+pub extern "C" fn fmaxf(x: f32, y: f32) -> f32 {
     // IEEE754 says: maxNum(x, y) is the canonicalized number y if x < y, x if y < x, the
     // canonicalized number if one operand is a number and the other a quiet NaN. Otherwise it
     // is either x or y, canonicalized (this means results might differ among implementations).
diff --git a/src/math/fmin.rs b/crates/libm/src/math/fmin.rs
similarity index 93%
rename from src/math/fmin.rs
rename to crates/libm/src/math/fmin.rs
index d1ccc3a46..537143bc2 100644
--- a/src/math/fmin.rs
+++ b/crates/libm/src/math/fmin.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fmin(x: f64, y: f64) -> f64 {
+pub extern "C" fn fmin(x: f64, y: f64) -> f64 {
     // IEEE754 says: minNum(x, y) is the canonicalized number x if x < y, y if y < x, the
     // canonicalized number if one operand is a number and the other a quiet NaN. Otherwise it
     // is either x or y, canonicalized (this means results might differ among implementations).
diff --git a/src/math/fminf.rs b/crates/libm/src/math/fminf.rs
similarity index 93%
rename from src/math/fminf.rs
rename to crates/libm/src/math/fminf.rs
index 43ec97cb5..1cef90e5b 100644
--- a/src/math/fminf.rs
+++ b/crates/libm/src/math/fminf.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fminf(x: f32, y: f32) -> f32 {
+pub extern "C" fn fminf(x: f32, y: f32) -> f32 {
     // IEEE754 says: minNum(x, y) is the canonicalized number x if x < y, y if y < x, the
     // canonicalized number if one operand is a number and the other a quiet NaN. Otherwise it
     // is either x or y, canonicalized (this means results might differ among implementations).
diff --git a/src/math/fmod.rs b/crates/libm/src/math/fmod.rs
similarity index 97%
rename from src/math/fmod.rs
rename to crates/libm/src/math/fmod.rs
index 2cdd8a9ba..cb0c6cf5a 100644
--- a/src/math/fmod.rs
+++ b/crates/libm/src/math/fmod.rs
@@ -2,7 +2,7 @@ use core::u64;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fmod(x: f64, y: f64) -> f64 {
+pub extern "C" fn fmod(x: f64, y: f64) -> f64 {
     let mut uxi = x.to_bits();
     let mut uyi = y.to_bits();
     let mut ex = (uxi >> 52 & 0x7ff) as i64;
diff --git a/src/math/fmodf.rs b/crates/libm/src/math/fmodf.rs
similarity index 97%
rename from src/math/fmodf.rs
rename to crates/libm/src/math/fmodf.rs
index 3e6779a93..25cabe1d4 100644
--- a/src/math/fmodf.rs
+++ b/crates/libm/src/math/fmodf.rs
@@ -3,7 +3,7 @@ use core::u32;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn fmodf(x: f32, y: f32) -> f32 {
+pub extern "C" fn fmodf(x: f32, y: f32) -> f32 {
     let mut uxi = x.to_bits();
     let mut uyi = y.to_bits();
     let mut ex = (uxi >> 23 & 0xff) as i32;
diff --git a/src/math/frexp.rs b/crates/libm/src/math/frexp.rs
similarity index 79%
rename from src/math/frexp.rs
rename to crates/libm/src/math/frexp.rs
index badad786a..676efa28b 100644
--- a/src/math/frexp.rs
+++ b/crates/libm/src/math/frexp.rs
@@ -1,4 +1,6 @@
-pub fn frexp(x: f64) -> (f64, i32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn frexp(x: f64) -> (f64, i32) {
     let mut y = x.to_bits();
     let ee = ((y >> 52) & 0x7ff) as i32;
 
diff --git a/src/math/frexpf.rs b/crates/libm/src/math/frexpf.rs
similarity index 78%
rename from src/math/frexpf.rs
rename to crates/libm/src/math/frexpf.rs
index 2919c0ab0..2d63f6646 100644
--- a/src/math/frexpf.rs
+++ b/crates/libm/src/math/frexpf.rs
@@ -1,4 +1,6 @@
-pub fn frexpf(x: f32) -> (f32, i32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn frexpf(x: f32) -> (f32, i32) {
     let mut y = x.to_bits();
     let ee: i32 = ((y >> 23) & 0xff) as i32;
 
diff --git a/src/math/hypot.rs b/crates/libm/src/math/hypot.rs
similarity index 96%
rename from src/math/hypot.rs
rename to crates/libm/src/math/hypot.rs
index e53baf539..d77eac945 100644
--- a/src/math/hypot.rs
+++ b/crates/libm/src/math/hypot.rs
@@ -20,7 +20,7 @@ fn sq(x: f64) -> (f64, f64) {
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn hypot(mut x: f64, mut y: f64) -> f64 {
+pub extern "C" fn hypot(mut x: f64, mut y: f64) -> f64 {
     let x1p700 = f64::from_bits(0x6bb0000000000000); // 0x1p700 === 2 ^ 700
     let x1p_700 = f64::from_bits(0x1430000000000000); // 0x1p-700 === 2 ^ -700
 
diff --git a/src/math/hypotf.rs b/crates/libm/src/math/hypotf.rs
similarity index 94%
rename from src/math/hypotf.rs
rename to crates/libm/src/math/hypotf.rs
index 4636b8f1d..a15735cac 100644
--- a/src/math/hypotf.rs
+++ b/crates/libm/src/math/hypotf.rs
@@ -4,7 +4,7 @@ use super::sqrtf;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn hypotf(mut x: f32, mut y: f32) -> f32 {
+pub extern "C" fn hypotf(mut x: f32, mut y: f32) -> f32 {
     let x1p90 = f32::from_bits(0x6c800000); // 0x1p90f === 2 ^ 90
     let x1p_90 = f32::from_bits(0x12800000); // 0x1p-90f === 2 ^ -90
 
diff --git a/src/math/ilogb.rs b/crates/libm/src/math/ilogb.rs
similarity index 85%
rename from src/math/ilogb.rs
rename to crates/libm/src/math/ilogb.rs
index 0a380b7ef..3c2fed309 100644
--- a/src/math/ilogb.rs
+++ b/crates/libm/src/math/ilogb.rs
@@ -1,7 +1,9 @@
 const FP_ILOGBNAN: i32 = -1 - 0x7fffffff;
 const FP_ILOGB0: i32 = FP_ILOGBNAN;
 
-pub fn ilogb(x: f64) -> i32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn ilogb(x: f64) -> i32 {
     let mut i: u64 = x.to_bits();
     let e = ((i >> 52) & 0x7ff) as i32;
 
diff --git a/src/math/ilogbf.rs b/crates/libm/src/math/ilogbf.rs
similarity index 84%
rename from src/math/ilogbf.rs
rename to crates/libm/src/math/ilogbf.rs
index b384fa4b2..09230007c 100644
--- a/src/math/ilogbf.rs
+++ b/crates/libm/src/math/ilogbf.rs
@@ -1,7 +1,9 @@
 const FP_ILOGBNAN: i32 = -1 - 0x7fffffff;
 const FP_ILOGB0: i32 = FP_ILOGBNAN;
 
-pub fn ilogbf(x: f32) -> i32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn ilogbf(x: f32) -> i32 {
     let mut i = x.to_bits();
     let e = ((i >> 23) & 0xff) as i32;
 
diff --git a/src/math/j0.rs b/crates/libm/src/math/j0.rs
similarity index 98%
rename from src/math/j0.rs
rename to crates/libm/src/math/j0.rs
index c4258ccca..516d0048c 100644
--- a/src/math/j0.rs
+++ b/crates/libm/src/math/j0.rs
@@ -109,7 +109,9 @@ const S02: f64 = 1.16926784663337450260e-04; /* 0x3F1EA6D2, 0xDD57DBF4 */
 const S03: f64 = 5.13546550207318111446e-07; /* 0x3EA13B54, 0xCE84D5A9 */
 const S04: f64 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
 
-pub fn j0(mut x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn j0(mut x: f64) -> f64 {
     let z: f64;
     let r: f64;
     let s: f64;
@@ -162,7 +164,9 @@ const V02: f64 = 7.60068627350353253702e-05; /* 0x3F13ECBB, 0xF578C6C1 */
 const V03: f64 = 2.59150851840457805467e-07; /* 0x3E91642D, 0x7FF202FD */
 const V04: f64 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
 
-pub fn y0(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn y0(x: f64) -> f64 {
     let z: f64;
     let u: f64;
     let v: f64;
diff --git a/src/math/j0f.rs b/crates/libm/src/math/j0f.rs
similarity index 97%
rename from src/math/j0f.rs
rename to crates/libm/src/math/j0f.rs
index 91c03dbbc..a93fbf9b8 100644
--- a/src/math/j0f.rs
+++ b/crates/libm/src/math/j0f.rs
@@ -62,7 +62,9 @@ const S02: f32 = 1.1692678527e-04; /* 0x38f53697 */
 const S03: f32 = 5.1354652442e-07; /* 0x3509daa6 */
 const S04: f32 = 1.1661400734e-09; /* 0x30a045e8 */
 
-pub fn j0f(mut x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn j0f(mut x: f32) -> f32 {
     let z: f32;
     let r: f32;
     let s: f32;
@@ -106,8 +108,9 @@ const V01: f32 = 1.2730483897e-02; /* 0x3c509385 */
 const V02: f32 = 7.6006865129e-05; /* 0x389f65e0 */
 const V03: f32 = 2.5915085189e-07; /* 0x348b216c */
 const V04: f32 = 4.4111031494e-10; /* 0x2ff280c2 */
-
-pub fn y0f(x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn y0f(x: f32) -> f32 {
     let z: f32;
     let u: f32;
     let v: f32;
diff --git a/src/math/j1.rs b/crates/libm/src/math/j1.rs
similarity index 98%
rename from src/math/j1.rs
rename to crates/libm/src/math/j1.rs
index 02a65ca5a..23a01e3ad 100644
--- a/src/math/j1.rs
+++ b/crates/libm/src/math/j1.rs
@@ -113,7 +113,9 @@ const S03: f64 = 1.17718464042623683263e-06; /* 0x3EB3BFF8, 0x333F8498 */
 const S04: f64 = 5.04636257076217042715e-09; /* 0x3E35AC88, 0xC97DFF2C */
 const S05: f64 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
 
-pub fn j1(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn j1(x: f64) -> f64 {
     let mut z: f64;
     let r: f64;
     let s: f64;
@@ -158,7 +160,9 @@ const V0: [f64; 5] = [
     1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
 ];
 
-pub fn y1(x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn y1(x: f64) -> f64 {
     let z: f64;
     let u: f64;
     let v: f64;
diff --git a/src/math/j1f.rs b/crates/libm/src/math/j1f.rs
similarity index 96%
rename from src/math/j1f.rs
rename to crates/libm/src/math/j1f.rs
index 5095894d7..377c1cd3d 100644
--- a/src/math/j1f.rs
+++ b/crates/libm/src/math/j1f.rs
@@ -63,7 +63,9 @@ const S03: f32 = 1.1771846857e-06; /* 0x359dffc2 */
 const S04: f32 = 5.0463624390e-09; /* 0x31ad6446 */
 const S05: f32 = 1.2354227016e-11; /* 0x2d59567e */
 
-pub fn j1f(x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn j1f(x: f32) -> f32 {
     let mut z: f32;
     let r: f32;
     let s: f32;
@@ -107,7 +109,9 @@ const V0: [f32; 5] = [
     1.6655924903e-11, /* 0x2d9281cf */
 ];
 
-pub fn y1f(x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn y1f(x: f32) -> f32 {
     let z: f32;
     let u: f32;
     let v: f32;
@@ -356,17 +360,3 @@ fn qonef(x: f32) -> f32 {
     s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
     return (0.375 + r / s) / x;
 }
-
-#[cfg(test)]
-mod tests {
-    use super::{j1f, y1f};
-    #[test]
-    fn test_j1f_2488() {
-        // 0x401F3E49
-        assert_eq!(j1f(2.4881766_f32), 0.49999475_f32);
-    }
-    #[test]
-    fn test_y1f_2002() {
-        assert_eq!(y1f(2.0000002_f32), -0.10703229_f32);
-    }
-}
diff --git a/src/math/jn.rs b/crates/libm/src/math/jn.rs
similarity index 97%
rename from src/math/jn.rs
rename to crates/libm/src/math/jn.rs
index 1be167f84..682b868e9 100644
--- a/src/math/jn.rs
+++ b/crates/libm/src/math/jn.rs
@@ -38,7 +38,9 @@ use super::{cos, fabs, get_high_word, get_low_word, j0, j1, log, sin, sqrt, y0,
 
 const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
 
-pub fn jn(n: i32, mut x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn jn(n: i32, mut x: f64) -> f64 {
     let mut ix: u32;
     let lx: u32;
     let nm1: i32;
@@ -251,7 +253,9 @@ pub fn jn(n: i32, mut x: f64) -> f64 {
     }
 }
 
-pub fn yn(n: i32, x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn yn(n: i32, x: f64) -> f64 {
     let mut ix: u32;
     let lx: u32;
     let mut ib: u32;
diff --git a/src/math/jnf.rs b/crates/libm/src/math/jnf.rs
similarity index 96%
rename from src/math/jnf.rs
rename to crates/libm/src/math/jnf.rs
index 360f62e20..efb4a1df9 100644
--- a/src/math/jnf.rs
+++ b/crates/libm/src/math/jnf.rs
@@ -15,7 +15,9 @@
 
 use super::{fabsf, j0f, j1f, logf, y0f, y1f};
 
-pub fn jnf(n: i32, mut x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn jnf(n: i32, mut x: f32) -> f32 {
     let mut ix: u32;
     let mut nm1: i32;
     let mut sign: bool;
@@ -195,7 +197,9 @@ pub fn jnf(n: i32, mut x: f32) -> f32 {
     }
 }
 
-pub fn ynf(n: i32, x: f32) -> f32 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn ynf(n: i32, x: f32) -> f32 {
     let mut ix: u32;
     let mut ib: u32;
     let nm1: i32;
diff --git a/src/math/k_cos.rs b/crates/libm/src/math/k_cos.rs
similarity index 100%
rename from src/math/k_cos.rs
rename to crates/libm/src/math/k_cos.rs
diff --git a/src/math/k_cosf.rs b/crates/libm/src/math/k_cosf.rs
similarity index 100%
rename from src/math/k_cosf.rs
rename to crates/libm/src/math/k_cosf.rs
diff --git a/src/math/k_expo2.rs b/crates/libm/src/math/k_expo2.rs
similarity index 100%
rename from src/math/k_expo2.rs
rename to crates/libm/src/math/k_expo2.rs
diff --git a/src/math/k_expo2f.rs b/crates/libm/src/math/k_expo2f.rs
similarity index 100%
rename from src/math/k_expo2f.rs
rename to crates/libm/src/math/k_expo2f.rs
diff --git a/src/math/k_sin.rs b/crates/libm/src/math/k_sin.rs
similarity index 100%
rename from src/math/k_sin.rs
rename to crates/libm/src/math/k_sin.rs
diff --git a/src/math/k_sinf.rs b/crates/libm/src/math/k_sinf.rs
similarity index 100%
rename from src/math/k_sinf.rs
rename to crates/libm/src/math/k_sinf.rs
diff --git a/src/math/k_tan.rs b/crates/libm/src/math/k_tan.rs
similarity index 100%
rename from src/math/k_tan.rs
rename to crates/libm/src/math/k_tan.rs
diff --git a/src/math/k_tanf.rs b/crates/libm/src/math/k_tanf.rs
similarity index 100%
rename from src/math/k_tanf.rs
rename to crates/libm/src/math/k_tanf.rs
diff --git a/src/math/ldexp.rs b/crates/libm/src/math/ldexp.rs
similarity index 66%
rename from src/math/ldexp.rs
rename to crates/libm/src/math/ldexp.rs
index 780ddfc11..70dabbc54 100644
--- a/src/math/ldexp.rs
+++ b/crates/libm/src/math/ldexp.rs
@@ -1,5 +1,5 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn ldexp(x: f64, n: i32) -> f64 {
+pub extern "C" fn ldexp(x: f64, n: i32) -> f64 {
     super::scalbn(x, n)
 }
diff --git a/src/math/ldexpf.rs b/crates/libm/src/math/ldexpf.rs
similarity index 65%
rename from src/math/ldexpf.rs
rename to crates/libm/src/math/ldexpf.rs
index 70935a002..bb4902670 100644
--- a/src/math/ldexpf.rs
+++ b/crates/libm/src/math/ldexpf.rs
@@ -1,5 +1,5 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn ldexpf(x: f32, n: i32) -> f32 {
+pub extern "C" fn ldexpf(x: f32, n: i32) -> f32 {
     super::scalbnf(x, n)
 }
diff --git a/crates/libm/src/math/lgamma.rs b/crates/libm/src/math/lgamma.rs
new file mode 100644
index 000000000..d1bfa8dfb
--- /dev/null
+++ b/crates/libm/src/math/lgamma.rs
@@ -0,0 +1,7 @@
+use super::lgamma_r;
+
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn lgamma(x: f64) -> f64 {
+    lgamma_r(x).0
+}
diff --git a/src/math/lgamma_r.rs b/crates/libm/src/math/lgamma_r.rs
similarity index 99%
rename from src/math/lgamma_r.rs
rename to crates/libm/src/math/lgamma_r.rs
index 382a501fc..8d9bfcc85 100644
--- a/src/math/lgamma_r.rs
+++ b/crates/libm/src/math/lgamma_r.rs
@@ -164,7 +164,9 @@ fn sin_pi(mut x: f64) -> f64 {
     }
 }
 
-pub fn lgamma_r(mut x: f64) -> (f64, i32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn lgamma_r(mut x: f64) -> (f64, i32) {
     let u: u64 = x.to_bits();
     let mut t: f64;
     let y: f64;
diff --git a/crates/libm/src/math/lgammaf.rs b/crates/libm/src/math/lgammaf.rs
new file mode 100644
index 000000000..2eabb59d6
--- /dev/null
+++ b/crates/libm/src/math/lgammaf.rs
@@ -0,0 +1,7 @@
+use super::lgammaf_r;
+
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn lgammaf(x: f32) -> f32 {
+    lgammaf_r(x).0
+}
diff --git a/src/math/lgammaf_r.rs b/crates/libm/src/math/lgammaf_r.rs
similarity index 98%
rename from src/math/lgammaf_r.rs
rename to crates/libm/src/math/lgammaf_r.rs
index 0745359a2..416cb54ed 100644
--- a/src/math/lgammaf_r.rs
+++ b/crates/libm/src/math/lgammaf_r.rs
@@ -98,8 +98,9 @@ fn sin_pi(mut x: f32) -> f32 {
         0 | _ => k_sinf(y),
     }
 }
-
-pub fn lgammaf_r(mut x: f32) -> (f32, i32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn lgammaf_r(mut x: f32) -> (f32, i32) {
     let u = x.to_bits();
     let mut t: f32;
     let y: f32;
diff --git a/src/math/log.rs b/crates/libm/src/math/log.rs
similarity index 99%
rename from src/math/log.rs
rename to crates/libm/src/math/log.rs
index 4126e413b..51cf810a0 100644
--- a/src/math/log.rs
+++ b/crates/libm/src/math/log.rs
@@ -72,7 +72,7 @@ const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn log(mut x: f64) -> f64 {
+pub extern "C" fn log(mut x: f64) -> f64 {
     let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
 
     let mut ui = x.to_bits();
diff --git a/src/math/log10.rs b/crates/libm/src/math/log10.rs
similarity index 98%
rename from src/math/log10.rs
rename to crates/libm/src/math/log10.rs
index c99696040..e1caec64b 100644
--- a/src/math/log10.rs
+++ b/crates/libm/src/math/log10.rs
@@ -33,7 +33,7 @@ const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn log10(mut x: f64) -> f64 {
+pub extern "C" fn log10(mut x: f64) -> f64 {
     let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
 
     let mut ui: u64 = x.to_bits();
diff --git a/src/math/log10f.rs b/crates/libm/src/math/log10f.rs
similarity index 98%
rename from src/math/log10f.rs
rename to crates/libm/src/math/log10f.rs
index 9845cda5d..e97ecf2a5 100644
--- a/src/math/log10f.rs
+++ b/crates/libm/src/math/log10f.rs
@@ -27,7 +27,7 @@ const LG4: f32 = 0.24279078841; /* 0xf89e26.0p-26 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn log10f(mut x: f32) -> f32 {
+pub extern "C" fn log10f(mut x: f32) -> f32 {
     let x1p25f = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
 
     let mut ui: u32 = x.to_bits();
diff --git a/src/math/log1p.rs b/crates/libm/src/math/log1p.rs
similarity index 99%
rename from src/math/log1p.rs
rename to crates/libm/src/math/log1p.rs
index cd7045ac9..e16e11f1c 100644
--- a/src/math/log1p.rs
+++ b/crates/libm/src/math/log1p.rs
@@ -67,7 +67,7 @@ const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn log1p(x: f64) -> f64 {
+pub extern "C" fn log1p(x: f64) -> f64 {
     let mut ui: u64 = x.to_bits();
     let hfsq: f64;
     let mut f: f64 = 0.;
diff --git a/src/math/log1pf.rs b/crates/libm/src/math/log1pf.rs
similarity index 98%
rename from src/math/log1pf.rs
rename to crates/libm/src/math/log1pf.rs
index 8e9651357..717295b90 100644
--- a/src/math/log1pf.rs
+++ b/crates/libm/src/math/log1pf.rs
@@ -22,7 +22,7 @@ const LG4: f32 = 0.24279078841; /* 0xf89e26.0p-26 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn log1pf(x: f32) -> f32 {
+pub extern "C" fn log1pf(x: f32) -> f32 {
     let mut ui: u32 = x.to_bits();
     let hfsq: f32;
     let mut f: f32 = 0.;
diff --git a/src/math/log2.rs b/crates/libm/src/math/log2.rs
similarity index 98%
rename from src/math/log2.rs
rename to crates/libm/src/math/log2.rs
index a3d43e55c..8b304855c 100644
--- a/src/math/log2.rs
+++ b/crates/libm/src/math/log2.rs
@@ -31,7 +31,7 @@ const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn log2(mut x: f64) -> f64 {
+pub extern "C" fn log2(mut x: f64) -> f64 {
     let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
 
     let mut ui: u64 = x.to_bits();
diff --git a/src/math/log2f.rs b/crates/libm/src/math/log2f.rs
similarity index 98%
rename from src/math/log2f.rs
rename to crates/libm/src/math/log2f.rs
index 53a37e503..cf4ebed2a 100644
--- a/src/math/log2f.rs
+++ b/crates/libm/src/math/log2f.rs
@@ -25,7 +25,7 @@ const LG4: f32 = 0.24279078841; /* 0xf89e26.0p-26 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn log2f(mut x: f32) -> f32 {
+pub extern "C" fn log2f(mut x: f32) -> f32 {
     let x1p25f = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
 
     let mut ui: u32 = x.to_bits();
diff --git a/src/math/logf.rs b/crates/libm/src/math/logf.rs
similarity index 97%
rename from src/math/logf.rs
rename to crates/libm/src/math/logf.rs
index 95195601c..6a0a9f7a6 100644
--- a/src/math/logf.rs
+++ b/crates/libm/src/math/logf.rs
@@ -23,7 +23,7 @@ const LG4: f32 = 0.24279078841; /*  0xf89e26.0p-26 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn logf(mut x: f32) -> f32 {
+pub extern "C" fn logf(mut x: f32) -> f32 {
     let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
 
     let mut ix = x.to_bits();
diff --git a/src/math/mod.rs b/crates/libm/src/math/mod.rs
similarity index 86%
rename from src/math/mod.rs
rename to crates/libm/src/math/mod.rs
index 48b400a92..87424f8ce 100644
--- a/src/math/mod.rs
+++ b/crates/libm/src/math/mod.rs
@@ -1,5 +1,42 @@
+#![allow(
+    clippy::many_single_char_names,
+    clippy::similar_names,
+    clippy::unreadable_literal,
+    clippy::unseparated_literal_suffix,
+    clippy::doc_markdown,
+    clippy::needless_return,
+    clippy::if_not_else,
+    clippy::int_plus_one,
+    clippy::large_digit_groups,
+    clippy::items_after_statements,
+    clippy::collapsible_if,
+    clippy::mixed_case_hex_literals,
+    clippy::precedence,
+    clippy::suspicious_else_formatting,
+    clippy::cast_lossless,
+    clippy::cast_possible_truncation,
+    clippy::cast_precision_loss,
+    clippy::cast_sign_loss,
+    clippy::cast_possible_wrap,
+    clippy::excessive_precision,
+    clippy::range_minus_one,
+    clippy::cognitive_complexity,
+    clippy::eq_op,
+    clippy::shadow_unrelated,
+    clippy::needless_range_loop,
+    clippy::float_cmp,
+    clippy::assign_op_pattern,
+    clippy::zero_divided_by_zero,
+    clippy::misrefactored_assign_op,
+    clippy::if_same_then_else,
+    clippy::verbose_bit_mask,
+    clippy::replace_consts,
+    clippy::used_underscore_binding,
+)]
+
 macro_rules! force_eval {
     ($e:expr) => {
+        #[allow(unused_unsafe)]
         unsafe {
             ::core::ptr::read_volatile(&$e);
         }
@@ -146,6 +183,8 @@ mod log2f;
 mod logf;
 mod modf;
 mod modff;
+mod nextafter;
+mod nextafterf;
 mod pow;
 mod powf;
 mod remainder;
@@ -258,6 +297,8 @@ pub use self::log2f::log2f;
 pub use self::logf::logf;
 pub use self::modf::modf;
 pub use self::modff::modff;
+pub use self::nextafter::nextafter;
+pub use self::nextafterf::nextafterf;
 pub use self::pow::pow;
 pub use self::powf::powf;
 pub use self::remainder::remainder;
diff --git a/src/math/modf.rs b/crates/libm/src/math/modf.rs
similarity index 85%
rename from src/math/modf.rs
rename to crates/libm/src/math/modf.rs
index bcab33a81..e564c46a6 100644
--- a/src/math/modf.rs
+++ b/crates/libm/src/math/modf.rs
@@ -1,4 +1,6 @@
-pub fn modf(x: f64) -> (f64, f64) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn modf(x: f64) -> (f64, f64) {
     let rv2: f64;
     let mut u = x.to_bits();
     let mask: u64;
diff --git a/src/math/modff.rs b/crates/libm/src/math/modff.rs
similarity index 85%
rename from src/math/modff.rs
rename to crates/libm/src/math/modff.rs
index 56ece12e3..f79b87fe2 100644
--- a/src/math/modff.rs
+++ b/crates/libm/src/math/modff.rs
@@ -1,4 +1,6 @@
-pub fn modff(x: f32) -> (f32, f32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn modff(x: f32) -> (f32, f32) {
     let rv2: f32;
     let mut u: u32 = x.to_bits();
     let mask: u32;
diff --git a/crates/libm/src/math/nextafter.rs b/crates/libm/src/math/nextafter.rs
new file mode 100644
index 000000000..99c0db63a
--- /dev/null
+++ b/crates/libm/src/math/nextafter.rs
@@ -0,0 +1,38 @@
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn nextafter(x: f64, y: f64) -> f64 {
+    if x.is_nan() || y.is_nan() {
+        return x + y;
+    }
+
+    let mut ux_i = x.to_bits();
+    let uy_i = y.to_bits();
+    if ux_i == uy_i {
+        return y;
+    }
+
+    let ax = ux_i & !1_u64 / 2;
+    let ay = uy_i & !1_u64 / 2;
+    if ax == 0 {
+        if ay == 0 {
+            return y;
+        }
+        ux_i = (uy_i & 1_u64 << 63) | 1;
+    } else if ax > ay || ((ux_i ^ uy_i) & 1_u64 << 63) != 0 {
+        ux_i -= 1;
+    } else {
+        ux_i += 1;
+    }
+
+    let e = ux_i.wrapping_shr(52 & 0x7ff);
+    // raise overflow if ux.f is infinite and x is finite
+    if e == 0x7ff {
+        force_eval!(x + x);
+    }
+    let ux_f = f64::from_bits(ux_i);
+    // raise underflow if ux.f is subnormal or zero
+    if e == 0 {
+        force_eval!(x * x + ux_f * ux_f);
+    }
+    ux_f
+}
diff --git a/crates/libm/src/math/nextafterf.rs b/crates/libm/src/math/nextafterf.rs
new file mode 100644
index 000000000..9d0be54b1
--- /dev/null
+++ b/crates/libm/src/math/nextafterf.rs
@@ -0,0 +1,38 @@
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn nextafterf(x: f32, y: f32) -> f32 {
+    if x.is_nan() || y.is_nan() {
+        return x + y;
+    }
+
+    let mut ux_i = x.to_bits();
+    let uy_i = y.to_bits();
+    if ux_i == uy_i {
+        return y;
+    }
+
+    let ax = ux_i & 0x7fff_ffff_u32;
+    let ay = uy_i & 0x7fff_ffff_u32;
+    if ax == 0 {
+        if ay == 0 {
+            return y;
+        }
+        ux_i = (uy_i & 0x8000_0000_u32) | 1;
+    } else if ax > ay || ((ux_i ^ uy_i) & 0x8000_0000_u32) != 0 {
+        ux_i -= 1;
+    } else {
+        ux_i += 1;
+    }
+
+    let e = ux_i.wrapping_shr(0x7f80_0000_u32);
+    // raise overflow if ux_f is infinite and x is finite
+    if e == 0x7f80_0000_u32 {
+        force_eval!(x + x);
+    }
+    let ux_f = f32::from_bits(ux_i);
+    // raise underflow if ux_f is subnormal or zero
+    if e == 0 {
+        force_eval!(x * x + ux_f * ux_f);
+    }
+    ux_f
+}
diff --git a/src/math/pow.rs b/crates/libm/src/math/pow.rs
similarity index 63%
rename from src/math/pow.rs
rename to crates/libm/src/math/pow.rs
index 068a4ec47..362a62442 100644
--- a/src/math/pow.rs
+++ b/crates/libm/src/math/pow.rs
@@ -91,7 +91,7 @@ const IVLN2_L: f64 = 1.92596299112661746887e-08; /* 0x3e54ae0b_f85ddf44 =1/ln2 t
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn pow(x: f64, y: f64) -> f64 {
+pub extern "C" fn pow(x: f64, y: f64) -> f64 {
     let t1: f64;
     let t2: f64;
 
@@ -408,227 +408,3 @@ pub fn pow(x: f64, y: f64) -> f64 {
 
     s * z
 }
-
-#[cfg(test)]
-mod tests {
-    extern crate core;
-
-    use self::core::f64::consts::{E, PI};
-    use self::core::f64::{EPSILON, INFINITY, MAX, MIN, MIN_POSITIVE, NAN, NEG_INFINITY};
-    use super::pow;
-
-    const POS_ZERO: &[f64] = &[0.0];
-    const NEG_ZERO: &[f64] = &[-0.0];
-    const POS_ONE: &[f64] = &[1.0];
-    const NEG_ONE: &[f64] = &[-1.0];
-    const POS_FLOATS: &[f64] = &[99.0 / 70.0, E, PI];
-    const NEG_FLOATS: &[f64] = &[-99.0 / 70.0, -E, -PI];
-    const POS_SMALL_FLOATS: &[f64] = &[(1.0 / 2.0), MIN_POSITIVE, EPSILON];
-    const NEG_SMALL_FLOATS: &[f64] = &[-(1.0 / 2.0), -MIN_POSITIVE, -EPSILON];
-    const POS_EVENS: &[f64] = &[2.0, 6.0, 8.0, 10.0, 22.0, 100.0, MAX];
-    const NEG_EVENS: &[f64] = &[MIN, -100.0, -22.0, -10.0, -8.0, -6.0, -2.0];
-    const POS_ODDS: &[f64] = &[3.0, 7.0];
-    const NEG_ODDS: &[f64] = &[-7.0, -3.0];
-    const NANS: &[f64] = &[NAN];
-    const POS_INF: &[f64] = &[INFINITY];
-    const NEG_INF: &[f64] = &[NEG_INFINITY];
-
-    const ALL: &[&[f64]] = &[
-        POS_ZERO,
-        NEG_ZERO,
-        NANS,
-        NEG_SMALL_FLOATS,
-        POS_SMALL_FLOATS,
-        NEG_FLOATS,
-        POS_FLOATS,
-        NEG_EVENS,
-        POS_EVENS,
-        NEG_ODDS,
-        POS_ODDS,
-        NEG_INF,
-        POS_INF,
-        NEG_ONE,
-        POS_ONE,
-    ];
-    const POS: &[&[f64]] = &[POS_ZERO, POS_ODDS, POS_ONE, POS_FLOATS, POS_EVENS, POS_INF];
-    const NEG: &[&[f64]] = &[NEG_ZERO, NEG_ODDS, NEG_ONE, NEG_FLOATS, NEG_EVENS, NEG_INF];
-
-    fn pow_test(base: f64, exponent: f64, expected: f64) {
-        let res = pow(base, exponent);
-        assert!(
-            if expected.is_nan() {
-                res.is_nan()
-            } else {
-                pow(base, exponent) == expected
-            },
-            "{} ** {} was {} instead of {}",
-            base,
-            exponent,
-            res,
-            expected
-        );
-    }
-
-    fn test_sets_as_base(sets: &[&[f64]], exponent: f64, expected: f64) {
-        sets.iter()
-            .for_each(|s| s.iter().for_each(|val| pow_test(*val, exponent, expected)));
-    }
-
-    fn test_sets_as_exponent(base: f64, sets: &[&[f64]], expected: f64) {
-        sets.iter()
-            .for_each(|s| s.iter().for_each(|val| pow_test(base, *val, expected)));
-    }
-
-    fn test_sets(sets: &[&[f64]], computed: &dyn Fn(f64) -> f64, expected: &dyn Fn(f64) -> f64) {
-        sets.iter().for_each(|s| {
-            s.iter().for_each(|val| {
-                let exp = expected(*val);
-                let res = computed(*val);
-
-                assert!(
-                    if exp.is_nan() {
-                        res.is_nan()
-                    } else {
-                        exp == res
-                    },
-                    "test for {} was {} instead of {}",
-                    val,
-                    res,
-                    exp
-                );
-            })
-        });
-    }
-
-    #[test]
-    fn zero_as_exponent() {
-        test_sets_as_base(ALL, 0.0, 1.0);
-        test_sets_as_base(ALL, -0.0, 1.0);
-    }
-
-    #[test]
-    fn one_as_base() {
-        test_sets_as_exponent(1.0, ALL, 1.0);
-    }
-
-    #[test]
-    fn nan_inputs() {
-        // NAN as the base:
-        // (NAN ^ anything *but 0* should be NAN)
-        test_sets_as_exponent(NAN, &ALL[2..], NAN);
-
-        // NAN as the exponent:
-        // (anything *but 1* ^ NAN should be NAN)
-        test_sets_as_base(&ALL[..(ALL.len() - 2)], NAN, NAN);
-    }
-
-    #[test]
-    fn infinity_as_base() {
-        // Positive Infinity as the base:
-        // (+Infinity ^ positive anything but 0 and NAN should be +Infinity)
-        test_sets_as_exponent(INFINITY, &POS[1..], INFINITY);
-
-        // (+Infinity ^ negative anything except 0 and NAN should be 0.0)
-        test_sets_as_exponent(INFINITY, &NEG[1..], 0.0);
-
-        // Negative Infinity as the base:
-        // (-Infinity ^ positive odd ints should be -Infinity)
-        test_sets_as_exponent(NEG_INFINITY, &[POS_ODDS], NEG_INFINITY);
-
-        // (-Infinity ^ anything but odd ints should be == -0 ^ (-anything))
-        // We can lump in pos/neg odd ints here because they don't seem to
-        // cause panics (div by zero) in release mode (I think).
-        test_sets(ALL, &|v: f64| pow(NEG_INFINITY, v), &|v: f64| pow(-0.0, -v));
-    }
-
-    #[test]
-    fn infinity_as_exponent() {
-        // Positive/Negative base greater than 1:
-        // (pos/neg > 1 ^ Infinity should be Infinity - note this excludes NAN as the base)
-        test_sets_as_base(&ALL[5..(ALL.len() - 2)], INFINITY, INFINITY);
-
-        // (pos/neg > 1 ^ -Infinity should be 0.0)
-        test_sets_as_base(&ALL[5..ALL.len() - 2], NEG_INFINITY, 0.0);
-
-        // Positive/Negative base less than 1:
-        let base_below_one = &[POS_ZERO, NEG_ZERO, NEG_SMALL_FLOATS, POS_SMALL_FLOATS];
-
-        // (pos/neg < 1 ^ Infinity should be 0.0 - this also excludes NAN as the base)
-        test_sets_as_base(base_below_one, INFINITY, 0.0);
-
-        // (pos/neg < 1 ^ -Infinity should be Infinity)
-        test_sets_as_base(base_below_one, NEG_INFINITY, INFINITY);
-
-        // Positive/Negative 1 as the base:
-        // (pos/neg 1 ^ Infinity should be 1)
-        test_sets_as_base(&[NEG_ONE, POS_ONE], INFINITY, 1.0);
-
-        // (pos/neg 1 ^ -Infinity should be 1)
-        test_sets_as_base(&[NEG_ONE, POS_ONE], NEG_INFINITY, 1.0);
-    }
-
-    #[test]
-    fn zero_as_base() {
-        // Positive Zero as the base:
-        // (+0 ^ anything positive but 0 and NAN should be +0)
-        test_sets_as_exponent(0.0, &POS[1..], 0.0);
-
-        // (+0 ^ anything negative but 0 and NAN should be Infinity)
-        // (this should panic because we're dividing by zero)
-        test_sets_as_exponent(0.0, &NEG[1..], INFINITY);
-
-        // Negative Zero as the base:
-        // (-0 ^ anything positive but 0, NAN, and odd ints should be +0)
-        test_sets_as_exponent(-0.0, &POS[3..], 0.0);
-
-        // (-0 ^ anything negative but 0, NAN, and odd ints should be Infinity)
-        // (should panic because of divide by zero)
-        test_sets_as_exponent(-0.0, &NEG[3..], INFINITY);
-
-        // (-0 ^ positive odd ints should be -0)
-        test_sets_as_exponent(-0.0, &[POS_ODDS], -0.0);
-
-        // (-0 ^ negative odd ints should be -Infinity)
-        // (should panic because of divide by zero)
-        test_sets_as_exponent(-0.0, &[NEG_ODDS], NEG_INFINITY);
-    }
-
-    #[test]
-    fn special_cases() {
-        // One as the exponent:
-        // (anything ^ 1 should be anything - i.e. the base)
-        test_sets(ALL, &|v: f64| pow(v, 1.0), &|v: f64| v);
-
-        // Negative One as the exponent:
-        // (anything ^ -1 should be 1/anything)
-        test_sets(ALL, &|v: f64| pow(v, -1.0), &|v: f64| 1.0 / v);
-
-        // Factoring -1 out:
-        // (negative anything ^ integer should be (-1 ^ integer) * (positive anything ^ integer))
-        &[POS_ZERO, NEG_ZERO, POS_ONE, NEG_ONE, POS_EVENS, NEG_EVENS]
-            .iter()
-            .for_each(|int_set| {
-                int_set.iter().for_each(|int| {
-                    test_sets(ALL, &|v: f64| pow(-v, *int), &|v: f64| {
-                        pow(-1.0, *int) * pow(v, *int)
-                    });
-                })
-            });
-
-        // Negative base (imaginary results):
-        // (-anything except 0 and Infinity ^ non-integer should be NAN)
-        &NEG[1..(NEG.len() - 1)].iter().for_each(|set| {
-            set.iter().for_each(|val| {
-                test_sets(&ALL[3..7], &|v: f64| pow(*val, v), &|_| NAN);
-            })
-        });
-    }
-
-    #[test]
-    fn normal_cases() {
-        assert_eq!(pow(2.0, 20.0), (1 << 20) as f64);
-        assert_eq!(pow(-1.0, 9.0), -1.0);
-        assert!(pow(-1.0, 2.2).is_nan());
-        assert!(pow(-1.0, -1.14).is_nan());
-    }
-}
diff --git a/src/math/powf.rs b/crates/libm/src/math/powf.rs
similarity index 99%
rename from src/math/powf.rs
rename to crates/libm/src/math/powf.rs
index 015bade86..45ab1e4de 100644
--- a/src/math/powf.rs
+++ b/crates/libm/src/math/powf.rs
@@ -45,7 +45,7 @@ const IVLN2_L: f32 = 7.0526075433e-06;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn powf(x: f32, y: f32) -> f32 {
+pub extern "C" fn powf(x: f32, y: f32) -> f32 {
     let mut z: f32;
     let mut ax: f32;
     let z_h: f32;
diff --git a/src/math/rem_pio2.rs b/crates/libm/src/math/rem_pio2.rs
similarity index 100%
rename from src/math/rem_pio2.rs
rename to crates/libm/src/math/rem_pio2.rs
diff --git a/src/math/rem_pio2_large.rs b/crates/libm/src/math/rem_pio2_large.rs
similarity index 100%
rename from src/math/rem_pio2_large.rs
rename to crates/libm/src/math/rem_pio2_large.rs
diff --git a/src/math/rem_pio2f.rs b/crates/libm/src/math/rem_pio2f.rs
similarity index 100%
rename from src/math/rem_pio2f.rs
rename to crates/libm/src/math/rem_pio2f.rs
diff --git a/src/math/remainder.rs b/crates/libm/src/math/remainder.rs
similarity index 70%
rename from src/math/remainder.rs
rename to crates/libm/src/math/remainder.rs
index 7ce895004..57fd5aa38 100644
--- a/src/math/remainder.rs
+++ b/crates/libm/src/math/remainder.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn remainder(x: f64, y: f64) -> f64 {
+pub extern "C" fn remainder(x: f64, y: f64) -> f64 {
     let (result, _) = super::remquo(x, y);
     result
 }
diff --git a/src/math/remainderf.rs b/crates/libm/src/math/remainderf.rs
similarity index 70%
rename from src/math/remainderf.rs
rename to crates/libm/src/math/remainderf.rs
index 8b2aa5aab..3f5f92700 100644
--- a/src/math/remainderf.rs
+++ b/crates/libm/src/math/remainderf.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn remainderf(x: f32, y: f32) -> f32 {
+pub extern "C" fn remainderf(x: f32, y: f32) -> f32 {
     let (result, _) = super::remquof(x, y);
     result
 }
diff --git a/src/math/remquo.rs b/crates/libm/src/math/remquo.rs
similarity index 88%
rename from src/math/remquo.rs
rename to crates/libm/src/math/remquo.rs
index c72c8f187..480216baf 100644
--- a/src/math/remquo.rs
+++ b/crates/libm/src/math/remquo.rs
@@ -1,4 +1,6 @@
-pub fn remquo(mut x: f64, mut y: f64) -> (f64, i32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn remquo(mut x: f64, mut y: f64) -> (f64, i32) {
     let ux: u64 = x.to_bits();
     let mut uy: u64 = y.to_bits();
     let mut ex = ((ux >> 52) & 0x7ff) as i32;
@@ -96,14 +98,3 @@ pub fn remquo(mut x: f64, mut y: f64) -> (f64, i32) {
         (x, quo)
     }
 }
-
-#[cfg(test)]
-mod tests {
-    use super::remquo;
-
-    #[test]
-    fn test_q_overflow() {
-        // 0xc000000000000001, 0x04c0000000000004
-        let _ = remquo(-2.0000000000000004, 8.406091369059082e-286);
-    }
-}
diff --git a/src/math/remquof.rs b/crates/libm/src/math/remquof.rs
similarity index 93%
rename from src/math/remquof.rs
rename to crates/libm/src/math/remquof.rs
index 871d0c7d6..b2b732a96 100644
--- a/src/math/remquof.rs
+++ b/crates/libm/src/math/remquof.rs
@@ -1,4 +1,6 @@
-pub fn remquof(mut x: f32, mut y: f32) -> (f32, i32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn remquof(mut x: f32, mut y: f32) -> (f32, i32) {
     let ux: u32 = x.to_bits();
     let mut uy: u32 = y.to_bits();
     let mut ex = ((ux >> 23) & 0xff) as i32;
diff --git a/src/math/round.rs b/crates/libm/src/math/round.rs
similarity index 76%
rename from src/math/round.rs
rename to crates/libm/src/math/round.rs
index 67590d2c1..d467adb0e 100644
--- a/src/math/round.rs
+++ b/crates/libm/src/math/round.rs
@@ -4,7 +4,7 @@ const TOINT: f64 = 1.0 / f64::EPSILON;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn round(mut x: f64) -> f64 {
+pub extern "C" fn round(mut x: f64) -> f64 {
     let i = x.to_bits();
     let e: u64 = i >> 52 & 0x7ff;
     let mut y: f64;
@@ -35,13 +35,3 @@ pub fn round(mut x: f64) -> f64 {
         y
     }
 }
-
-#[cfg(test)]
-mod tests {
-    use super::round;
-
-    #[test]
-    fn negative_zero() {
-        assert_eq!(round(-0.0_f64).to_bits(), (-0.0_f64).to_bits());
-    }
-}
diff --git a/src/math/roundf.rs b/crates/libm/src/math/roundf.rs
similarity index 75%
rename from src/math/roundf.rs
rename to crates/libm/src/math/roundf.rs
index 85114be4b..cb03ac7b5 100644
--- a/src/math/roundf.rs
+++ b/crates/libm/src/math/roundf.rs
@@ -4,7 +4,7 @@ const TOINT: f32 = 1.0 / f32::EPSILON;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn roundf(mut x: f32) -> f32 {
+pub extern "C" fn roundf(mut x: f32) -> f32 {
     let i = x.to_bits();
     let e: u32 = i >> 23 & 0xff;
     let mut y: f32;
@@ -33,13 +33,3 @@ pub fn roundf(mut x: f32) -> f32 {
         y
     }
 }
-
-#[cfg(test)]
-mod tests {
-    use super::roundf;
-
-    #[test]
-    fn negative_zero() {
-        assert_eq!(roundf(-0.0_f32).to_bits(), (-0.0_f32).to_bits());
-    }
-}
diff --git a/src/math/scalbn.rs b/crates/libm/src/math/scalbn.rs
similarity index 94%
rename from src/math/scalbn.rs
rename to crates/libm/src/math/scalbn.rs
index d8c8409ac..03072f934 100644
--- a/src/math/scalbn.rs
+++ b/crates/libm/src/math/scalbn.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn scalbn(x: f64, mut n: i32) -> f64 {
+pub extern "C" fn scalbn(x: f64, mut n: i32) -> f64 {
     let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
     let x1p53 = f64::from_bits(0x4340000000000000); // 0x1p53 === 2 ^ 53
     let x1p_1022 = f64::from_bits(0x0010000000000000); // 0x1p-1022 === 2 ^ (-1022)
diff --git a/src/math/scalbnf.rs b/crates/libm/src/math/scalbnf.rs
similarity index 92%
rename from src/math/scalbnf.rs
rename to crates/libm/src/math/scalbnf.rs
index 4e9771175..b11af2544 100644
--- a/src/math/scalbnf.rs
+++ b/crates/libm/src/math/scalbnf.rs
@@ -1,6 +1,6 @@
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn scalbnf(mut x: f32, mut n: i32) -> f32 {
+pub extern "C" fn scalbnf(mut x: f32, mut n: i32) -> f32 {
     let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
     let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126
     let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
diff --git a/src/math/sin.rs b/crates/libm/src/math/sin.rs
similarity index 91%
rename from src/math/sin.rs
rename to crates/libm/src/math/sin.rs
index 51aed88a8..8af0f0263 100644
--- a/src/math/sin.rs
+++ b/crates/libm/src/math/sin.rs
@@ -42,7 +42,7 @@ use super::{k_cos, k_sin, rem_pio2};
 //      TRIG(x) returns trig(x) nearly rounded
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn sin(x: f64) -> f64 {
+pub extern "C" fn sin(x: f64) -> f64 {
     let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120
 
     /* High word of x. */
@@ -77,10 +77,3 @@ pub fn sin(x: f64) -> f64 {
         _ => -k_cos(y0, y1),
     }
 }
-
-#[test]
-fn test_near_pi() {
-    let x = f64::from_bits(0x400921fb000FD5DD); // 3.141592026217707
-    let sx = f64::from_bits(0x3ea50d15ced1a4a2); // 6.273720864039205e-7
-    assert_eq!(sin(x), sx);
-}
diff --git a/src/math/sincos.rs b/crates/libm/src/math/sincos.rs
similarity index 93%
rename from src/math/sincos.rs
rename to crates/libm/src/math/sincos.rs
index 750908df4..c06ce7db0 100644
--- a/src/math/sincos.rs
+++ b/crates/libm/src/math/sincos.rs
@@ -12,7 +12,9 @@
 
 use super::{get_high_word, k_cos, k_sin, rem_pio2};
 
-pub fn sincos(x: f64) -> (f64, f64) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn sincos(x: f64) -> (f64, f64) {
     let s: f64;
     let c: f64;
     let mut ix: u32;
diff --git a/src/math/sincosf.rs b/crates/libm/src/math/sincosf.rs
similarity index 96%
rename from src/math/sincosf.rs
rename to crates/libm/src/math/sincosf.rs
index bb9a00392..7d92dc10d 100644
--- a/src/math/sincosf.rs
+++ b/crates/libm/src/math/sincosf.rs
@@ -23,7 +23,9 @@ const S2PIO2: f32 = 2.0 * PI_2; /* 0x400921FB, 0x54442D18 */
 const S3PIO2: f32 = 3.0 * PI_2; /* 0x4012D97C, 0x7F3321D2 */
 const S4PIO2: f32 = 4.0 * PI_2; /* 0x401921FB, 0x54442D18 */
 
-pub fn sincosf(x: f32) -> (f32, f32) {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn sincosf(x: f32) -> (f32, f32) {
     let s: f32;
     let c: f32;
     let mut ix: u32;
diff --git a/src/math/sinf.rs b/crates/libm/src/math/sinf.rs
similarity index 98%
rename from src/math/sinf.rs
rename to crates/libm/src/math/sinf.rs
index 0c31099ed..9b91ebb5a 100644
--- a/src/math/sinf.rs
+++ b/crates/libm/src/math/sinf.rs
@@ -26,7 +26,7 @@ const S4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn sinf(x: f32) -> f32 {
+pub extern "C" fn sinf(x: f32) -> f32 {
     let x64 = x as f64;
 
     let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
diff --git a/src/math/sinh.rs b/crates/libm/src/math/sinh.rs
similarity index 96%
rename from src/math/sinh.rs
rename to crates/libm/src/math/sinh.rs
index d36de66c1..0c09f32d2 100644
--- a/src/math/sinh.rs
+++ b/crates/libm/src/math/sinh.rs
@@ -6,7 +6,7 @@ use super::{expm1, expo2};
 //
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn sinh(x: f64) -> f64 {
+pub extern "C" fn sinh(x: f64) -> f64 {
     // union {double f; uint64_t i;} u = {.f = x};
     // uint32_t w;
     // double t, h, absx;
diff --git a/src/math/sinhf.rs b/crates/libm/src/math/sinhf.rs
similarity index 94%
rename from src/math/sinhf.rs
rename to crates/libm/src/math/sinhf.rs
index fd0b2bfc8..b553c6b7b 100644
--- a/src/math/sinhf.rs
+++ b/crates/libm/src/math/sinhf.rs
@@ -3,7 +3,7 @@ use super::k_expo2f;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn sinhf(x: f32) -> f32 {
+pub extern "C" fn sinhf(x: f32) -> f32 {
     let mut h = 0.5f32;
     let mut ix = x.to_bits();
     if (ix >> 31) != 0 {
diff --git a/src/math/sqrt.rs b/crates/libm/src/math/sqrt.rs
similarity index 98%
rename from src/math/sqrt.rs
rename to crates/libm/src/math/sqrt.rs
index 14404d4eb..e71fcde9c 100644
--- a/src/math/sqrt.rs
+++ b/crates/libm/src/math/sqrt.rs
@@ -83,12 +83,12 @@ const TINY: f64 = 1.0e-300;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn sqrt(x: f64) -> f64 {
+pub extern "C" fn sqrt(x: f64) -> f64 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f64.sqrt` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return if x < 0.0 {
                 f64::NAN
             } else {
diff --git a/src/math/sqrtf.rs b/crates/libm/src/math/sqrtf.rs
similarity index 96%
rename from src/math/sqrtf.rs
rename to crates/libm/src/math/sqrtf.rs
index b9365c617..c6e978ccc 100644
--- a/src/math/sqrtf.rs
+++ b/crates/libm/src/math/sqrtf.rs
@@ -17,12 +17,12 @@ const TINY: f32 = 1.0e-30;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn sqrtf(x: f32) -> f32 {
+pub extern "C" fn sqrtf(x: f32) -> f32 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f32.sqrt` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return if x < 0.0 {
                 ::core::f32::NAN
             } else {
diff --git a/src/math/tan.rs b/crates/libm/src/math/tan.rs
similarity index 98%
rename from src/math/tan.rs
rename to crates/libm/src/math/tan.rs
index e5c94cbb1..eac452070 100644
--- a/src/math/tan.rs
+++ b/crates/libm/src/math/tan.rs
@@ -41,7 +41,7 @@ use super::{k_tan, rem_pio2};
 //      TRIG(x) returns trig(x) nearly rounded
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn tan(x: f64) -> f64 {
+pub extern "C" fn tan(x: f64) -> f64 {
     let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
 
     let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
diff --git a/src/math/tanf.rs b/crates/libm/src/math/tanf.rs
similarity index 98%
rename from src/math/tanf.rs
rename to crates/libm/src/math/tanf.rs
index c286cdeb4..1e593b116 100644
--- a/src/math/tanf.rs
+++ b/crates/libm/src/math/tanf.rs
@@ -26,7 +26,7 @@ const T4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn tanf(x: f32) -> f32 {
+pub extern "C" fn tanf(x: f32) -> f32 {
     let x64 = x as f64;
 
     let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
diff --git a/src/math/tanh.rs b/crates/libm/src/math/tanh.rs
similarity index 96%
rename from src/math/tanh.rs
rename to crates/libm/src/math/tanh.rs
index 75d695cf7..4cad2a998 100644
--- a/src/math/tanh.rs
+++ b/crates/libm/src/math/tanh.rs
@@ -6,7 +6,7 @@ use super::expm1;
  */
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn tanh(mut x: f64) -> f64 {
+pub extern "C" fn tanh(mut x: f64) -> f64 {
     let mut uf: f64 = x;
     let mut ui: u64 = f64::to_bits(uf);
 
diff --git a/src/math/tanhf.rs b/crates/libm/src/math/tanhf.rs
similarity index 95%
rename from src/math/tanhf.rs
rename to crates/libm/src/math/tanhf.rs
index ac4657b5a..ca051af4c 100644
--- a/src/math/tanhf.rs
+++ b/crates/libm/src/math/tanhf.rs
@@ -2,7 +2,7 @@ use super::expm1f;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn tanhf(mut x: f32) -> f32 {
+pub extern "C" fn tanhf(mut x: f32) -> f32 {
     /* x = |x| */
     let mut ix = x.to_bits();
     let sign = (ix >> 31) != 0;
diff --git a/src/math/tgamma.rs b/crates/libm/src/math/tgamma.rs
similarity index 97%
rename from src/math/tgamma.rs
rename to crates/libm/src/math/tgamma.rs
index f8ccf669a..0e0de0ab0 100644
--- a/src/math/tgamma.rs
+++ b/crates/libm/src/math/tgamma.rs
@@ -130,7 +130,9 @@ fn s(x: f64) -> f64 {
     return num / den;
 }
 
-pub fn tgamma(mut x: f64) -> f64 {
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn tgamma(mut x: f64) -> f64 {
     let u: u64 = x.to_bits();
     let absx: f64;
     let mut y: f64;
diff --git a/crates/libm/src/math/tgammaf.rs b/crates/libm/src/math/tgammaf.rs
new file mode 100644
index 000000000..cc9529290
--- /dev/null
+++ b/crates/libm/src/math/tgammaf.rs
@@ -0,0 +1,7 @@
+use super::tgamma;
+
+#[inline]
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub extern "C" fn tgammaf(x: f32) -> f32 {
+    tgamma(x as f64) as f32
+}
diff --git a/src/math/trunc.rs b/crates/libm/src/math/trunc.rs
similarity index 80%
rename from src/math/trunc.rs
rename to crates/libm/src/math/trunc.rs
index 1ee46fc7d..0d8387a74 100644
--- a/src/math/trunc.rs
+++ b/crates/libm/src/math/trunc.rs
@@ -2,12 +2,12 @@ use core::f64;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn trunc(x: f64) -> f64 {
+pub extern "C" fn trunc(x: f64) -> f64 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f64.trunc` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::truncf64(x) }
         }
     }
@@ -31,11 +31,3 @@ pub fn trunc(x: f64) -> f64 {
     i &= !m;
     f64::from_bits(i)
 }
-
-#[cfg(test)]
-mod tests {
-    #[test]
-    fn sanity_check() {
-        assert_eq!(super::trunc(1.1), 1.0);
-    }
-}
diff --git a/src/math/truncf.rs b/crates/libm/src/math/truncf.rs
similarity index 80%
rename from src/math/truncf.rs
rename to crates/libm/src/math/truncf.rs
index f93383269..5c3646741 100644
--- a/src/math/truncf.rs
+++ b/crates/libm/src/math/truncf.rs
@@ -2,12 +2,12 @@ use core::f32;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn truncf(x: f32) -> f32 {
+pub extern "C" fn truncf(x: f32) -> f32 {
     // On wasm32 we know that LLVM's intrinsic will compile to an optimized
     // `f32.trunc` native instruction, so we can leverage this for both code size
     // and speed.
     llvm_intrinsically_optimized! {
-        #[cfg(target_arch = "wasm32")] {
+        #[cfg(all(target_arch = "wasm32", not(feature = "stable")))] {
             return unsafe { ::core::intrinsics::truncf32(x) }
         }
     }
@@ -31,11 +31,3 @@ pub fn truncf(x: f32) -> f32 {
     i &= !m;
     f32::from_bits(i)
 }
-
-#[cfg(test)]
-mod tests {
-    #[test]
-    fn sanity_check() {
-        assert_eq!(super::truncf(1.1), 1.0);
-    }
-}
diff --git a/src/lib.rs b/src/lib.rs
deleted file mode 100644
index b15857dbe..000000000
--- a/src/lib.rs
+++ /dev/null
@@ -1,49 +0,0 @@
-//! libm in pure Rust
-#![deny(warnings)]
-#![no_std]
-#![cfg_attr(
-    all(target_arch = "wasm32", not(feature = "stable")),
-    feature(core_intrinsics)
-)]
-
-mod math;
-
-use core::{f32, f64};
-
-pub use self::math::*;
-
-/// Approximate equality with 1 ULP of tolerance
-#[doc(hidden)]
-#[inline]
-pub fn _eqf(a: f32, b: f32) -> Result<(), u32> {
-    if a.is_nan() && b.is_nan() {
-        Ok(())
-    } else {
-        let err = (a.to_bits() as i32).wrapping_sub(b.to_bits() as i32).abs();
-
-        if err <= 1 {
-            Ok(())
-        } else {
-            Err(err as u32)
-        }
-    }
-}
-
-#[doc(hidden)]
-#[inline]
-pub fn _eq(a: f64, b: f64) -> Result<(), u64> {
-    if a.is_nan() && b.is_nan() {
-        Ok(())
-    } else {
-        let err = (a.to_bits() as i64).wrapping_sub(b.to_bits() as i64).abs();
-
-        if err <= 1 {
-            Ok(())
-        } else {
-            Err(err as u64)
-        }
-    }
-}
-
-#[cfg(all(test, feature = "musl-reference-tests"))]
-include!(concat!(env!("OUT_DIR"), "/musl-tests.rs"));
diff --git a/src/math/fenv.rs b/src/math/fenv.rs
deleted file mode 100644
index 63bb20368..000000000
--- a/src/math/fenv.rs
+++ /dev/null
@@ -1,33 +0,0 @@
-// src: musl/src/fenv/fenv.c
-/* Dummy functions for archs lacking fenv implementation */
-
-pub const FE_UNDERFLOW: i32 = 0;
-pub const FE_INEXACT: i32 = 0;
-
-pub const FE_TONEAREST: i32 = 0;
-pub const FE_TOWARDZERO: i32 = 0;
-
-#[inline]
-pub fn feclearexcept(_mask: i32) -> i32 {
-    0
-}
-
-#[inline]
-pub fn feraiseexcept(_mask: i32) -> i32 {
-    0
-}
-
-#[inline]
-pub fn fetestexcept(_mask: i32) -> i32 {
-    0
-}
-
-#[inline]
-pub fn fegetround() -> i32 {
-    FE_TONEAREST
-}
-
-#[inline]
-pub fn fesetround(_r: i32) -> i32 {
-    0
-}
diff --git a/src/math/lgamma.rs b/src/math/lgamma.rs
deleted file mode 100644
index 5bc87e85e..000000000
--- a/src/math/lgamma.rs
+++ /dev/null
@@ -1,5 +0,0 @@
-use super::lgamma_r;
-
-pub fn lgamma(x: f64) -> f64 {
-    lgamma_r(x).0
-}
diff --git a/src/math/lgammaf.rs b/src/math/lgammaf.rs
deleted file mode 100644
index dfdc87f96..000000000
--- a/src/math/lgammaf.rs
+++ /dev/null
@@ -1,5 +0,0 @@
-use super::lgammaf_r;
-
-pub fn lgammaf(x: f32) -> f32 {
-    lgammaf_r(x).0
-}
diff --git a/src/math/tgammaf.rs b/src/math/tgammaf.rs
deleted file mode 100644
index a8f161f0c..000000000
--- a/src/math/tgammaf.rs
+++ /dev/null
@@ -1,5 +0,0 @@
-use super::tgamma;
-
-pub fn tgammaf(x: f32) -> f32 {
-    tgamma(x as f64) as f32
-}