diff --git a/CHANGELOG.md b/CHANGELOG.md
index f8c7722..e2793a9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,10 @@ Possible sections are:
 <!-- next-header -->
 ## [Unreleased] - ReleaseDate
 
+### Added
+
+- add `no_std` support (`std` is still enabled by default though)
+
 ## [0.5.1] - 2024-01-27
 
 ### Fixed
diff --git a/Cargo.toml b/Cargo.toml
index 060f2e3..cb92152 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -10,7 +10,7 @@ homepage = "https://docs.rs/h3o"
 repository = "https://github.com/HydroniumLabs/h3o"
 license = "BSD-3-Clause"
 keywords = ["geography", "geospatial", "gis", "h3", "spatial-index"]
-categories = ["science::geo"]
+categories = ["science::geo", "no-std"]
 
 [package.metadata.docs.rs]
 all-features = true
@@ -24,21 +24,22 @@ pre-release-replacements = [
 ]
 
 [features]
-default = []
+default = ["std"]
+std = ["dep:ahash"]
 geo = ["dep:geo", "dep:geojson"]
 serde = ["dep:serde", "dep:serde_repr"]
 tools = ["polyfit-rs"]
 
 [dependencies]
-ahash = { version = "0.8", default-features = false, features = ["std", "compile-time-rng"] }
+ahash = { version = "0.8", optional = true, default-features = false, features = ["std", "compile-time-rng"] }
 arbitrary = { version = "1.0", optional = true, default-features = false }
-auto_ops = { version = "0.3", default-features = false }
 konst = { version = "0.3", default-features = false, features = ["parsing"] }
 either = { version = "1.0", default-features = false }
 float_eq = { version = "1.0", default-features = false }
 geo = { version = "0.27", optional = true, default-features = false }
 geojson = { version = "0.24", optional = true, default-features = false, features = ["geo-types"] }
 h3o-bit = { version = "0.1", default-features = false }
+libm = { version = "0.2", default-features = false }
 polyfit-rs = { version = "0.2", optional = true, default-features = false }
 serde = { version = "1.0", optional = true, default-features = false, features = ["derive"] }
 serde_repr = { version = "0.1", optional = true, default-features = false }
diff --git a/src/base_cell.rs b/src/base_cell.rs
index a5e84cd..cc2ff6b 100644
--- a/src/base_cell.rs
+++ b/src/base_cell.rs
@@ -2,7 +2,7 @@ use crate::{
     coord::{CoordIJK, FaceIJK},
     error, Direction, Face, NUM_PENTAGONS, NUM_PENT_VERTS,
 };
-use std::fmt;
+use core::fmt;
 
 /// Maximum value for a base cell.
 pub const MAX: u8 = 121;
diff --git a/src/boundary.rs b/src/boundary.rs
index 2d72fd3..53e5ce0 100644
--- a/src/boundary.rs
+++ b/src/boundary.rs
@@ -1,5 +1,5 @@
 use crate::LatLng;
-use std::{fmt, ops::Deref};
+use core::{fmt, ops::Deref};
 
 /// Maximum number of cell boundary vertices.
 ///
@@ -42,14 +42,14 @@ impl Deref for Boundary {
 
 impl fmt::Display for Boundary {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        write!(
-            f,
-            "[{}]",
-            self.iter()
-                .map(ToString::to_string)
-                .collect::<Vec<_>>()
-                .join("-")
-        )
+        write!(f, "[",)?;
+        for (i, ll) in self.iter().enumerate() {
+            if i != 0 {
+                write!(f, "-")?;
+            }
+            write!(f, "{ll}")?;
+        }
+        write!(f, "]",)
     }
 }
 
diff --git a/src/coord/cube.rs b/src/coord/cube.rs
index fd8f372..a2a751f 100644
--- a/src/coord/cube.rs
+++ b/src/coord/cube.rs
@@ -1,4 +1,5 @@
 use super::CoordIJK;
+use crate::math::{abs, round};
 
 /// Cube coordinates.
 ///
@@ -29,11 +30,11 @@ impl CoordCube {
 
         #[allow(clippy::cast_possible_truncation)] // on purpose
         let (mut ri, mut rj, mut rk) =
-            { (i.round() as i32, j.round() as i32, k.round() as i32) };
+            { (round(i) as i32, round(j) as i32, round(k) as i32) };
 
-        let i_diff = (f64::from(ri) - i).abs();
-        let j_diff = (f64::from(rj) - j).abs();
-        let k_diff = (f64::from(rk) - k).abs();
+        let i_diff = abs(f64::from(ri) - i);
+        let j_diff = abs(f64::from(rj) - j);
+        let k_diff = abs(f64::from(rk) - k);
 
         // Round, maintaining valid cube coords.
         if i_diff > j_diff && i_diff > k_diff {
diff --git a/src/coord/ijk.rs b/src/coord/ijk.rs
index a400f95..c567d4b 100644
--- a/src/coord/ijk.rs
+++ b/src/coord/ijk.rs
@@ -9,12 +9,16 @@
 //! a unique address consisting of the minimal positive `IJK` components; this
 //! always results in at most two non-zero components.
 
-#![allow(clippy::use_self)] // False positive with `auto_ops::impl_op_ex`
-
 use super::{CoordCube, Vec2d, SQRT3_2};
-use crate::{error::HexGridError, Direction};
-use auto_ops::impl_op_ex;
-use std::{cmp, fmt};
+use crate::{
+    error::HexGridError,
+    math::{mul_add, round},
+    Direction,
+};
+use core::{
+    cmp, fmt,
+    ops::{Add, MulAssign, Sub},
+};
 
 // -----------------------------------------------------------------------------
 
@@ -116,7 +120,7 @@ impl CoordIJK {
     }
 
     pub fn distance(&self, other: &Self) -> i32 {
-        let diff = (self - other).normalize();
+        let diff = (*self - *other).normalize();
 
         cmp::max(diff.i.abs(), cmp::max(diff.j.abs(), diff.k.abs()))
     }
@@ -133,7 +137,7 @@ impl CoordIJK {
             (f64::from(2 * i + j) / 7., f64::from(3 * j - i) / 7.)
         };
 
-        Self::new(i.round() as i32, j.round() as i32, 0).normalize()
+        Self::new(round(i) as i32, round(j) as i32, 0).normalize()
     }
 
     /// Returns the normalized `IJK` coordinates of the indexing parent of a
@@ -154,7 +158,7 @@ impl CoordIJK {
             )
         };
 
-        Self::new(i.round() as i32, j.round() as i32, 0).checked_normalize()
+        Self::new(round(i) as i32, round(j) as i32, 0).checked_normalize()
     }
 
     /// Returns the normalized `IJK` coordinates of the hex centered on the
@@ -194,7 +198,7 @@ impl CoordIJK {
     /// Returns the normalized `IJK` coordinates of the hex in the specified
     /// direction from the current position.
     pub fn neighbor(&self, direction: Direction) -> Self {
-        (self + direction.coordinate()).normalize()
+        (*self + direction.coordinate()).normalize()
     }
 
     /// Returns the `IJK` coordinates after a 60 degrees rotation.
@@ -213,27 +217,41 @@ impl CoordIJK {
     }
 }
 
-impl_op_ex!(+ |lhs: &CoordIJK, rhs: &CoordIJK| -> CoordIJK {
-    CoordIJK{
-        i: lhs.i + rhs.i,
-        j: lhs.j + rhs.j,
-        k: lhs.k + rhs.k,
+// -----------------------------------------------------------------------------
+
+impl Add for CoordIJK {
+    type Output = Self;
+
+    fn add(self, other: Self) -> Self {
+        Self {
+            i: self.i + other.i,
+            j: self.j + other.j,
+            k: self.k + other.k,
+        }
     }
-});
+}
+
+impl Sub for CoordIJK {
+    type Output = Self;
 
-impl_op_ex!(-|lhs: &CoordIJK, rhs: &CoordIJK| -> CoordIJK {
-    CoordIJK {
-        i: lhs.i - rhs.i,
-        j: lhs.j - rhs.j,
-        k: lhs.k - rhs.k,
+    fn sub(self, other: Self) -> Self {
+        Self {
+            i: self.i - other.i,
+            j: self.j - other.j,
+            k: self.k - other.k,
+        }
+    }
+}
+
+impl MulAssign<i32> for CoordIJK {
+    fn mul_assign(&mut self, rhs: i32) {
+        self.i *= rhs;
+        self.j *= rhs;
+        self.k *= rhs;
     }
-});
+}
 
-impl_op_ex!(*= |lhs: &mut CoordIJK, rhs: i32| {
-    lhs.i *= rhs;
-    lhs.j *= rhs;
-    lhs.k *= rhs;
-});
+// -----------------------------------------------------------------------------
 
 impl From<CoordIJK> for Vec2d {
     // Returns the center point in 2D cartesian coordinates of a hex.
@@ -241,7 +259,7 @@ impl From<CoordIJK> for Vec2d {
         let i = f64::from(value.i - value.k);
         let j = f64::from(value.j - value.k);
 
-        Self::new(0.5_f64.mul_add(-j, i), j * SQRT3_2)
+        Self::new(mul_add(0.5, -j, i), j * SQRT3_2)
     }
 }
 
diff --git a/src/coord/latlng.rs b/src/coord/latlng.rs
index 0662781..7c8e3ce 100644
--- a/src/coord/latlng.rs
+++ b/src/coord/latlng.rs
@@ -3,14 +3,16 @@ use super::{
     RES0_U_GNOMONIC, SQRT7_POWERS,
 };
 use crate::{
-    error::InvalidLatLng, face, CellIndex, Face, Resolution, EARTH_RADIUS_KM,
-    TWO_PI,
+    error::InvalidLatLng,
+    face,
+    math::{acos, asin, atan2, cos, mul_add, sin, sqrt, tan},
+    CellIndex, Face, Resolution, EARTH_RADIUS_KM, TWO_PI,
 };
-use float_eq::float_eq;
-use std::{
+use core::{
     f64::consts::{FRAC_PI_2, PI},
     fmt,
 };
+use float_eq::float_eq;
 
 /// Epsilon of ~0.1mm in degrees.
 const EPSILON_DEG: f64 = 0.000000001;
@@ -153,15 +155,16 @@ impl LatLng {
     /// ```
     #[must_use]
     pub fn distance_rads(self, other: Self) -> f64 {
-        let sin_lat = ((other.lat - self.lat) / 2.).sin();
-        let sin_lng = ((other.lng - self.lng) / 2.).sin();
+        let sin_lat = sin((other.lat - self.lat) / 2.);
+        let sin_lng = sin((other.lng - self.lng) / 2.);
 
-        let a = sin_lat.mul_add(
+        let a = mul_add(
+            sin_lat,
             sin_lat,
-            self.lat.cos() * other.lat.cos() * sin_lng * sin_lng,
+            cos(self.lat) * cos(other.lat) * sin_lng * sin_lng,
         );
 
-        2. * a.sqrt().atan2((1. - a).sqrt())
+        2. * atan2(sqrt(a), sqrt(1. - a))
     }
 
     /// The great circle distance, in kilometers, between two spherical
@@ -231,7 +234,7 @@ impl LatLng {
 
         let r = {
             // cos(r) = 1 - 2 * sin^2(r/2) = 1 - 2 * (sqd / 4) = 1 - sqd/2
-            let r = (1. - distance / 2.).acos();
+            let r = acos(1. - distance / 2.);
 
             if r < EPSILON {
                 return Vec2d::new(0., 0.);
@@ -239,7 +242,7 @@ impl LatLng {
 
             // Perform gnomonic scaling of `r` (`tan(r)`) and scale for current
             // resolution length `u`.
-            (r.tan() / RES0_U_GNOMONIC) * SQRT7_POWERS[usize::from(resolution)]
+            (tan(r) / RES0_U_GNOMONIC) * SQRT7_POWERS[usize::from(resolution)]
         };
 
         let theta = {
@@ -255,7 +258,7 @@ impl LatLng {
         };
 
         // Convert to local x, y.
-        Vec2d::new(r * theta.cos(), r * theta.sin())
+        Vec2d::new(r * cos(theta), r * sin(theta))
     }
 
     /// Finds the closest icosahedral face from the current coordinate.
@@ -290,12 +293,12 @@ impl LatLng {
     /// Computes the azimuth to `other` from `self`, in radians.
     #[must_use]
     pub(crate) fn azimuth(self, other: &Self) -> f64 {
-        (other.lat.cos() * (other.lng - self.lng).sin()).atan2(
-            self.lat.cos().mul_add(
-                other.lat.sin(),
-                -self.lat.sin()
-                    * other.lat.cos()
-                    * (other.lng - self.lng).cos(),
+        atan2(
+            cos(other.lat) * sin(other.lng - self.lng),
+            mul_add(
+                cos(self.lat),
+                sin(other.lat),
+                -sin(self.lat) * cos(other.lat) * cos(other.lng - self.lng),
             ),
         )
     }
@@ -319,14 +322,14 @@ impl LatLng {
                 self.lat - distance // Due South.
             }
         } else {
-            self.lat
-                .sin()
-                .mul_add(
-                    distance.cos(),
-                    self.lat.cos() * distance.sin() * azimuth.cos(),
+            asin(
+                mul_add(
+                    sin(self.lat),
+                    cos(distance),
+                    cos(self.lat) * sin(distance) * cos(azimuth),
                 )
-                .clamp(-1., 1.)
-                .asin()
+                .clamp(-1., 1.),
+            )
         };
 
         // Handle poles.
@@ -341,11 +344,11 @@ impl LatLng {
             self.lng
         } else {
             let sinlng =
-                (azimuth.sin() * distance.sin() / lat.cos()).clamp(-1., 1.);
-            let coslng = self.lat.sin().mul_add(-lat.sin(), distance.cos())
-                / self.lat.cos()
-                / lat.cos();
-            self.lng + sinlng.atan2(coslng)
+                (sin(azimuth) * sin(distance) / cos(lat)).clamp(-1., 1.);
+            let coslng = mul_add(sin(self.lat), sin(-lat), cos(distance))
+                / cos(self.lat)
+                / cos(lat);
+            self.lng + atan2(sinlng, coslng)
         };
 
         // XXX: make sure longitudes are in the proper bounds.
@@ -401,11 +404,11 @@ impl From<LatLng> for Vec3d {
     /// Computes the 3D coordinate on unit sphere from the latitude and
     /// longitude.
     fn from(value: LatLng) -> Self {
-        let r = value.lat.cos();
+        let r = cos(value.lat);
 
-        let z = value.lat.sin();
-        let x = value.lng.cos() * r;
-        let y = value.lng.sin() * r;
+        let z = sin(value.lat);
+        let x = cos(value.lng) * r;
+        let y = sin(value.lng) * r;
 
         Self::new(x, y, z)
     }
diff --git a/src/coord/localij.rs b/src/coord/localij.rs
index 738fe02..cde7703 100644
--- a/src/coord/localij.rs
+++ b/src/coord/localij.rs
@@ -22,7 +22,7 @@ use crate::{
     index::bits,
     BaseCell, CellIndex, Direction, Resolution, CCW, CW, DEFAULT_CELL_INDEX,
 };
-use std::{fmt, num::NonZeroU8};
+use core::{fmt, num::NonZeroU8};
 
 // -----------------------------------------------------------------------------
 
diff --git a/src/coord/vec2d.rs b/src/coord/vec2d.rs
index f97abe1..811bda5 100644
--- a/src/coord/vec2d.rs
+++ b/src/coord/vec2d.rs
@@ -14,7 +14,12 @@ use super::{
     to_positive_angle, CoordIJK, AP7_ROT_RADS, EPSILON, RES0_U_GNOMONIC,
     SQRT7_POWERS,
 };
-use crate::{face, resolution::ExtendedResolution, Face, LatLng};
+use crate::{
+    face,
+    math::{abs, atan, atan2, hypot, mul_add},
+    resolution::ExtendedResolution,
+    Face, LatLng,
+};
 use float_eq::float_eq;
 
 /// sin(60')
@@ -48,7 +53,7 @@ impl Vec2d {
 
     /// Calculates the magnitude.
     pub fn magnitude(self) -> f64 {
-        self.x.hypot(self.y)
+        hypot(self.x, self.y)
     }
 
     /// Finds the intersection between two lines.
@@ -65,14 +70,15 @@ impl Vec2d {
             y: line2.1.y - line2.0.y,
         };
 
-        let t = s2
-            .x
-            .mul_add(line1.0.y - line2.0.y, -s2.y * (line1.0.x - line2.0.x))
-            / (-s2.x).mul_add(s1.y, s1.x * s2.y);
+        let t = mul_add(
+            s2.x,
+            line1.0.y - line2.0.y,
+            -s2.y * (line1.0.x - line2.0.x),
+        ) / mul_add(-s2.x, s1.y, s1.x * s2.y);
 
         Self {
-            x: t.mul_add(s1.x, line1.0.x),
-            y: t.mul_add(s1.y, line1.0.y),
+            x: mul_add(t, s1.x, line1.0.x),
+            y: mul_add(t, s1.y, line1.0.y),
         }
     }
 
@@ -111,11 +117,11 @@ impl Vec2d {
             }
 
             // Perform inverse gnomonic scaling of `r`.
-            (r * RES0_U_GNOMONIC).atan()
+            atan(r * RES0_U_GNOMONIC)
         };
 
         let theta = {
-            let mut theta = self.y.atan2(self.x);
+            let mut theta = atan2(self.y, self.x);
 
             // Adjust theta for Class III.
             // If a substrate grid, then it's already adjusted for Class III.
@@ -139,8 +145,8 @@ impl From<Vec2d> for CoordIJK {
         // Quantize into the IJ system and then normalize.
         let k = 0;
 
-        let a1 = value.x.abs();
-        let a2 = value.y.abs();
+        let a1 = abs(value.x);
+        let a2 = abs(value.y);
 
         // First do a reverse conversion.
         let x2 = a2 / SIN60;
@@ -168,8 +174,8 @@ impl From<Vec2d> for CoordIJK {
             }
         } else if r1 < 2. / 3. {
             let j = m2 + i32::from(r2 >= (1. - r1));
-            let i = m1
-                + i32::from(2.0_f64.mul_add(r1, -1.) >= r2 || r2 >= (1. - r1));
+            let i =
+                m1 + i32::from(mul_add(2.0, r1, -1.) >= r2 || r2 >= (1. - r1));
             (i, j)
         } else {
             let i = m1 + 1;
diff --git a/src/coord/vec3d.rs b/src/coord/vec3d.rs
index 1a4efa2..8e61f68 100644
--- a/src/coord/vec3d.rs
+++ b/src/coord/vec3d.rs
@@ -1,3 +1,5 @@
+use crate::math::mul_add;
+
 /// 3D floating-point vector.
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub struct Vec3d {
@@ -21,7 +23,7 @@ impl Vec3d {
         let y_diff = self.y - other.y;
         let z_diff = self.z - other.z;
 
-        x_diff.mul_add(x_diff, y_diff.mul_add(y_diff, z_diff * z_diff))
+        mul_add(x_diff, x_diff, mul_add(y_diff, y_diff, z_diff * z_diff))
     }
 }
 
diff --git a/src/direction.rs b/src/direction.rs
index 837329a..38d77e5 100644
--- a/src/direction.rs
+++ b/src/direction.rs
@@ -2,7 +2,7 @@ use crate::{
     coord::CoordIJK, error, CellIndex, Edge, Vertex, NUM_HEX_VERTS,
     NUM_PENT_VERTS,
 };
-use std::{fmt, num::NonZeroU8};
+use core::{fmt, num::NonZeroU8};
 
 /// Maximum value for a direction.
 const MAX: u8 = 6;
@@ -111,7 +111,7 @@ impl Direction {
     pub(crate) const fn new_unchecked(value: u8) -> Self {
         assert!(value <= MAX, "direction out of range");
         // SAFETY: range checked above.
-        unsafe { std::mem::transmute::<u8, Self>(value) }
+        unsafe { core::mem::transmute::<u8, Self>(value) }
     }
 
     /// Returns a direction rotated `count` time, by 60 degrees step.
diff --git a/src/error/compaction.rs b/src/error/compaction.rs
index 294a644..133eedc 100644
--- a/src/error/compaction.rs
+++ b/src/error/compaction.rs
@@ -1,4 +1,4 @@
-use std::{error::Error, fmt};
+use core::fmt;
 
 /// Errors occurring while compacting a set of cell indices.
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
@@ -21,8 +21,9 @@ impl fmt::Display for CompactionError {
     }
 }
 
-impl Error for CompactionError {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
+#[cfg(feature = "std")]
+impl std::error::Error for CompactionError {
+    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
         None
     }
 }
diff --git a/src/error/geometry.rs b/src/error/geometry.rs
index a54188b..96957df 100644
--- a/src/error/geometry.rs
+++ b/src/error/geometry.rs
@@ -1,4 +1,4 @@
-use std::{error::Error, fmt};
+use core::fmt;
 
 /// Errors related to the geometries.
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
@@ -19,8 +19,9 @@ impl fmt::Display for InvalidGeometry {
     }
 }
 
-impl Error for InvalidGeometry {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
+#[cfg(feature = "std")]
+impl std::error::Error for InvalidGeometry {
+    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
         None
     }
 }
diff --git a/src/error/hex_grid.rs b/src/error/hex_grid.rs
index e088707..102fa2e 100644
--- a/src/error/hex_grid.rs
+++ b/src/error/hex_grid.rs
@@ -1,4 +1,4 @@
-use std::{error::Error, fmt};
+use core::fmt;
 
 /// Errors related to the `IJK` coordinate system and its variants (e.g.
 /// [`LocalIJ`](crate::LocalIJ)).
@@ -20,8 +20,9 @@ impl fmt::Display for HexGridError {
     }
 }
 
-impl Error for HexGridError {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
+#[cfg(feature = "std")]
+impl std::error::Error for HexGridError {
+    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
         None
     }
 }
diff --git a/src/error/invalid_value.rs b/src/error/invalid_value.rs
index cc958a9..9b9e727 100644
--- a/src/error/invalid_value.rs
+++ b/src/error/invalid_value.rs
@@ -1,4 +1,4 @@
-use std::{error::Error, fmt};
+use core::fmt;
 
 // Macro to declare type-specific InvalidValue error type.
 macro_rules! invalid_value_error {
@@ -33,8 +33,9 @@ macro_rules! invalid_value_error {
             }
         }
 
-        impl Error for $error {
-            fn source(&self) -> Option<&(dyn Error + 'static)> {
+        #[cfg(feature = "std")]
+        impl std::error::Error for $error {
+            fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
                 None
             }
         }
diff --git a/src/error/localij.rs b/src/error/localij.rs
index 6a57f60..60aa2a0 100644
--- a/src/error/localij.rs
+++ b/src/error/localij.rs
@@ -1,5 +1,5 @@
 use super::HexGridError;
-use std::{error::Error, fmt};
+use core::fmt;
 
 /// Errors occurring during [`LocalIJ`](crate::LocalIJ) coordinate system
 /// conversions.
@@ -26,8 +26,9 @@ impl fmt::Display for LocalIjError {
     }
 }
 
-impl Error for LocalIjError {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
+#[cfg(feature = "std")]
+impl std::error::Error for LocalIjError {
+    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
         match *self {
             Self::ResolutionMismatch | Self::Pentagon => None,
             Self::HexGrid(ref err) => Some(err),
diff --git a/src/error/outliner.rs b/src/error/outliner.rs
index 7b79bda..317869f 100644
--- a/src/error/outliner.rs
+++ b/src/error/outliner.rs
@@ -1,4 +1,4 @@
-use std::{error::Error, fmt};
+use core::fmt;
 
 /// Errors occurring during the outline computation of a set of cell indices.
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
@@ -21,8 +21,9 @@ impl fmt::Display for OutlinerError {
     }
 }
 
-impl Error for OutlinerError {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
+#[cfg(feature = "std")]
+impl std::error::Error for OutlinerError {
+    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
         None
     }
 }
diff --git a/src/error/resolution_mismatch.rs b/src/error/resolution_mismatch.rs
index 8441a0a..2ba30e9 100644
--- a/src/error/resolution_mismatch.rs
+++ b/src/error/resolution_mismatch.rs
@@ -1,4 +1,4 @@
-use std::{error::Error, fmt};
+use core::fmt;
 
 /// Resolution mismatch between two cell indexes.
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
@@ -10,8 +10,9 @@ impl fmt::Display for ResolutionMismatch {
     }
 }
 
-impl Error for ResolutionMismatch {
-    fn source(&self) -> Option<&(dyn Error + 'static)> {
+#[cfg(feature = "std")]
+impl std::error::Error for ResolutionMismatch {
+    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
         None
     }
 }
diff --git a/src/error/tests.rs b/src/error/tests.rs
index 3bb9f40..50df4c1 100644
--- a/src/error/tests.rs
+++ b/src/error/tests.rs
@@ -6,7 +6,7 @@ use crate::error::{
 };
 #[cfg(feature = "geo")]
 use crate::error::{InvalidGeometry, OutlinerError};
-use std::error::Error;
+use alloc::string::ToString;
 
 // All error must have a non-empty display.
 #[test]
@@ -57,8 +57,11 @@ fn display() {
 }
 
 // All errors are root errors.
+#[cfg(feature = "std")]
 #[test]
 fn source() {
+    use std::error::Error;
+
     let hex_grid_error = HexGridError::new("error");
 
     assert!(CompactionError::HeterogeneousResolution.source().is_none());
diff --git a/src/face.rs b/src/face.rs
index 8395da7..7580b8b 100644
--- a/src/face.rs
+++ b/src/face.rs
@@ -4,7 +4,7 @@ use crate::{
     coord::{CoordIJK, LatLng, Vec3d},
     error, NUM_ICOSA_FACES,
 };
-use std::fmt;
+use core::fmt;
 
 // -----------------------------------------------------------------------------
 
@@ -143,14 +143,14 @@ impl FaceSet {
 
 impl fmt::Display for FaceSet {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        write!(
-            f,
-            "[{}]",
-            self.iter()
-                .map(|face| face.to_string())
-                .collect::<Vec<_>>()
-                .join("-")
-        )
+        write!(f, "[",)?;
+        for (i, face) in self.iter().enumerate() {
+            if i != 0 {
+                write!(f, "-")?;
+            }
+            write!(f, "{face}")?;
+        }
+        write!(f, "]",)
     }
 }
 
diff --git a/src/geom/geometry/bbox.rs b/src/geom/geometry/bbox.rs
index e0b76a3..f92b478 100644
--- a/src/geom/geometry/bbox.rs
+++ b/src/geom/geometry/bbox.rs
@@ -96,7 +96,7 @@ pub fn hex_estimate(bbox: &Rect, resolution: Resolution) -> usize {
     #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
     let estimate = estimate as usize;
 
-    std::cmp::max(estimate, 1)
+    core::cmp::max(estimate, 1)
 }
 
 fn get_min_max(value: f64, min: f64, max: f64) -> (f64, f64) {
diff --git a/src/geom/geometry/geometrycollection.rs b/src/geom/geometry/geometrycollection.rs
index c70a96d..1994540 100644
--- a/src/geom/geometry/geometrycollection.rs
+++ b/src/geom/geometry/geometrycollection.rs
@@ -4,7 +4,7 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex,
 };
-use std::boxed::Box;
+use alloc::{boxed::Box, vec::Vec};
 
 /// A collection of [`geo::Geometry`].
 #[derive(Clone, Debug, PartialEq)]
diff --git a/src/geom/geometry/line.rs b/src/geom/geometry/line.rs
index 1c29c4a..1240cee 100644
--- a/src/geom/geometry/line.rs
+++ b/src/geom/geometry/line.rs
@@ -3,8 +3,8 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex, LatLng, Resolution,
 };
+use alloc::{boxed::Box, vec::Vec};
 use geo::Coord;
-use std::boxed::Box;
 
 /// A line segment made up of exactly two [`geo::Coord`]s.
 ///
diff --git a/src/geom/geometry/linestring.rs b/src/geom/geometry/linestring.rs
index 3203eff..6d0eb73 100644
--- a/src/geom/geometry/linestring.rs
+++ b/src/geom/geometry/linestring.rs
@@ -4,7 +4,7 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex,
 };
-use std::boxed::Box;
+use alloc::boxed::Box;
 
 /// An ordered collection of two or more [`geo::Coord`]s, representing a
 /// path between locations.
diff --git a/src/geom/geometry/mod.rs b/src/geom/geometry/mod.rs
index 202de60..11dadb6 100644
--- a/src/geom/geometry/mod.rs
+++ b/src/geom/geometry/mod.rs
@@ -3,7 +3,8 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex, TWO_PI,
 };
-use std::{boxed::Box, f64::consts::PI};
+use alloc::boxed::Box;
+use core::f64::consts::PI;
 
 mod bbox;
 mod geometrycollection;
diff --git a/src/geom/geometry/multilinestring.rs b/src/geom/geometry/multilinestring.rs
index bc0fb06..62af34f 100644
--- a/src/geom/geometry/multilinestring.rs
+++ b/src/geom/geometry/multilinestring.rs
@@ -4,7 +4,7 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex,
 };
-use std::boxed::Box;
+use alloc::{boxed::Box, vec::Vec};
 
 /// A collection of [`geo::LineString`].
 ///
diff --git a/src/geom/geometry/multipoint.rs b/src/geom/geometry/multipoint.rs
index ee23703..3a6905e 100644
--- a/src/geom/geometry/multipoint.rs
+++ b/src/geom/geometry/multipoint.rs
@@ -4,7 +4,7 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex,
 };
-use std::boxed::Box;
+use alloc::{boxed::Box, vec::Vec};
 
 /// A collection of [`geo::Point`]s.
 #[derive(Clone, Debug, PartialEq)]
diff --git a/src/geom/geometry/multipolygon.rs b/src/geom/geometry/multipolygon.rs
index ca5e469..d0332db 100644
--- a/src/geom/geometry/multipolygon.rs
+++ b/src/geom/geometry/multipolygon.rs
@@ -4,7 +4,7 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex,
 };
-use std::boxed::Box;
+use alloc::{boxed::Box, vec::Vec};
 
 /// A collection of [`geo::Polygon`].
 #[derive(Clone, Debug, PartialEq)]
diff --git a/src/geom/geometry/point.rs b/src/geom/geometry/point.rs
index 5291575..0018ecd 100644
--- a/src/geom/geometry/point.rs
+++ b/src/geom/geometry/point.rs
@@ -3,7 +3,7 @@ use crate::{
     geom::{PolyfillConfig, ToCells},
     CellIndex, LatLng,
 };
-use std::boxed::Box;
+use alloc::boxed::Box;
 
 /// A single point in 2D space.
 #[derive(Clone, Copy, Debug, PartialEq)]
@@ -89,6 +89,6 @@ impl ToCells for Point {
         config: PolyfillConfig,
     ) -> Box<dyn Iterator<Item = CellIndex> + '_> {
         let ll = LatLng::try_from(*self).expect("valid coordinate");
-        Box::new(std::iter::once(ll.to_cell(config.resolution)))
+        Box::new(core::iter::once(ll.to_cell(config.resolution)))
     }
 }
diff --git a/src/geom/geometry/polygon.rs b/src/geom/geometry/polygon.rs
index e10db34..e520c5e 100644
--- a/src/geom/geometry/polygon.rs
+++ b/src/geom/geometry/polygon.rs
@@ -4,12 +4,22 @@ use crate::{
     geom::{ContainmentMode, PolyfillConfig, ToCells},
     CellIndex, LatLng, Resolution,
 };
-use ahash::{HashSet, HashSetExt};
+use alloc::{borrow::Cow, boxed::Box, vec::Vec};
+use core::cmp;
 use geo::{coord, CoordsIter};
-use std::{borrow::Cow, boxed::Box, cmp};
+
+#[cfg(feature = "std")]
+use ahash::{HashSet, HashSetExt};
+#[cfg(not(feature = "std"))]
+use alloc::collections::BTreeSet;
 
 type ContainmentPredicate = fn(polygon: &Polygon, cell: CellIndex) -> bool;
 
+#[cfg(not(feature = "std"))]
+type Set<K> = BTreeSet<K>;
+#[cfg(feature = "std")]
+type Set<K> = HashSet<K>;
+
 /// A bounded two-dimensional area.
 #[derive(Clone, Debug, PartialEq)]
 pub struct Polygon {
@@ -147,7 +157,7 @@ impl Polygon {
     fn hex_outline(
         &self,
         resolution: Resolution,
-        already_seen: &mut HashSet<CellIndex>,
+        already_seen: &mut Set<CellIndex>,
         scratchpad: &mut [u64],
         contains: ContainmentPredicate,
     ) -> Vec<CellIndex> {
@@ -157,7 +167,7 @@ impl Polygon {
         // Compute the set of cells making the outlines of the polygon.
         let outlines = self
             .interiors()
-            .chain(std::iter::once(self.exterior()))
+            .chain(core::iter::once(self.exterior()))
             .flat_map(|ring| get_edge_cells(ring, resolution))
             .filter(|cell| already_seen.insert(*cell))
             .collect::<Vec<_>>();
@@ -192,7 +202,7 @@ impl Polygon {
     fn outermost_inner_cells(
         &self,
         outlines: &[CellIndex],
-        already_seen: &mut HashSet<CellIndex>,
+        already_seen: &mut Set<CellIndex>,
         scratchpad: &mut [u64],
         contains: ContainmentPredicate,
     ) -> Vec<CellIndex> {
@@ -276,7 +286,7 @@ impl ToCells for Polygon {
         };
 
         // Set used for dedup.
-        let mut seen = HashSet::new();
+        let mut seen = Set::new();
         // Scratchpad memory to store a cell and its immediate neighbors.
         // Cell itself + at most 6 neighbors = 7.
         let mut scratchpad = [0; 7];
@@ -291,7 +301,7 @@ impl ToCells for Polygon {
 
         if outlines.is_empty() && config.containment == ContainmentMode::Covers
         {
-            return Box::new(std::iter::once(
+            return Box::new(core::iter::once(
                 self.exterior.centroid().to_cell(config.resolution),
             ));
         }
@@ -305,7 +315,10 @@ impl ToCells for Polygon {
             contains,
         );
         let mut next_gen = Vec::with_capacity(candidates.len() * 7);
-        let mut new_seen = HashSet::with_capacity(seen.len());
+        #[cfg(not(feature = "std"))]
+        let mut new_seen = Set::new();
+        #[cfg(feature = "std")]
+        let mut new_seen = Set::with_capacity(seen.len());
 
         if config.containment == ContainmentMode::ContainsBoundary {
             outlines.retain(|&cell| !intersects_boundary(self, cell));
@@ -313,7 +326,7 @@ impl ToCells for Polygon {
         }
 
         // Last step: inward propagation from the outermost layers.
-        let inward_propagation = std::iter::from_fn(move || {
+        let inward_propagation = core::iter::from_fn(move || {
             if candidates.is_empty() {
                 return None;
             }
@@ -337,10 +350,10 @@ impl ToCells for Polygon {
 
             let curr_gen = candidates.clone();
 
-            std::mem::swap(&mut next_gen, &mut candidates);
+            core::mem::swap(&mut next_gen, &mut candidates);
             next_gen.clear();
 
-            std::mem::swap(&mut new_seen, &mut seen);
+            core::mem::swap(&mut new_seen, &mut seen);
             new_seen.clear();
 
             Some(curr_gen.into_iter())
diff --git a/src/geom/geometry/rect.rs b/src/geom/geometry/rect.rs
index d81e332..e1003eb 100644
--- a/src/geom/geometry/rect.rs
+++ b/src/geom/geometry/rect.rs
@@ -3,7 +3,7 @@ use crate::{
     geom::{PolyfillConfig, Polygon, ToCells},
     CellIndex,
 };
-use std::boxed::Box;
+use alloc::boxed::Box;
 
 /// An axis-aligned bounded 2D rectangle whose area is defined by minimum and
 /// maximum [`geo::Coord`]s.
diff --git a/src/geom/geometry/ring.rs b/src/geom/geometry/ring.rs
index 895cc4c..41ae65c 100644
--- a/src/geom/geometry/ring.rs
+++ b/src/geom/geometry/ring.rs
@@ -1,5 +1,7 @@
 use super::bbox;
 use crate::{error::InvalidGeometry, LatLng, TWO_PI};
+use alloc::{borrow::Cow, vec::Vec};
+use core::f64::consts::PI;
 use geo::{
     algorithm::{
         centroid::Centroid,
@@ -7,7 +9,6 @@ use geo::{
     },
     Contains, Coord, Intersects, Polygon,
 };
-use std::{borrow::Cow, f64::consts::PI};
 
 /// A closed ring and its bounding box.
 #[derive(Clone, Debug, PartialEq)]
diff --git a/src/geom/geometry/triangle.rs b/src/geom/geometry/triangle.rs
index 4deef7f..162185c 100644
--- a/src/geom/geometry/triangle.rs
+++ b/src/geom/geometry/triangle.rs
@@ -3,8 +3,8 @@ use crate::{
     geom::{PolyfillConfig, Polygon, ToCells},
     CellIndex,
 };
+use alloc::boxed::Box;
 use geo::CoordsIter;
-use std::boxed::Box;
 
 /// A bounded 2D area whose three vertices are defined by [`geo::Coord`]s.
 #[derive(Clone, Debug, PartialEq)]
diff --git a/src/geom/json.rs b/src/geom/json.rs
index b7fc613..37a3e44 100644
--- a/src/geom/json.rs
+++ b/src/geom/json.rs
@@ -1,5 +1,6 @@
 use super::{Geometry, GeometryCollection};
 use crate::error::InvalidGeometry;
+use alloc::vec::Vec;
 
 impl TryFrom<&geojson::Geometry> for Geometry {
     type Error = InvalidGeometry;
diff --git a/src/geom/ring_hierarchy.rs b/src/geom/ring_hierarchy.rs
index 4bdfd9b..fda8995 100644
--- a/src/geom/ring_hierarchy.rs
+++ b/src/geom/ring_hierarchy.rs
@@ -1,5 +1,6 @@
+use alloc::{vec, vec::Vec};
+use core::iter::Peekable;
 use geo::{Contains, Coord, LineString, MultiPolygon, Polygon};
-use std::iter::Peekable;
 
 /// A rings hierarchy.
 pub struct RingHierarchy {
@@ -86,12 +87,12 @@ impl RingHierarchy {
     /// Consumes the hierarchy into a stream of Polygon.
     pub fn into_iter(mut self) -> impl Iterator<Item = Polygon<f64>> {
         type OuterRingIterator =
-            Peekable<std::vec::IntoIter<(usize, LineString<f64>)>>;
+            Peekable<alloc::vec::IntoIter<(usize, LineString<f64>)>>;
 
         // Outers ring at the current nesting level.
         let mut outers: Option<OuterRingIterator> = None;
 
-        std::iter::from_fn(move || {
+        core::iter::from_fn(move || {
             // If the current layer is exhausted, peel the next one.
             if outers.as_mut().map_or(true, |rings| rings.peek().is_none()) {
                 outers = self
@@ -125,7 +126,7 @@ impl RingHierarchy {
                 (!self.is_assigned[i] && self.is_outer(i)).then(|| {
                     // Extract the ring in place to preserve `rings` size and
                     // ordering.
-                    let ring = std::mem::replace(
+                    let ring = core::mem::replace(
                         &mut self.rings[i],
                         LineString(vec![]),
                     );
@@ -152,7 +153,7 @@ impl RingHierarchy {
                 .then(|| {
                     // Extract the ring in place to preserve `rings` size and
                     // ordering.
-                    let ring = std::mem::replace(
+                    let ring = core::mem::replace(
                         &mut self.rings[inner_id],
                         LineString(vec![]),
                     );
diff --git a/src/geom/to_geo.rs b/src/geom/to_geo.rs
index a71f4d3..0f7b2ef 100644
--- a/src/geom/to_geo.rs
+++ b/src/geom/to_geo.rs
@@ -2,8 +2,9 @@ use super::VertexGraph;
 use crate::{
     error::OutlinerError, CellIndex, DirectedEdgeIndex, LatLng, VertexIndex,
 };
+use alloc::vec::Vec;
+use core::convert::Infallible;
 use geo::{Coord, Line, LineString, MultiPolygon, Point, Polygon};
-use std::convert::Infallible;
 
 /// A trait to trace the outline of an H3 object.
 pub trait ToGeo
diff --git a/src/geom/to_h3.rs b/src/geom/to_h3.rs
index 61eb0ae..1fe2164 100644
--- a/src/geom/to_h3.rs
+++ b/src/geom/to_h3.rs
@@ -1,5 +1,5 @@
 use crate::{CellIndex, Resolution};
-use std::boxed::Box;
+use alloc::boxed::Box;
 
 /// A trait to convert a geometry (or a collection of geometries) into a list of
 /// cell indexes of the specified resolution.
diff --git a/src/geom/vertex_graph.rs b/src/geom/vertex_graph.rs
index d524c4d..624c8fd 100644
--- a/src/geom/vertex_graph.rs
+++ b/src/geom/vertex_graph.rs
@@ -1,11 +1,25 @@
 use super::RingHierarchy;
 use crate::{error::OutlinerError, CellIndex, LatLng, Resolution, VertexIndex};
-use ahash::{HashMap, HashMapExt};
+use alloc::{vec, vec::Vec};
 use geo::{LineString, MultiPolygon, Polygon};
-use std::collections::hash_map::Entry;
+
+#[cfg(not(feature = "std"))]
+use alloc::collections::{btree_map::Entry, BTreeMap};
+#[cfg(feature = "std")]
+use {
+    ahash::{HashMap, HashMapExt},
+    std::collections::hash_map::Entry,
+};
+
+#[cfg(not(feature = "std"))]
+type Map<K, V> = BTreeMap<K, V>;
+#[cfg(feature = "std")]
+type Map<K, V> = HashMap<K, V>;
 
 /// A single node in a vertex graph.
-#[derive(Debug, Eq, PartialEq, Hash)]
+#[derive(Debug, Eq, PartialEq)]
+#[cfg_attr(feature = "std", derive(Hash))]
+#[cfg_attr(not(feature = "std"), derive(Ord, PartialOrd))]
 pub struct Node {
     from: VertexIndex,
     to: VertexIndex,
@@ -14,8 +28,8 @@ pub struct Node {
 /// A data structure to store a graph of vertices.
 #[derive(Default)]
 pub struct VertexGraph {
-    nodes: HashMap<VertexIndex, Vec<VertexIndex>>,
-    distortions: HashMap<Node, LatLng>,
+    nodes: Map<VertexIndex, Vec<VertexIndex>>,
+    distortions: Map<Node, LatLng>,
     is_class3: bool,
 }
 
@@ -41,8 +55,8 @@ impl VertexGraph {
             .copied()
             .map_or_else(|| Resolution::Zero, CellIndex::resolution);
         let mut graph = Self {
-            nodes: HashMap::new(),
-            distortions: HashMap::new(),
+            nodes: Map::new(),
+            distortions: Map::new(),
             is_class3: resolution.is_class3(),
         };
 
diff --git a/src/grid/iterator.rs b/src/grid/iterator.rs
index ca8cf72..25e3956 100644
--- a/src/grid/iterator.rs
+++ b/src/grid/iterator.rs
@@ -1,6 +1,15 @@
 use crate::{CellIndex, Direction};
+use alloc::collections::VecDeque;
+
+#[cfg(feature = "std")]
 use ahash::{HashSet, HashSetExt};
-use std::collections::VecDeque;
+#[cfg(not(feature = "std"))]
+use alloc::collections::BTreeSet;
+
+#[cfg(not(feature = "std"))]
+type Set<K> = BTreeSet<K>;
+#[cfg(feature = "std")]
+type Set<K> = HashSet<K>;
 
 /// Direction to the next ring.
 const NEXT_RING_DIRECTION: Direction = Direction::I;
@@ -32,7 +41,7 @@ pub struct DiskDistancesSafe {
     k: u32,
 
     /// Already visited neighbors.
-    seen: HashSet<CellIndex>,
+    seen: Set<CellIndex>,
     /// Next set of neighbors to visit.
     candidates: VecDeque<(CellIndex, u32)>,
 }
@@ -49,7 +58,10 @@ impl DiskDistancesSafe {
         candidates.push_back((origin, 0));
         Self {
             k,
-            seen: HashSet::with_capacity(size),
+            #[cfg(feature = "std")]
+            seen: Set::with_capacity(size),
+            #[cfg(not(feature = "std"))]
+            seen: Set::new(),
             candidates,
         }
     }
diff --git a/src/grid/iterator_tests.rs b/src/grid/iterator_tests.rs
index 44d22ed..7e20b54 100644
--- a/src/grid/iterator_tests.rs
+++ b/src/grid/iterator_tests.rs
@@ -1,4 +1,4 @@
-use std::convert::TryFrom;
+use core::convert::TryFrom;
 
 use super::*;
 
diff --git a/src/index/bits.rs b/src/index/bits.rs
index 2a3b679..386f35f 100644
--- a/src/index/bits.rs
+++ b/src/index/bits.rs
@@ -2,7 +2,7 @@
 
 use super::IndexMode;
 use crate::{Direction, Edge, Resolution, Vertex};
-use std::{cmp, num::NonZeroU8};
+use core::{cmp, num::NonZeroU8};
 
 /// Offset (in bits) of the mode in an H3 index.
 const MODE_OFFSET: usize = 59;
diff --git a/src/index/cell.rs b/src/index/cell.rs
index afb0518..35f5d05 100644
--- a/src/index/cell.rs
+++ b/src/index/cell.rs
@@ -11,13 +11,13 @@ use crate::{
     FaceSet, LatLng, LocalIJ, Resolution, Vertex, VertexIndex, CCW, CW,
     DEFAULT_CELL_INDEX, EARTH_RADIUS_KM, NUM_HEX_VERTS, NUM_PENT_VERTS,
 };
-use either::Either;
-use std::{
+use core::{
     cmp::Ordering,
     fmt, iter,
     num::{NonZeroU64, NonZeroU8},
     str::FromStr,
 };
+use either::Either;
 
 /// Lookup table for number of children for hexagonal cells.
 // 7.pow(resolution_delta)
diff --git a/src/index/cell_tests.rs b/src/index/cell_tests.rs
index 3538c75..170527f 100644
--- a/src/index/cell_tests.rs
+++ b/src/index/cell_tests.rs
@@ -1,4 +1,5 @@
 use super::*;
+use alloc::{format, vec};
 
 #[test]
 fn direction_at() {
diff --git a/src/index/edge.rs b/src/index/edge.rs
index 5d616f8..70cb0de 100644
--- a/src/index/edge.rs
+++ b/src/index/edge.rs
@@ -3,7 +3,7 @@ use crate::{
     coord::FaceIJK, error, grid, Boundary, CellIndex, Direction,
     EARTH_RADIUS_KM,
 };
-use std::{cmp::Ordering, fmt, num::NonZeroU64, str::FromStr};
+use core::{cmp::Ordering, fmt, num::NonZeroU64, str::FromStr};
 
 /// Minimum value for a cell edge.
 const MIN: u8 = 1;
diff --git a/src/index/edge_tests.rs b/src/index/edge_tests.rs
index 2db6d0b..3c39f90 100644
--- a/src/index/edge_tests.rs
+++ b/src/index/edge_tests.rs
@@ -1,4 +1,5 @@
 use super::*;
+use alloc::{format, vec};
 
 #[test]
 fn zero_is_invalid() {
diff --git a/src/index/iterator.rs b/src/index/iterator.rs
index 09c629e..7db0bcf 100644
--- a/src/index/iterator.rs
+++ b/src/index/iterator.rs
@@ -5,7 +5,8 @@ use crate::{
     index::bits,
     Direction, Resolution,
 };
-use std::cmp::max;
+use alloc::vec::Vec;
+use core::cmp::max;
 
 /// Iterator over a children cell index at a given resolution.
 pub struct Children {
diff --git a/src/index/mode.rs b/src/index/mode.rs
index 5a83a18..2e12fa6 100644
--- a/src/index/mode.rs
+++ b/src/index/mode.rs
@@ -1,4 +1,4 @@
-use std::fmt;
+use core::fmt;
 
 /// H3 index modes.
 #[derive(Clone, Copy, Debug, Eq, PartialEq, Hash, Ord, PartialOrd)]
diff --git a/src/index/triangle.rs b/src/index/triangle.rs
index 9fe0893..c8b800b 100644
--- a/src/index/triangle.rs
+++ b/src/index/triangle.rs
@@ -1,4 +1,7 @@
-use crate::LatLng;
+use crate::{
+    math::{atan, sqrt, tan},
+    LatLng,
+};
 
 /// A triangle on unit sphere.
 pub struct Triangle {
@@ -43,5 +46,5 @@ fn area_from_edges(mut a: f64, mut b: f64, mut c: f64) -> f64 {
     c = (s - c) / 2.;
     s /= 2.;
 
-    4. * ((s.tan() * a.tan() * b.tan() * c.tan()).sqrt()).atan()
+    4. * atan(sqrt(tan(s) * tan(a) * tan(b) * tan(c)))
 }
diff --git a/src/index/vertex.rs b/src/index/vertex.rs
index 60b337f..6db8be1 100644
--- a/src/index/vertex.rs
+++ b/src/index/vertex.rs
@@ -3,7 +3,7 @@ use crate::{
     coord::FaceIJK, error, CellIndex, Direction, LatLng, NUM_HEX_VERTS,
     NUM_PENT_VERTS,
 };
-use std::{cmp::Ordering, fmt, num::NonZeroU64, str::FromStr};
+use core::{cmp::Ordering, fmt, num::NonZeroU64, str::FromStr};
 
 /// Maximum value for a cell vertex.
 const MAX: u8 = 5;
diff --git a/src/index/vertex_tests.rs b/src/index/vertex_tests.rs
index cb90155..e7ef1c3 100644
--- a/src/index/vertex_tests.rs
+++ b/src/index/vertex_tests.rs
@@ -1,4 +1,5 @@
 use super::*;
+use alloc::{format, vec};
 
 #[test]
 fn zero_is_invalid() {
diff --git a/src/lib.rs b/src/lib.rs
index c43a286..1d3a3cc 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,9 +1,24 @@
+#![cfg_attr(not(feature = "std"), no_std)]
 //! The `h3o` library implements the H3 geospatial indexing system.
 //!
 //! H3 is a geospatial indexing system using a hexagonal grid that can be
 //! (approximately) subdivided into finer and finer hexagonal grids, combining
 //! the benefits of a hexagonal grid with S2's hierarchical subdivisions.
 //!
+//! ## Crate features
+//!
+//! * **std** -
+//! When enabled, this will cause `h3o` to use the standard library. In terms of
+//! APIs, `std` causes error types to implement the `std::error::Error` trait.
+//! Enabling `std` will also result in performance optimizations.
+//!
+//! * **geo** -
+//! When enabled, you'll be able to convert lists of H3 cell indexes from and
+//! into geometric shapes. Also enables the `GeoJSON` support. Requires `std`.
+//!
+//! * **serde** -
+//! When enabled, H3 index types (cell, vertex and edge) derive serde traits.
+//!
 //! ## H3 to H3O mapping
 //!
 //! For people used to the H3 API, here is the mapping to H3O.
@@ -119,7 +134,6 @@
 //! | `greatCircleDistanceRads`   | [`LatLng::distance_rads`]          |
 
 // Lints {{{
-
 #![deny(
     nonstandard_style,
     rust_2018_idioms,
@@ -230,9 +244,10 @@
     // Too many irrelevant warning (about internal invariants).
     clippy::missing_panics_doc,
 )]
-
 // }}}
 
+extern crate alloc;
+
 use konst::{primitive::parse_u8 as as_u8, result::unwrap_ctx as unwrap};
 
 mod base_cell;
@@ -247,6 +262,13 @@ mod grid;
 mod index;
 mod resolution;
 
+#[cfg(not(feature = "std"))]
+#[path = "math-libm.rs"]
+mod math;
+#[cfg(feature = "std")]
+#[path = "math-std.rs"]
+mod math;
+
 pub use base_cell::BaseCell;
 pub use boundary::Boundary;
 pub use coord::{CoordIJ, LatLng, LocalIJ};
@@ -290,7 +312,7 @@ const NUM_PENTAGONS: u8 = 12;
 const DEFAULT_CELL_INDEX: u64 = 0x0800_1fff_ffff_ffff;
 
 // 2π
-const TWO_PI: f64 = 2. * std::f64::consts::PI;
+const TWO_PI: f64 = 2. * core::f64::consts::PI;
 
 // -----------------------------------------------------------------------------
 
diff --git a/src/math-libm.rs b/src/math-libm.rs
new file mode 100644
index 0000000..b7d7107
--- /dev/null
+++ b/src/math-libm.rs
@@ -0,0 +1,59 @@
+#[inline]
+pub fn abs(x: f64) -> f64 {
+    libm::fabs(x)
+}
+
+#[inline]
+pub fn sin(x: f64) -> f64 {
+    libm::sin(x)
+}
+
+#[inline]
+pub fn cos(x: f64) -> f64 {
+    libm::cos(x)
+}
+
+#[inline]
+pub fn tan(x: f64) -> f64 {
+    libm::tan(x)
+}
+
+#[inline]
+pub fn asin(x: f64) -> f64 {
+    libm::asin(x)
+}
+
+#[inline]
+pub fn acos(x: f64) -> f64 {
+    libm::acos(x)
+}
+
+#[inline]
+pub fn atan(x: f64) -> f64 {
+    libm::atan(x)
+}
+
+#[inline]
+pub fn atan2(y: f64, x: f64) -> f64 {
+    libm::atan2(y, x)
+}
+
+#[inline]
+pub fn hypot(x: f64, y: f64) -> f64 {
+    libm::hypot(x, y)
+}
+
+#[inline]
+pub fn sqrt(x: f64) -> f64 {
+    libm::sqrt(x)
+}
+
+#[inline]
+pub fn round(x: f64) -> f64 {
+    libm::round(x)
+}
+
+#[inline]
+pub fn mul_add(a: f64, b: f64, c: f64) -> f64 {
+    (a * b) + c
+}
diff --git a/src/math-std.rs b/src/math-std.rs
new file mode 100644
index 0000000..94e2b23
--- /dev/null
+++ b/src/math-std.rs
@@ -0,0 +1,59 @@
+#[inline]
+pub fn abs(x: f64) -> f64 {
+    x.abs()
+}
+
+#[inline]
+pub fn sin(x: f64) -> f64 {
+    x.sin()
+}
+
+#[inline]
+pub fn cos(x: f64) -> f64 {
+    x.cos()
+}
+
+#[inline]
+pub fn tan(x: f64) -> f64 {
+    x.tan()
+}
+
+#[inline]
+pub fn asin(x: f64) -> f64 {
+    x.asin()
+}
+
+#[inline]
+pub fn acos(x: f64) -> f64 {
+    x.acos()
+}
+
+#[inline]
+pub fn atan(x: f64) -> f64 {
+    x.atan()
+}
+
+#[inline]
+pub fn atan2(y: f64, x: f64) -> f64 {
+    y.atan2(x)
+}
+
+#[inline]
+pub fn hypot(x: f64, y: f64) -> f64 {
+    x.hypot(y)
+}
+
+#[inline]
+pub fn sqrt(x: f64) -> f64 {
+    x.sqrt()
+}
+
+#[inline]
+pub fn round(x: f64) -> f64 {
+    x.round()
+}
+
+#[inline]
+pub fn mul_add(a: f64, b: f64, c: f64) -> f64 {
+    a.mul_add(b, c)
+}
diff --git a/src/resolution.rs b/src/resolution.rs
index 4132f8a..38cdb25 100644
--- a/src/resolution.rs
+++ b/src/resolution.rs
@@ -1,5 +1,5 @@
 use crate::{error, index::bits, BaseCell, CellIndex, NUM_PENTAGONS};
-use std::{ffi::c_int, fmt, iter::DoubleEndedIterator, str::FromStr};
+use core::{ffi::c_int, fmt, iter::DoubleEndedIterator, str::FromStr};
 
 /// Cell resolution, from 0 to 15.
 #[derive(Clone, Copy, Debug, Eq, PartialEq, Hash, Ord, PartialOrd)]
@@ -132,7 +132,7 @@ impl Resolution {
     ) -> impl Iterator<Item = Self> + DoubleEndedIterator {
         (u8::from(start)..=u8::from(end))
             // SAFETY: values between two resolutions are valid resolutions.
-            .map(|value| unsafe { std::mem::transmute::<u8, Self>(value) })
+            .map(|value| unsafe { core::mem::transmute::<u8, Self>(value) })
     }
 
     /// Returns the average hexagon area, in square radians, at this
@@ -390,7 +390,7 @@ impl Resolution {
     pub(crate) const fn new_unchecked(value: u8) -> Self {
         assert!(value <= h3o_bit::MAX_RESOLUTION, "resolution out of range");
         // SAFETY: range is checked above!
-        unsafe { std::mem::transmute::<u8, Self>(value) }
+        unsafe { core::mem::transmute::<u8, Self>(value) }
     }
 
     /// Returns the bit offset of the direction at this resolution in an H3
diff --git a/tests/api/geom/to_cells.rs b/tests/api/geom/to_cells.rs
index 4013f89..72e9f99 100644
--- a/tests/api/geom/to_cells.rs
+++ b/tests/api/geom/to_cells.rs
@@ -1,4 +1,4 @@
-use ahash::HashSet;
+use alloc::collections::BTreeSet;
 use geo::{coord, polygon, LineString};
 use h3o::{
     geom::{ContainmentMode, PolyfillConfig, Polygon, ToCells},
@@ -30,7 +30,7 @@ macro_rules! world_test {
             ];
             let poly1 = Polygon::from_radians(shape1).expect("poly 1");
             let count1 = poly1.max_cells_count(config);
-            let cells1 = poly1.to_cells(config).collect::<HashSet<_>>();
+            let cells1 = poly1.to_cells(config).collect::<BTreeSet<_>>();
 
             assert_eq!(count1, $expected);
             assert!(count1 >= cells1.len());
@@ -44,7 +44,7 @@ macro_rules! world_test {
             ];
             let poly2 = Polygon::from_radians(shape2).expect("poly 2");
             let count2 = poly2.max_cells_count(config);
-            let cells2 = poly2.to_cells(config).collect::<HashSet<_>>();
+            let cells2 = poly2.to_cells(config).collect::<BTreeSet<_>>();
 
             assert_eq!(count2, $expected);
             assert!(count2 >= cells2.len());
@@ -258,9 +258,9 @@ macro_rules! exhaustive_test {
                 let shape = geo::Polygon::new(ring, Vec::new());
                 let polygon = Polygon::from_radians(shape).expect("polygon");
 
-                let result = polygon.to_cells(config).collect::<HashSet<_>>();
+                let result = polygon.to_cells(config).collect::<BTreeSet<_>>();
                 let expected =
-                    index.children(resolution).collect::<HashSet<_>>();
+                    index.children(resolution).collect::<BTreeSet<_>>();
                 assert_eq!(
                     result, expected,
                     "cell {index} at given resolution"
@@ -269,8 +269,9 @@ macro_rules! exhaustive_test {
                 let next_res = Resolution::try_from($resolution + 1)
                     .expect("next resolution");
                 let config = PolyfillConfig::new(next_res);
-                let result = polygon.to_cells(config).collect::<HashSet<_>>();
-                let expected = index.children(next_res).collect::<HashSet<_>>();
+                let result = polygon.to_cells(config).collect::<BTreeSet<_>>();
+                let expected =
+                    index.children(next_res).collect::<BTreeSet<_>>();
                 assert_eq!(result, expected, "cell {index} at next resolution");
             }
         }
diff --git a/tests/integration_tests.rs b/tests/integration_tests.rs
index 387fa9f..d67b57d 100644
--- a/tests/integration_tests.rs
+++ b/tests/integration_tests.rs
@@ -1,3 +1,5 @@
+extern crate alloc;
+
 mod api;
 // Test against the reference implementation.
 mod h3;