Merge branch 'master' into rebased

70 files changed, 2083 insertions, 773 deletions

diff --git a/std/math/acos.zig b/std/math/acos.zig
index 763d9d8abd..de07da8fe0 100644
--- a/std/math/acos.zig
+++ b/std/math/acos.zig
@@ -1,11 +1,17 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - acos(x)   = nan if x < -1 or x > 1
+// https://git.musl-libc.org/cgit/musl/tree/src/math/acosf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/acos.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the arc-cosine of x.
+///
+/// Special cases:
+///  - acos(x)   = nan if x < -1 or x > 1
 pub fn acos(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/acosh.zig b/std/math/acosh.zig
index a2a9926863..503c0433fc 100644
--- a/std/math/acosh.zig
+++ b/std/math/acosh.zig
@@ -1,13 +1,19 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - acosh(x)   = snan if x < 1
-// - acosh(nan) = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/acoshf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/acosh.c
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the hyperbolic arc-cosine of x.
+///
+/// Special cases:
+///  - acosh(x)   = snan if x < 1
+///  - acosh(nan) = nan
 pub fn acosh(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/asin.zig b/std/math/asin.zig
index b01a1ac49e..2db9f86ff1 100644
--- a/std/math/asin.zig
+++ b/std/math/asin.zig
@@ -1,12 +1,18 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - asin(+-0) = +-0
-// - asin(x)   = nan if x < -1 or x > 1
+// https://git.musl-libc.org/cgit/musl/tree/src/math/asinf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/asin.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the arc-sin of x.
+///
+/// Special Cases:
+///  - asin(+-0) = +-0
+///  - asin(x)   = nan if x < -1 or x > 1
 pub fn asin(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/asinh.zig b/std/math/asinh.zig
index cef9bfb23f..0fb51d1b43 100644
--- a/std/math/asinh.zig
+++ b/std/math/asinh.zig
@@ -1,14 +1,20 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - asinh(+-0)   = +-0
-// - asinh(+-inf) = +-inf
-// - asinh(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/asinhf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/asinh.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns the hyperbolic arc-sin of x.
+///
+/// Special Cases:
+///  - asinh(+-0)   = +-0
+///  - asinh(+-inf) = +-inf
+///  - asinh(nan)   = nan
 pub fn asinh(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/atan.zig b/std/math/atan.zig
index ba5a11dd10..cc4cad0dd4 100644
--- a/std/math/atan.zig
+++ b/std/math/atan.zig
@@ -1,12 +1,18 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - atan(+-0)   = +-0
-// - atan(+-inf) = +-pi/2
+// https://git.musl-libc.org/cgit/musl/tree/src/math/atanf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/atan.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the arc-tangent of x.
+///
+/// Special Cases:
+///  - atan(+-0)   = +-0
+///  - atan(+-inf) = +-pi/2
 pub fn atan(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/atan2.zig b/std/math/atan2.zig
index 7f13507402..68e381607d 100644
--- a/std/math/atan2.zig
+++ b/std/math/atan2.zig
@@ -1,27 +1,33 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-//  atan2(y, nan)     = nan
-//  atan2(nan, x)     = nan
-//  atan2(+0, x>=0)   = +0
-//  atan2(-0, x>=0)   = -0
-//  atan2(+0, x<=-0)  = +pi
-//  atan2(-0, x<=-0)  = -pi
-//  atan2(y>0, 0)     = +pi/2
-//  atan2(y<0, 0)     = -pi/2
-//  atan2(+inf, +inf) = +pi/4
-//  atan2(-inf, +inf) = -pi/4
-//  atan2(+inf, -inf) = 3pi/4
-//  atan2(-inf, -inf) = -3pi/4
-//  atan2(y, +inf)    = 0
-//  atan2(y>0, -inf)  = +pi
-//  atan2(y<0, -inf)  = -pi
-//  atan2(+inf, x)    = +pi/2
-//  atan2(-inf, x)    = -pi/2
+// https://git.musl-libc.org/cgit/musl/tree/src/math/atan2f.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/atan2.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the arc-tangent of y/x.
+///
+/// Special Cases:
+///  - atan2(y, nan)     = nan
+///  - atan2(nan, x)     = nan
+///  - atan2(+0, x>=0)   = +0
+///  - atan2(-0, x>=0)   = -0
+///  - atan2(+0, x<=-0)  = +pi
+///  - atan2(-0, x<=-0)  = -pi
+///  - atan2(y>0, 0)     = +pi/2
+///  - atan2(y<0, 0)     = -pi/2
+///  - atan2(+inf, +inf) = +pi/4
+///  - atan2(-inf, +inf) = -pi/4
+///  - atan2(+inf, -inf) = 3pi/4
+///  - atan2(-inf, -inf) = -3pi/4
+///  - atan2(y, +inf)    = 0
+///  - atan2(y>0, -inf)  = +pi
+///  - atan2(y<0, -inf)  = -pi
+///  - atan2(+inf, x)    = +pi/2
+///  - atan2(-inf, x)    = -pi/2
 pub fn atan2(comptime T: type, y: T, x: T) T {
     return switch (T) {
         f32 => atan2_32(y, x),
diff --git a/std/math/atanh.zig b/std/math/atanh.zig
index 9056064f5a..8ba29be761 100644
--- a/std/math/atanh.zig
+++ b/std/math/atanh.zig
@@ -1,14 +1,20 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - atanh(+-1) = +-inf with signal
-// - atanh(x)   = nan if |x| > 1 with signal
-// - atanh(nan) = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/atanhf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/atanh.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns the hyperbolic arc-tangent of x.
+///
+/// Special Cases:
+///  - atanh(+-1) = +-inf with signal
+///  - atanh(x)   = nan if |x| > 1 with signal
+///  - atanh(nan) = nan
 pub fn atanh(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/big.zig b/std/math/big.zig
index 94b6d864e7..44b5ce675f 100644
--- a/std/math/big.zig
+++ b/std/math/big.zig
@@ -1,5 +1,7 @@
 pub use @import("big/int.zig");
+pub use @import("big/rational.zig");
 
 test "math.big" {
     _ = @import("big/int.zig");
+    _ = @import("big/rational.zig");
 }
diff --git a/std/math/big/int.zig b/std/math/big/int.zig
index 8800c2c7a9..beac3c85fe 100644
--- a/std/math/big/int.zig
+++ b/std/math/big/int.zig
@@ -21,78 +21,160 @@ comptime {
     debug.assert(Limb.is_signed == false);
 }
 
+/// An arbitrary-precision big integer.
+///
+/// Memory is allocated by an Int as needed to ensure operations never overflow. The range of an
+/// Int is bounded only by available memory.
 pub const Int = struct {
-    allocator: *Allocator,
-    positive: bool,
-    //  - little-endian ordered
-    //  - len >= 1 always
-    //  - zero value -> len == 1 with limbs[0] == 0
+    const sign_bit: usize = 1 << (usize.bit_count - 1);
+
+    /// Default number of limbs to allocate on creation of an Int.
+    pub const default_capacity = 4;
+
+    /// Allocator used by the Int when requesting memory.
+    allocator: ?*Allocator,
+
+    /// Raw digits. These are:
+    ///
+    /// * Little-endian ordered
+    /// * limbs.len >= 1
+    /// * Zero is represent as Int.len() == 1 with limbs[0] == 0.
+    ///
+    /// Accessing limbs directly should be avoided.
     limbs: []Limb,
-    len: usize,
 
-    const default_capacity = 4;
+    /// High bit is the sign bit. If set, Int is negative, else Int is positive.
+    /// The remaining bits represent the number of limbs used by Int.
+    metadata: usize,
 
+    /// Creates a new Int. default_capacity limbs will be allocated immediately.
+    /// Int will be zeroed.
     pub fn init(allocator: *Allocator) !Int {
         return try Int.initCapacity(allocator, default_capacity);
     }
 
+    /// Creates a new Int. Int will be set to `value`.
+    ///
+    /// This is identical to an `init`, followed by a `set`.
     pub fn initSet(allocator: *Allocator, value: var) !Int {
         var s = try Int.init(allocator);
         try s.set(value);
         return s;
     }
 
+    /// Creates a new Int with a specific capacity. If capacity < default_capacity then the
+    /// default capacity will be used instead.
     pub fn initCapacity(allocator: *Allocator, capacity: usize) !Int {
         return Int{
             .allocator = allocator,
-            .positive = true,
+            .metadata = 1,
             .limbs = block: {
                 var limbs = try allocator.alloc(Limb, math.max(default_capacity, capacity));
                 limbs[0] = 0;
                 break :block limbs;
             },
-            .len = 1,
         };
     }
 
+    /// Returns the number of limbs currently in use.
+    pub fn len(self: Int) usize {
+        return self.metadata & ~sign_bit;
+    }
+
+    /// Returns whether an Int is positive.
+    pub fn isPositive(self: Int) bool {
+        return self.metadata & sign_bit == 0;
+    }
+
+    /// Sets the sign of an Int.
+    pub fn setSign(self: *Int, positive: bool) void {
+        if (positive) {
+            self.metadata &= ~sign_bit;
+        } else {
+            self.metadata |= sign_bit;
+        }
+    }
+
+    /// Sets the length of an Int.
+    ///
+    /// If setLen is used, then the Int must be normalized to suit.
+    pub fn setLen(self: *Int, new_len: usize) void {
+        self.metadata &= sign_bit;
+        self.metadata |= new_len;
+    }
+
+    /// Returns an Int backed by a fixed set of limb values.
+    /// This is read-only and cannot be used as a result argument. If the Int tries to allocate
+    /// memory a runtime panic will occur.
+    pub fn initFixed(limbs: []const Limb) Int {
+        var self = Int{
+            .allocator = null,
+            .metadata = limbs.len,
+            // Cast away the const, invalid use to pass as a pointer argument.
+            .limbs = @intToPtr([*]Limb, @ptrToInt(limbs.ptr))[0..limbs.len],
+        };
+
+        self.normalize(limbs.len);
+        return self;
+    }
+
+    /// Ensures an Int has enough space allocated for capacity limbs. If the Int does not have
+    /// sufficient capacity, the exact amount will be allocated. This occurs even if the requested
+    /// capacity is only greater than the current capacity by one limb.
     pub fn ensureCapacity(self: *Int, capacity: usize) !void {
+        self.assertWritable();
         if (capacity <= self.limbs.len) {
             return;
         }
 
-        self.limbs = try self.allocator.realloc(self.limbs, capacity);
+        self.limbs = try self.allocator.?.realloc(self.limbs, capacity);
     }
 
+    fn assertWritable(self: Int) void {
+        if (self.allocator == null) {
+            @panic("provided Int value is read-only but must be writable");
+        }
+    }
+
+    /// Frees all memory associated with an Int.
     pub fn deinit(self: *Int) void {
-        self.allocator.free(self.limbs);
+        self.assertWritable();
+        self.allocator.?.free(self.limbs);
         self.* = undefined;
     }
 
+    /// Clones an Int and returns a new Int with the same value. The new Int is a deep copy and
+    /// can be modified separately from the original.
     pub fn clone(other: Int) !Int {
+        other.assertWritable();
         return Int{
             .allocator = other.allocator,
-            .positive = other.positive,
+            .metadata = other.metadata,
             .limbs = block: {
-                var limbs = try other.allocator.alloc(Limb, other.len);
-                mem.copy(Limb, limbs[0..], other.limbs[0..other.len]);
+                var limbs = try other.allocator.?.alloc(Limb, other.len());
+                mem.copy(Limb, limbs[0..], other.limbs[0..other.len()]);
                 break :block limbs;
             },
-            .len = other.len,
         };
     }
 
+    /// Copies the value of an Int to an existing Int so that they both have the same value.
+    /// Extra memory will be allocated if the receiver does not have enough capacity.
     pub fn copy(self: *Int, other: Int) !void {
-        if (self == &other) {
+        self.assertWritable();
+        if (self.limbs.ptr == other.limbs.ptr) {
             return;
         }
 
-        self.positive = other.positive;
-        try self.ensureCapacity(other.len);
-        mem.copy(Limb, self.limbs[0..], other.limbs[0..other.len]);
-        self.len = other.len;
+        try self.ensureCapacity(other.len());
+        mem.copy(Limb, self.limbs[0..], other.limbs[0..other.len()]);
+        self.metadata = other.metadata;
     }
 
+    /// Efficiently swap an Int with another. This swaps the limb pointers and a full copy is not
+    /// performed. The address of the limbs field will not be the same after this function.
     pub fn swap(self: *Int, other: *Int) void {
+        self.assertWritable();
         mem.swap(Int, self, other);
     }
 
@@ -103,45 +185,49 @@ pub const Int = struct {
         debug.warn("\n");
     }
 
-    pub fn negate(r: *Int) void {
-        r.positive = !r.positive;
+    /// Negate the sign of an Int.
+    pub fn negate(self: *Int) void {
+        self.metadata ^= sign_bit;
     }
 
-    pub fn abs(r: *Int) void {
-        r.positive = true;
+    /// Make an Int positive.
+    pub fn abs(self: *Int) void {
+        self.metadata &= ~sign_bit;
     }
 
-    pub fn isOdd(r: Int) bool {
-        return r.limbs[0] & 1 != 0;
+    /// Returns true if an Int is odd.
+    pub fn isOdd(self: Int) bool {
+        return self.limbs[0] & 1 != 0;
     }
 
-    pub fn isEven(r: Int) bool {
-        return !r.isOdd();
+    /// Returns true if an Int is even.
+    pub fn isEven(self: Int) bool {
+        return !self.isOdd();
     }
 
-    // Returns the number of bits required to represent the absolute value of self.
+    /// Returns the number of bits required to represent the absolute value an Int.
     fn bitCountAbs(self: Int) usize {
-        return (self.len - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len - 1]));
+        return (self.len() - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len() - 1]));
     }
 
-    // Returns the number of bits required to represent the integer in twos-complement form.
-    //
-    // If the integer is negative the value returned is the number of bits needed by a signed
-    // integer to represent the value. If positive the value is the number of bits for an
-    // unsigned integer. Any unsigned integer will fit in the signed integer with bitcount
-    // one greater than the returned value.
-    //
-    // e.g. -127 returns 8 as it will fit in an i8. 127 returns 7 since it fits in a u7.
+    /// Returns the number of bits required to represent the integer in twos-complement form.
+    ///
+    /// If the integer is negative the value returned is the number of bits needed by a signed
+    /// integer to represent the value. If positive the value is the number of bits for an
+    /// unsigned integer. Any unsigned integer will fit in the signed integer with bitcount
+    /// one greater than the returned value.
+    ///
+    /// e.g. -127 returns 8 as it will fit in an i8. 127 returns 7 since it fits in a u7.
     fn bitCountTwosComp(self: Int) usize {
         var bits = self.bitCountAbs();
 
         // If the entire value has only one bit set (e.g. 0b100000000) then the negation in twos
         // complement requires one less bit.
-        if (!self.positive) block: {
+        if (!self.isPositive()) block: {
             bits += 1;
 
-            if (@popCount(self.limbs[self.len - 1]) == 1) {
-                for (self.limbs[0 .. self.len - 1]) |limb| {
+            if (@popCount(self.limbs[self.len() - 1]) == 1) {
+                for (self.limbs[0 .. self.len() - 1]) |limb| {
                     if (@popCount(limb) != 0) {
                         break :block;
                     }
@@ -154,31 +240,34 @@ pub const Int = struct {
         return bits;
     }
 
-    pub fn fitsInTwosComp(self: Int, is_signed: bool, bit_count: usize) bool {
+    fn fitsInTwosComp(self: Int, is_signed: bool, bit_count: usize) bool {
         if (self.eqZero()) {
             return true;
         }
-        if (!is_signed and !self.positive) {
+        if (!is_signed and !self.isPositive()) {
             return false;
         }
 
-        const req_bits = self.bitCountTwosComp() + @boolToInt(self.positive and is_signed);
+        const req_bits = self.bitCountTwosComp() + @boolToInt(self.isPositive() and is_signed);
         return bit_count >= req_bits;
     }
 
+    /// Returns whether self can fit into an integer of the requested type.
     pub fn fits(self: Int, comptime T: type) bool {
         return self.fitsInTwosComp(T.is_signed, T.bit_count);
     }
 
-    // Returns the approximate size of the integer in the given base. Negative values accommodate for
-    // the minus sign. This is used for determining the number of characters needed to print the
-    // value. It is inexact and will exceed the given value by 1-2 digits.
+    /// Returns the approximate size of the integer in the given base. Negative values accommodate for
+    /// the minus sign. This is used for determining the number of characters needed to print the
+    /// value. It is inexact and may exceed the given value by ~1-2 bytes.
     pub fn sizeInBase(self: Int, base: usize) usize {
-        const bit_count = usize(@boolToInt(!self.positive)) + self.bitCountAbs();
+        const bit_count = usize(@boolToInt(!self.isPositive())) + self.bitCountAbs();
         return (bit_count / math.log2(base)) + 1;
     }
 
+    /// Sets an Int to value. Value must be an primitive integer type.
     pub fn set(self: *Int, value: var) Allocator.Error!void {
+        self.assertWritable();
         const T = @typeOf(value);
 
         switch (@typeInfo(T)) {
@@ -186,19 +275,19 @@ pub const Int = struct {
                 const UT = if (T.is_signed) @IntType(false, T.bit_count - 1) else T;
 
                 try self.ensureCapacity(@sizeOf(UT) / @sizeOf(Limb));
-                self.positive = value >= 0;
-                self.len = 0;
+                self.metadata = 0;
+                self.setSign(value >= 0);
 
                 var w_value: UT = if (value < 0) @intCast(UT, -value) else @intCast(UT, value);
 
                 if (info.bits <= Limb.bit_count) {
                     self.limbs[0] = Limb(w_value);
-                    self.len = 1;
+                    self.metadata += 1;
                 } else {
                     var i: usize = 0;
                     while (w_value != 0) : (i += 1) {
                         self.limbs[i] = @truncate(Limb, w_value);
-                        self.len += 1;
+                        self.metadata += 1;
 
                         // TODO: shift == 64 at compile-time fails. Fails on u128 limbs.
                         w_value >>= Limb.bit_count / 2;
@@ -212,8 +301,8 @@ pub const Int = struct {
                 const req_limbs = @divFloor(math.log2(w_value), Limb.bit_count) + 1;
                 try self.ensureCapacity(req_limbs);
 
-                self.positive = value >= 0;
-                self.len = req_limbs;
+                self.metadata = req_limbs;
+                self.setSign(value >= 0);
 
                 if (w_value <= maxInt(Limb)) {
                     self.limbs[0] = w_value;
@@ -240,6 +329,9 @@ pub const Int = struct {
         TargetTooSmall,
     };
 
+    /// Convert self to type T.
+    ///
+    /// Returns an error if self cannot be narrowed into the requested type without truncation.
     pub fn to(self: Int, comptime T: type) ConvertError!T {
         switch (@typeId(T)) {
             TypeId.Int => {
@@ -254,17 +346,17 @@ pub const Int = struct {
                 if (@sizeOf(UT) <= @sizeOf(Limb)) {
                     r = @intCast(UT, self.limbs[0]);
                 } else {
-                    for (self.limbs[0..self.len]) |_, ri| {
-                        const limb = self.limbs[self.len - ri - 1];
+                    for (self.limbs[0..self.len()]) |_, ri| {
+                        const limb = self.limbs[self.len() - ri - 1];
                         r <<= Limb.bit_count;
                         r |= limb;
                     }
                 }
 
                 if (!T.is_signed) {
-                    return if (self.positive) @intCast(T, r) else error.NegativeIntoUnsigned;
+                    return if (self.isPositive()) @intCast(T, r) else error.NegativeIntoUnsigned;
                 } else {
-                    if (self.positive) {
+                    if (self.isPositive()) {
                         return @intCast(T, r);
                     } else {
                         if (math.cast(T, r)) |ok| {
@@ -303,7 +395,15 @@ pub const Int = struct {
         };
     }
 
+    /// Set self from the string representation `value`.
+    ///
+    /// value must contain only digits <= `base`. Base prefixes are not allowed (e.g. 0x43 should
+    /// simply be 43).
+    ///
+    /// Returns an error if memory could not be allocated or `value` has invalid digits for the
+    /// requested base.
     pub fn setString(self: *Int, base: u8, value: []const u8) !void {
+        self.assertWritable();
         if (base < 2 or base > 16) {
             return error.InvalidBase;
         }
@@ -315,27 +415,22 @@ pub const Int = struct {
             i += 1;
         }
 
-        // TODO values less than limb size should guarantee non allocating
-        var base_buffer: [512]u8 = undefined;
-        const base_al = &std.heap.FixedBufferAllocator.init(base_buffer[0..]).allocator;
-        const base_ap = try Int.initSet(base_al, base);
-
-        var d_buffer: [512]u8 = undefined;
-        var d_fba = std.heap.FixedBufferAllocator.init(d_buffer[0..]);
-        const d_al = &d_fba.allocator;
-
+        const ap_base = Int.initFixed(([]Limb{base})[0..]);
         try self.set(0);
+
         for (value[i..]) |ch| {
             const d = try charToDigit(ch, base);
-            d_fba.end_index = 0;
-            const d_ap = try Int.initSet(d_al, d);
 
-            try self.mul(self.*, base_ap);
-            try self.add(self.*, d_ap);
+            const ap_d = Int.initFixed(([]Limb{d})[0..]);
+
+            try self.mul(self.*, ap_base);
+            try self.add(self.*, ap_d);
         }
-        self.positive = positive;
+        self.setSign(positive);
     }
 
+    /// Converts self to a string in the requested base. Memory is allocated from the provided
+    /// allocator and not the one present in self.
     /// TODO make this call format instead of the other way around
     pub fn toString(self: Int, allocator: *Allocator, base: u8) ![]const u8 {
         if (base < 2 or base > 16) {
@@ -355,7 +450,7 @@ pub const Int = struct {
         if (base & (base - 1) == 0) {
             const base_shift = math.log2_int(Limb, base);
 
-            for (self.limbs[0..self.len]) |limb| {
+            for (self.limbs[0..self.len()]) |limb| {
                 var shift: usize = 0;
                 while (shift < Limb.bit_count) : (shift += base_shift) {
                     const r = @intCast(u8, (limb >> @intCast(Log2Limb, shift)) & Limb(base - 1));
@@ -382,11 +477,11 @@ pub const Int = struct {
             }
 
             var q = try self.clone();
-            q.positive = true;
+            q.abs();
             var r = try Int.init(allocator);
             var b = try Int.initSet(allocator, limb_base);
 
-            while (q.len >= 2) {
+            while (q.len() >= 2) {
                 try Int.divTrunc(&q, &r, q, b);
 
                 var r_word = r.limbs[0];
@@ -399,7 +494,7 @@ pub const Int = struct {
             }
 
             {
-                debug.assert(q.len == 1);
+                debug.assert(q.len() == 1);
 
                 var r_word = q.limbs[0];
                 while (r_word != 0) {
@@ -410,7 +505,7 @@ pub const Int = struct {
             }
         }
 
-        if (!self.positive) {
+        if (!self.isPositive()) {
             try digits.append('-');
         }
 
@@ -419,7 +514,7 @@ pub const Int = struct {
         return s;
     }
 
-    /// for the std lib format function
+    /// To allow `std.fmt.printf` to work with Int.
     /// TODO make this non-allocating
     pub fn format(
         self: Int,
@@ -428,22 +523,24 @@ pub const Int = struct {
         comptime FmtError: type,
         output: fn (@typeOf(context), []const u8) FmtError!void,
     ) FmtError!void {
+        self.assertWritable();
         // TODO look at fmt and support other bases
-        const str = self.toString(self.allocator, 10) catch @panic("TODO make this non allocating");
-        defer self.allocator.free(str);
+        // TODO support read-only fixed integers
+        const str = self.toString(self.allocator.?, 10) catch @panic("TODO make this non allocating");
+        defer self.allocator.?.free(str);
         return output(context, str);
     }
 
-    // returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.
+    /// Returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.
     pub fn cmpAbs(a: Int, b: Int) i8 {
-        if (a.len < b.len) {
+        if (a.len() < b.len()) {
             return -1;
         }
-        if (a.len > b.len) {
+        if (a.len() > b.len()) {
             return 1;
         }
 
-        var i: usize = a.len - 1;
+        var i: usize = a.len() - 1;
         while (i != 0) : (i -= 1) {
             if (a.limbs[i] != b.limbs[i]) {
                 break;
@@ -459,53 +556,37 @@ pub const Int = struct {
         }
     }
 
-    // returns -1, 0, 1 if a < b, a == b or a > b respectively.
+    /// Returns -1, 0, 1 if a < b, a == b or a > b respectively.
     pub fn cmp(a: Int, b: Int) i8 {
-        if (a.positive != b.positive) {
-            return if (a.positive) i8(1) else -1;
+        if (a.isPositive() != b.isPositive()) {
+            return if (a.isPositive()) i8(1) else -1;
         } else {
             const r = cmpAbs(a, b);
-            return if (a.positive) r else -r;
+            return if (a.isPositive()) r else -r;
         }
     }
 
-    // if a == 0
+    /// Returns true if a == 0.
     pub fn eqZero(a: Int) bool {
-        return a.len == 1 and a.limbs[0] == 0;
+        return a.len() == 1 and a.limbs[0] == 0;
     }
 
-    // if |a| == |b|
+    /// Returns true if |a| == |b|.
     pub fn eqAbs(a: Int, b: Int) bool {
         return cmpAbs(a, b) == 0;
     }
 
-    // if a == b
+    /// Returns true if a == b.
     pub fn eq(a: Int, b: Int) bool {
         return cmp(a, b) == 0;
     }
 
-    // Normalize for a possible single carry digit.
-    //
-    // [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
-    // [1, 2, 3, 4, 5] -> [1, 2, 3, 4, 5]
-    // [0]             -> [0]
-    fn norm1(r: *Int, length: usize) void {
-        debug.assert(length > 0);
-        debug.assert(length <= r.limbs.len);
-
-        if (r.limbs[length - 1] == 0) {
-            r.len = if (length > 1) length - 1 else 1;
-        } else {
-            r.len = length;
-        }
-    }
-
     // Normalize a possible sequence of leading zeros.
     //
     // [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
     // [1, 2, 0, 0, 0] -> [1, 2]
     // [0, 0, 0, 0, 0] -> [0]
-    fn normN(r: *Int, length: usize) void {
+    fn normalize(r: *Int, length: usize) void {
         debug.assert(length > 0);
         debug.assert(length <= r.limbs.len);
 
@@ -517,11 +598,25 @@ pub const Int = struct {
         }
 
         // Handle zero
-        r.len = if (j != 0) j else 1;
+        r.setLen(if (j != 0) j else 1);
     }
 
-    // r = a + b
+    // Cannot be used as a result argument to any function.
+    fn readOnlyPositive(a: Int) Int {
+        return Int{
+            .allocator = null,
+            .metadata = a.len(),
+            .limbs = a.limbs,
+        };
+    }
+
+    /// r = a + b
+    ///
+    /// r, a and b may be aliases.
+    ///
+    /// Returns an error if memory could not be allocated.
     pub fn add(r: *Int, a: Int, b: Int) Allocator.Error!void {
+        r.assertWritable();
         if (a.eqZero()) {
             try r.copy(b);
             return;
@@ -530,38 +625,26 @@ pub const Int = struct {
             return;
         }
 
-        if (a.positive != b.positive) {
-            if (a.positive) {
+        if (a.isPositive() != b.isPositive()) {
+            if (a.isPositive()) {
                 // (a) + (-b) => a - b
-                const bp = Int{
-                    .allocator = undefined,
-                    .positive = true,
-                    .limbs = b.limbs,
-                    .len = b.len,
-                };
-                try r.sub(a, bp);
+                try r.sub(a, readOnlyPositive(b));
             } else {
                 // (-a) + (b) => b - a
-                const ap = Int{
-                    .allocator = undefined,
-                    .positive = true,
-                    .limbs = a.limbs,
-                    .len = a.len,
-                };
-                try r.sub(b, ap);
+                try r.sub(b, readOnlyPositive(a));
             }
         } else {
-            if (a.len >= b.len) {
-                try r.ensureCapacity(a.len + 1);
-                lladd(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
-                r.norm1(a.len + 1);
+            if (a.len() >= b.len()) {
+                try r.ensureCapacity(a.len() + 1);
+                lladd(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
+                r.normalize(a.len() + 1);
             } else {
-                try r.ensureCapacity(b.len + 1);
-                lladd(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
-                r.norm1(b.len + 1);
+                try r.ensureCapacity(b.len() + 1);
+                lladd(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
+                r.normalize(b.len() + 1);
             }
 
-            r.positive = a.positive;
+            r.setSign(a.isPositive());
         }
     }
 
@@ -589,55 +672,48 @@ pub const Int = struct {
         r[i] = carry;
     }
 
-    // r = a - b
+    /// r = a - b
+    ///
+    /// r, a and b may be aliases.
+    ///
+    /// Returns an error if memory could not be allocated.
     pub fn sub(r: *Int, a: Int, b: Int) !void {
-        if (a.positive != b.positive) {
-            if (a.positive) {
+        r.assertWritable();
+        if (a.isPositive() != b.isPositive()) {
+            if (a.isPositive()) {
                 // (a) - (-b) => a + b
-                const bp = Int{
-                    .allocator = undefined,
-                    .positive = true,
-                    .limbs = b.limbs,
-                    .len = b.len,
-                };
-                try r.add(a, bp);
+                try r.add(a, readOnlyPositive(b));
             } else {
                 // (-a) - (b) => -(a + b)
-                const ap = Int{
-                    .allocator = undefined,
-                    .positive = true,
-                    .limbs = a.limbs,
-                    .len = a.len,
-                };
-                try r.add(ap, b);
-                r.positive = false;
+                try r.add(readOnlyPositive(a), b);
+                r.setSign(false);
             }
         } else {
-            if (a.positive) {
+            if (a.isPositive()) {
                 // (a) - (b) => a - b
                 if (a.cmp(b) >= 0) {
-                    try r.ensureCapacity(a.len + 1);
-                    llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
-                    r.normN(a.len);
-                    r.positive = true;
+                    try r.ensureCapacity(a.len() + 1);
+                    llsub(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
+                    r.normalize(a.len());
+                    r.setSign(true);
                 } else {
-                    try r.ensureCapacity(b.len + 1);
-                    llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
-                    r.normN(b.len);
-                    r.positive = false;
+                    try r.ensureCapacity(b.len() + 1);
+                    llsub(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
+                    r.normalize(b.len());
+                    r.setSign(false);
                 }
             } else {
                 // (-a) - (-b) => -(a - b)
                 if (a.cmp(b) < 0) {
-                    try r.ensureCapacity(a.len + 1);
-                    llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
-                    r.normN(a.len);
-                    r.positive = false;
+                    try r.ensureCapacity(a.len() + 1);
+                    llsub(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
+                    r.normalize(a.len());
+                    r.setSign(false);
                 } else {
-                    try r.ensureCapacity(b.len + 1);
-                    llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
-                    r.normN(b.len);
-                    r.positive = true;
+                    try r.ensureCapacity(b.len() + 1);
+                    llsub(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
+                    r.normalize(b.len());
+                    r.setSign(true);
                 }
             }
         }
@@ -667,16 +743,20 @@ pub const Int = struct {
         debug.assert(borrow == 0);
     }
 
-    // rma = a * b
-    //
-    // For greatest efficiency, ensure rma does not alias a or b.
+    /// rma = a * b
+    ///
+    /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
+    ///
+    /// Returns an error if memory could not be allocated.
     pub fn mul(rma: *Int, a: Int, b: Int) !void {
+        rma.assertWritable();
+
         var r = rma;
         var aliased = rma.limbs.ptr == a.limbs.ptr or rma.limbs.ptr == b.limbs.ptr;
 
         var sr: Int = undefined;
         if (aliased) {
-            sr = try Int.initCapacity(rma.allocator, a.len + b.len);
+            sr = try Int.initCapacity(rma.allocator.?, a.len() + b.len());
             r = &sr;
             aliased = true;
         }
@@ -685,16 +765,16 @@ pub const Int = struct {
             r.deinit();
         };
 
-        try r.ensureCapacity(a.len + b.len);
+        try r.ensureCapacity(a.len() + b.len());
 
-        if (a.len >= b.len) {
-            llmul(r.limbs, a.limbs[0..a.len], b.limbs[0..b.len]);
+        if (a.len() >= b.len()) {
+            llmul(r.limbs, a.limbs[0..a.len()], b.limbs[0..b.len()]);
         } else {
-            llmul(r.limbs, b.limbs[0..b.len], a.limbs[0..a.len]);
+            llmul(r.limbs, b.limbs[0..b.len()], a.limbs[0..a.len()]);
         }
 
-        r.positive = a.positive == b.positive;
-        r.normN(a.len + b.len);
+        r.normalize(a.len() + b.len());
+        r.setSign(a.isPositive() == b.isPositive());
     }
 
     // a + b * c + *carry, sets carry to the overflow bits
@@ -740,29 +820,34 @@ pub const Int = struct {
         }
     }
 
+    /// q = a / b (rem r)
+    ///
+    /// a / b are floored (rounded towards 0).
     pub fn divFloor(q: *Int, r: *Int, a: Int, b: Int) !void {
         try div(q, r, a, b);
 
         // Trunc -> Floor.
-        if (!q.positive) {
-            // TODO values less than limb size should guarantee non allocating
-            var one_buffer: [512]u8 = undefined;
-            const one_al = &std.heap.FixedBufferAllocator.init(one_buffer[0..]).allocator;
-            const one_ap = try Int.initSet(one_al, 1);
-
-            try q.sub(q.*, one_ap);
-            try r.add(q.*, one_ap);
+        if (!q.isPositive()) {
+            const one = Int.initFixed(([]Limb{1})[0..]);
+            try q.sub(q.*, one);
+            try r.add(q.*, one);
         }
-        r.positive = b.positive;
+        r.setSign(b.isPositive());
     }
 
+    /// q = a / b (rem r)
+    ///
+    /// a / b are truncated (rounded towards -inf).
     pub fn divTrunc(q: *Int, r: *Int, a: Int, b: Int) !void {
         try div(q, r, a, b);
-        r.positive = a.positive;
+        r.setSign(a.isPositive());
     }
 
     // Truncates by default.
     fn div(quo: *Int, rem: *Int, a: Int, b: Int) !void {
+        quo.assertWritable();
+        rem.assertWritable();
+
         if (b.eqZero()) {
             @panic("division by zero");
         }
@@ -773,36 +858,67 @@ pub const Int = struct {
         if (a.cmpAbs(b) < 0) {
             // quo may alias a so handle rem first
             try rem.copy(a);
-            rem.positive = a.positive == b.positive;
+            rem.setSign(a.isPositive() == b.isPositive());
 
-            quo.positive = true;
-            quo.len = 1;
+            quo.metadata = 1;
             quo.limbs[0] = 0;
             return;
         }
 
-        if (b.len == 1) {
-            try quo.ensureCapacity(a.len);
+        // Handle trailing zero-words of divisor/dividend. These are not handled in the following
+        // algorithms.
+        const a_zero_limb_count = blk: {
+            var i: usize = 0;
+            while (i < a.len()) : (i += 1) {
+                if (a.limbs[i] != 0) break;
+            }
+            break :blk i;
+        };
+        const b_zero_limb_count = blk: {
+            var i: usize = 0;
+            while (i < b.len()) : (i += 1) {
+                if (b.limbs[i] != 0) break;
+            }
+            break :blk i;
+        };
 
-            lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[0..a.len], b.limbs[0]);
-            quo.norm1(a.len);
-            quo.positive = a.positive == b.positive;
+        const ab_zero_limb_count = std.math.min(a_zero_limb_count, b_zero_limb_count);
 
-            rem.len = 1;
-            rem.positive = true;
+        if (b.len() - ab_zero_limb_count == 1) {
+            try quo.ensureCapacity(a.len());
+
+            lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.len()], b.limbs[b.len() - 1]);
+            quo.normalize(a.len() - ab_zero_limb_count);
+            quo.setSign(a.isPositive() == b.isPositive());
+
+            rem.metadata = 1;
         } else {
             // x and y are modified during division
-            var x = try a.clone();
+            var x = try Int.initCapacity(quo.allocator.?, a.len());
             defer x.deinit();
+            try x.copy(a);
 
-            var y = try b.clone();
+            var y = try Int.initCapacity(quo.allocator.?, b.len());
             defer y.deinit();
+            try y.copy(b);
 
             // x may grow one limb during normalization
-            try quo.ensureCapacity(a.len + y.len);
-            try divN(quo.allocator, quo, rem, &x, &y);
+            try quo.ensureCapacity(a.len() + y.len());
+
+            // Shrink x, y such that the trailing zero limbs shared between are removed.
+            if (ab_zero_limb_count != 0) {
+                std.mem.copy(Limb, x.limbs[0..], x.limbs[ab_zero_limb_count..]);
+                std.mem.copy(Limb, y.limbs[0..], y.limbs[ab_zero_limb_count..]);
+                x.metadata -= ab_zero_limb_count;
+                y.metadata -= ab_zero_limb_count;
+            }
 
-            quo.positive = a.positive == b.positive;
+            try divN(quo.allocator.?, quo, rem, &x, &y);
+            quo.setSign(a.isPositive() == b.isPositive());
+        }
+
+        if (ab_zero_limb_count != 0) {
+            try rem.shiftLeft(rem.*, ab_zero_limb_count * Limb.bit_count);
         }
     }
 
@@ -837,25 +953,28 @@ pub const Int = struct {
     //
     // x = qy + r where 0 <= r < y
     fn divN(allocator: *Allocator, q: *Int, r: *Int, x: *Int, y: *Int) !void {
-        debug.assert(y.len >= 2);
-        debug.assert(x.len >= y.len);
-        debug.assert(q.limbs.len >= x.len + y.len - 1);
+        debug.assert(y.len() >= 2);
+        debug.assert(x.len() >= y.len());
+        debug.assert(q.limbs.len >= x.len() + y.len() - 1);
         debug.assert(default_capacity >= 3); // see 3.2
 
         var tmp = try Int.init(allocator);
         defer tmp.deinit();
 
-        // Normalize so y > Limb.bit_count / 2 (i.e. leading bit is set)
-        const norm_shift = @clz(y.limbs[y.len - 1]);
+        // Normalize so y > Limb.bit_count / 2 (i.e. leading bit is set) and even
+        var norm_shift = @clz(y.limbs[y.len() - 1]);
+        if (norm_shift == 0 and y.isOdd()) {
+            norm_shift = Limb.bit_count;
+        }
         try x.shiftLeft(x.*, norm_shift);
         try y.shiftLeft(y.*, norm_shift);
 
-        const n = x.len - 1;
-        const t = y.len - 1;
+        const n = x.len() - 1;
+        const t = y.len() - 1;
 
         // 1.
-        q.len = n - t + 1;
-        mem.set(Limb, q.limbs[0..q.len], 0);
+        q.metadata = n - t + 1;
+        mem.set(Limb, q.limbs[0..q.len()], 0);
 
         // 2.
         try tmp.shiftLeft(y.*, Limb.bit_count * (n - t));
@@ -880,7 +999,7 @@ pub const Int = struct {
             tmp.limbs[0] = if (i >= 2) x.limbs[i - 2] else 0;
             tmp.limbs[1] = if (i >= 1) x.limbs[i - 1] else 0;
             tmp.limbs[2] = x.limbs[i];
-            tmp.normN(3);
+            tmp.normalize(3);
 
             while (true) {
                 // 2x1 limb multiplication unrolled against single-limb q[i-t-1]
@@ -888,7 +1007,7 @@ pub const Int = struct {
                 r.limbs[0] = addMulLimbWithCarry(0, if (t >= 1) y.limbs[t - 1] else 0, q.limbs[i - t - 1], &carry);
                 r.limbs[1] = addMulLimbWithCarry(0, y.limbs[t], q.limbs[i - t - 1], &carry);
                 r.limbs[2] = carry;
-                r.normN(3);
+                r.normalize(3);
 
                 if (r.cmpAbs(tmp) <= 0) {
                     break;
@@ -903,7 +1022,7 @@ pub const Int = struct {
             try tmp.shiftLeft(tmp, Limb.bit_count * (i - t - 1));
             try x.sub(x.*, tmp);
 
-            if (!x.positive) {
+            if (!x.isPositive()) {
                 try tmp.shiftLeft(y.*, Limb.bit_count * (i - t - 1));
                 try x.add(x.*, tmp);
                 q.limbs[i - t - 1] -= 1;
@@ -911,18 +1030,20 @@ pub const Int = struct {
         }
 
         // Denormalize
-        q.normN(q.len);
+        q.normalize(q.len());
 
         try r.shiftRight(x.*, norm_shift);
-        r.normN(r.len);
+        r.normalize(r.len());
     }
 
-    // r = a << shift, in other words, r = a * 2^shift
+    /// r = a << shift, in other words, r = a * 2^shift
     pub fn shiftLeft(r: *Int, a: Int, shift: usize) !void {
-        try r.ensureCapacity(a.len + (shift / Limb.bit_count) + 1);
-        llshl(r.limbs[0..], a.limbs[0..a.len], shift);
-        r.norm1(a.len + (shift / Limb.bit_count) + 1);
-        r.positive = a.positive;
+        r.assertWritable();
+
+        try r.ensureCapacity(a.len() + (shift / Limb.bit_count) + 1);
+        llshl(r.limbs[0..], a.limbs[0..a.len()], shift);
+        r.normalize(a.len() + (shift / Limb.bit_count) + 1);
+        r.setSign(a.isPositive());
     }
 
     fn llshl(r: []Limb, a: []const Limb, shift: usize) void {
@@ -948,19 +1069,20 @@ pub const Int = struct {
         mem.set(Limb, r[0 .. limb_shift - 1], 0);
     }
 
-    // r = a >> shift
+    /// r = a >> shift
     pub fn shiftRight(r: *Int, a: Int, shift: usize) !void {
-        if (a.len <= shift / Limb.bit_count) {
-            r.len = 1;
+        r.assertWritable();
+
+        if (a.len() <= shift / Limb.bit_count) {
+            r.metadata = 1;
             r.limbs[0] = 0;
-            r.positive = true;
             return;
         }
 
-        try r.ensureCapacity(a.len - (shift / Limb.bit_count));
-        const r_len = llshr(r.limbs[0..], a.limbs[0..a.len], shift);
-        r.len = a.len - (shift / Limb.bit_count);
-        r.positive = a.positive;
+        try r.ensureCapacity(a.len() - (shift / Limb.bit_count));
+        const r_len = llshr(r.limbs[0..], a.limbs[0..a.len()], shift);
+        r.metadata = a.len() - (shift / Limb.bit_count);
+        r.setSign(a.isPositive());
     }
 
     fn llshr(r: []Limb, a: []const Limb, shift: usize) void {
@@ -983,16 +1105,20 @@ pub const Int = struct {
         }
     }
 
-    // r = a | b
+    /// r = a | b
+    ///
+    /// a and b are zero-extended to the longer of a or b.
     pub fn bitOr(r: *Int, a: Int, b: Int) !void {
-        if (a.len > b.len) {
-            try r.ensureCapacity(a.len);
-            llor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
-            r.len = a.len;
+        r.assertWritable();
+
+        if (a.len() > b.len()) {
+            try r.ensureCapacity(a.len());
+            llor(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
+            r.setLen(a.len());
         } else {
-            try r.ensureCapacity(b.len);
-            llor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
-            r.len = b.len;
+            try r.ensureCapacity(b.len());
+            llor(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
+            r.setLen(b.len());
         }
     }
 
@@ -1010,16 +1136,18 @@ pub const Int = struct {
         }
     }
 
-    // r = a & b
+    /// r = a & b
     pub fn bitAnd(r: *Int, a: Int, b: Int) !void {
-        if (a.len > b.len) {
-            try r.ensureCapacity(b.len);
-            lland(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
-            r.normN(b.len);
+        r.assertWritable();
+
+        if (a.len() > b.len()) {
+            try r.ensureCapacity(b.len());
+            lland(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
+            r.normalize(b.len());
         } else {
-            try r.ensureCapacity(a.len);
-            lland(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
-            r.normN(a.len);
+            try r.ensureCapacity(a.len());
+            lland(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
+            r.normalize(a.len());
         }
     }
 
@@ -1034,16 +1162,18 @@ pub const Int = struct {
         }
     }
 
-    // r = a ^ b
+    /// r = a ^ b
     pub fn bitXor(r: *Int, a: Int, b: Int) !void {
-        if (a.len > b.len) {
-            try r.ensureCapacity(a.len);
-            llxor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
-            r.normN(a.len);
+        r.assertWritable();
+
+        if (a.len() > b.len()) {
+            try r.ensureCapacity(a.len());
+            llxor(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
+            r.normalize(a.len());
         } else {
-            try r.ensureCapacity(b.len);
-            llxor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
-            r.normN(b.len);
+            try r.ensureCapacity(b.len());
+            llxor(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
+            r.normalize(b.len());
         }
     }
 
@@ -1067,7 +1197,9 @@ pub const Int = struct {
 // They will still run on larger than this and should pass, but the multi-limb code-paths
 // may be untested in some cases.
 
-const al = debug.global_allocator;
+var buffer: [64 * 8192]u8 = undefined;
+var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]);
+const al = &fixed.allocator;
 
 test "big.int comptime_int set" {
     comptime var s = 0xefffffff00000001eeeeeeefaaaaaaab;
@@ -1088,14 +1220,14 @@ test "big.int comptime_int set negative" {
     var a = try Int.initSet(al, -10);
 
     testing.expect(a.limbs[0] == 10);
-    testing.expect(a.positive == false);
+    testing.expect(a.isPositive() == false);
 }
 
 test "big.int int set unaligned small" {
     var a = try Int.initSet(al, u7(45));
 
     testing.expect(a.limbs[0] == 45);
-    testing.expect(a.positive == true);
+    testing.expect(a.isPositive() == true);
 }
 
 test "big.int comptime_int to" {
@@ -1116,7 +1248,7 @@ test "big.int to target too small error" {
     testing.expectError(error.TargetTooSmall, a.to(u8));
 }
 
-test "big.int norm1" {
+test "big.int normalize" {
     var a = try Int.init(al);
     try a.ensureCapacity(8);
 
@@ -1124,26 +1256,26 @@ test "big.int norm1" {
     a.limbs[1] = 2;
     a.limbs[2] = 3;
     a.limbs[3] = 0;
-    a.norm1(4);
-    testing.expect(a.len == 3);
+    a.normalize(4);
+    testing.expect(a.len() == 3);
 
     a.limbs[0] = 1;
     a.limbs[1] = 2;
     a.limbs[2] = 3;
-    a.norm1(3);
-    testing.expect(a.len == 3);
+    a.normalize(3);
+    testing.expect(a.len() == 3);
 
     a.limbs[0] = 0;
     a.limbs[1] = 0;
-    a.norm1(2);
-    testing.expect(a.len == 1);
+    a.normalize(2);
+    testing.expect(a.len() == 1);
 
     a.limbs[0] = 0;
-    a.norm1(1);
-    testing.expect(a.len == 1);
+    a.normalize(1);
+    testing.expect(a.len() == 1);
 }
 
-test "big.int normN" {
+test "big.int normalize multi" {
     var a = try Int.init(al);
     try a.ensureCapacity(8);
 
@@ -1151,25 +1283,25 @@ test "big.int normN" {
     a.limbs[1] = 2;
     a.limbs[2] = 0;
     a.limbs[3] = 0;
-    a.normN(4);
-    testing.expect(a.len == 2);
+    a.normalize(4);
+    testing.expect(a.len() == 2);
 
     a.limbs[0] = 1;
     a.limbs[1] = 2;
     a.limbs[2] = 3;
-    a.normN(3);
-    testing.expect(a.len == 3);
+    a.normalize(3);
+    testing.expect(a.len() == 3);
 
     a.limbs[0] = 0;
     a.limbs[1] = 0;
     a.limbs[2] = 0;
     a.limbs[3] = 0;
-    a.normN(4);
-    testing.expect(a.len == 1);
+    a.normalize(4);
+    testing.expect(a.len() == 1);
 
     a.limbs[0] = 0;
-    a.normN(1);
-    testing.expect(a.len == 1);
+    a.normalize(1);
+    testing.expect(a.len() == 1);
 }
 
 test "big.int parity" {
@@ -1204,7 +1336,7 @@ test "big.int bitcount + sizeInBase" {
     try a.shiftLeft(a, 5000);
     testing.expect(a.bitCountAbs() == 5032);
     testing.expect(a.sizeInBase(2) >= 5032);
-    a.positive = false;
+    a.setSign(false);
 
     testing.expect(a.bitCountAbs() == 5032);
     testing.expect(a.sizeInBase(2) >= 5033);
@@ -1216,10 +1348,8 @@ test "big.int bitcount/to" {
     try a.set(0);
     testing.expect(a.bitCountTwosComp() == 0);
 
-    // TODO: stack smashing
-    // testing.expect((try a.to(u0)) == 0);
-    // TODO: sigsegv
-    // testing.expect((try a.to(i0)) == 0);
+    testing.expect((try a.to(u0)) == 0);
+    testing.expect((try a.to(i0)) == 0);
 
     try a.set(-1);
     testing.expect(a.bitCountTwosComp() == 1);
@@ -1980,6 +2110,98 @@ test "big.int div multi-multi (3.1/3.3 branch)" {
     testing.expect((try r.to(u256)) == 0x1111111111111111111110b12222222222222222282);
 }
 
+test "big.int div multi-single zero-limb trailing" {
+    var a = try Int.initSet(al, 0x60000000000000000000000000000000000000000000000000000000000000000);
+    var b = try Int.initSet(al, 0x10000000000000000);
+
+    var q = try Int.init(al);
+    var r = try Int.init(al);
+    try Int.divTrunc(&q, &r, a, b);
+
+    var expected = try Int.initSet(al, 0x6000000000000000000000000000000000000000000000000);
+    testing.expect(q.eq(expected));
+    testing.expect(r.eqZero());
+}
+
+test "big.int div multi-multi zero-limb trailing (with rem)" {
+    var a = try Int.initSet(al, 0x86666666555555558888888777777776111111111111111100000000000000000000000000000000);
+    var b = try Int.initSet(al, 0x8666666655555555444444443333333300000000000000000000000000000000);
+
+    var q = try Int.init(al);
+    var r = try Int.init(al);
+    try Int.divTrunc(&q, &r, a, b);
+
+    testing.expect((try q.to(u128)) == 0x10000000000000000);
+
+    const rs = try r.toString(al, 16);
+    testing.expect(std.mem.eql(u8, rs, "4444444344444443111111111111111100000000000000000000000000000000"));
+}
+
+test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count > divisor zero-limb count" {
+    var a = try Int.initSet(al, 0x8666666655555555888888877777777611111111111111110000000000000000);
+    var b = try Int.initSet(al, 0x8666666655555555444444443333333300000000000000000000000000000000);
+
+    var q = try Int.init(al);
+    var r = try Int.init(al);
+    try Int.divTrunc(&q, &r, a, b);
+
+    testing.expect((try q.to(u128)) == 0x1);
+
+    const rs = try r.toString(al, 16);
+    testing.expect(std.mem.eql(u8, rs, "444444434444444311111111111111110000000000000000"));
+}
+
+test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count < divisor zero-limb count" {
+    var a = try Int.initSet(al, 0x86666666555555558888888777777776111111111111111100000000000000000000000000000000);
+    var b = try Int.initSet(al, 0x866666665555555544444444333333330000000000000000);
+
+    var q = try Int.init(al);
+    var r = try Int.init(al);
+    try Int.divTrunc(&q, &r, a, b);
+
+    const qs = try q.toString(al, 16);
+    testing.expect(std.mem.eql(u8, qs, "10000000000000000820820803105186f"));
+
+    const rs = try r.toString(al, 16);
+    testing.expect(std.mem.eql(u8, rs, "4e11f2baa5896a321d463b543d0104e30000000000000000"));
+}
+
+test "big.int div multi-multi fuzz case #1" {
+    var a = try Int.init(al);
+    var b = try Int.init(al);
+
+    try a.setString(16, "ffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff80000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
+    try b.setString(16, "3ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe0000000000000000000000000000000000001ffffffffffffffffffffffffffffffffffffffffffffffffffc000000000000000000000000000000007fffffffffff");
+
+    var q = try Int.init(al);
+    var r = try Int.init(al);
+    try Int.divTrunc(&q, &r, a, b);
+
+    const qs = try q.toString(al, 16);
+    testing.expect(std.mem.eql(u8, qs, "3ffffffffffffffffffffffffffff0000000000000000000000000000000000001ffffffffffffffffffffffffffff7fffffffe000000000000000000000000000180000000000000000000003fffffbfffffffdfffffffffffffeffff800000100101000000100000000020003fffffdfbfffffe3ffffffffffffeffff7fffc00800a100000017ffe000002000400007efbfff7fe9f00000037ffff3fff7fffa004006100000009ffe00000190038200bf7d2ff7fefe80400060000f7d7f8fbf9401fe38e0403ffc0bdffffa51102c300d7be5ef9df4e5060007b0127ad3fa69f97d0f820b6605ff617ddf7f32ad7a05c0d03f2e7bc78a6000e087a8bbcdc59e07a5a079128a7861f553ddebed7e8e56701756f9ead39b48cd1b0831889ea6ec1fddf643d0565b075ff07e6caea4e2854ec9227fd635ed60a2f5eef2893052ffd54718fa08604acbf6a15e78a467c4a3c53c0278af06c4416573f925491b195e8fd79302cb1aaf7caf4ecfc9aec1254cc969786363ac729f914c6ddcc26738d6b0facd54eba026580aba2eb6482a088b0d224a8852420b91ec1"));
+
+    const rs = try r.toString(al, 16);
+    testing.expect(std.mem.eql(u8, rs, "310d1d4c414426b4836c2635bad1df3a424e50cbdd167ffccb4dfff57d36b4aae0d6ca0910698220171a0f3373c1060a046c2812f0027e321f72979daa5e7973214170d49e885de0c0ecc167837d44502430674a82522e5df6a0759548052420b91ec1"));
+}
+
+test "big.int div multi-multi fuzz case #2" {
+    var a = try Int.init(al);
+    var b = try Int.init(al);
+
+    try a.setString(16, "3ffffffffe00000000000000000000000000fffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000001fffffffffffffffff800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffc000000000000000000000000000000000000000000000000000000000000000");
+    try b.setString(16, "ffc0000000000000000000000000000000000000000000000000");
+
+    var q = try Int.init(al);
+    var r = try Int.init(al);
+    try Int.divTrunc(&q, &r, a, b);
+
+    const qs = try q.toString(al, 16);
+    testing.expect(std.mem.eql(u8, qs, "40100400fe3f8fe3f8fe3f8fe3f8fe3f8fe4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f91e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4992649926499264991e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4792e4b92e4b92e4b92e4b92a4a92a4a92a4"));
+
+    const rs = try r.toString(al, 16);
+    testing.expect(std.mem.eql(u8, rs, "a900000000000000000000000000000000000000000000000000"));
+}
+
 test "big.int shift-right single" {
     var a = try Int.initSet(al, 0xffff0000);
     try a.shiftRight(a, 16);
diff --git a/std/math/big/rational.zig b/std/math/big/rational.zig
new file mode 100644
index 0000000000..58a5e3ac76
--- /dev/null
+++ b/std/math/big/rational.zig
@@ -0,0 +1,938 @@
+const std = @import("../../std.zig");
+const builtin = @import("builtin");
+const debug = std.debug;
+const math = std.math;
+const mem = std.mem;
+const testing = std.testing;
+const Allocator = mem.Allocator;
+const ArrayList = std.ArrayList;
+
+const TypeId = builtin.TypeId;
+
+const bn = @import("int.zig");
+const Limb = bn.Limb;
+const DoubleLimb = bn.DoubleLimb;
+const Int = bn.Int;
+
+/// An arbitrary-precision rational number.
+///
+/// Memory is allocated as needed for operations to ensure full precision is kept. The precision
+/// of a Rational is only bounded by memory.
+///
+/// Rational's are always normalized. That is, for a Rational r = p/q where p and q are integers,
+/// gcd(p, q) = 1 always.
+pub const Rational = struct {
+    /// Numerator. Determines the sign of the Rational.
+    p: Int,
+
+    /// Denominator. Sign is ignored.
+    q: Int,
+
+    /// Create a new Rational. A small amount of memory will be allocated on initialization.
+    /// This will be 2 * Int.default_capacity.
+    pub fn init(a: *Allocator) !Rational {
+        return Rational{
+            .p = try Int.init(a),
+            .q = try Int.initSet(a, 1),
+        };
+    }
+
+    /// Frees all memory associated with a Rational.
+    pub fn deinit(self: *Rational) void {
+        self.p.deinit();
+        self.q.deinit();
+    }
+
+    /// Set a Rational from a primitive integer type.
+    pub fn setInt(self: *Rational, a: var) !void {
+        try self.p.set(a);
+        try self.q.set(1);
+    }
+
+    /// Set a Rational from a string of the form `A/B` where A and B are base-10 integers.
+    pub fn setFloatString(self: *Rational, str: []const u8) !void {
+        // TODO: Accept a/b fractions and exponent form
+        if (str.len == 0) {
+            return error.InvalidFloatString;
+        }
+
+        const State = enum {
+            Integer,
+            Fractional,
+        };
+
+        var state = State.Integer;
+        var point: ?usize = null;
+
+        var start: usize = 0;
+        if (str[0] == '-') {
+            start += 1;
+        }
+
+        for (str) |c, i| {
+            switch (state) {
+                State.Integer => {
+                    switch (c) {
+                        '.' => {
+                            state = State.Fractional;
+                            point = i;
+                        },
+                        '0'...'9' => {
+                            // okay
+                        },
+                        else => {
+                            return error.InvalidFloatString;
+                        },
+                    }
+                },
+                State.Fractional => {
+                    switch (c) {
+                        '0'...'9' => {
+                            // okay
+                        },
+                        else => {
+                            return error.InvalidFloatString;
+                        },
+                    }
+                },
+            }
+        }
+
+        // TODO: batch the multiplies by 10
+        if (point) |i| {
+            try self.p.setString(10, str[0..i]);
+
+            const base = Int.initFixed(([]Limb{10})[0..]);
+
+            var j: usize = start;
+            while (j < str.len - i - 1) : (j += 1) {
+                try self.p.mul(self.p, base);
+            }
+
+            try self.q.setString(10, str[i + 1 ..]);
+            try self.p.add(self.p, self.q);
+
+            try self.q.set(1);
+            var k: usize = i + 1;
+            while (k < str.len) : (k += 1) {
+                try self.q.mul(self.q, base);
+            }
+
+            try self.reduce();
+        } else {
+            try self.p.setString(10, str[0..]);
+            try self.q.set(1);
+        }
+    }
+
+    /// Set a Rational from a floating-point value. The rational will have enough precision to
+    /// completely represent the provided float.
+    pub fn setFloat(self: *Rational, comptime T: type, f: T) !void {
+        // Translated from golang.go/src/math/big/rat.go.
+        debug.assert(@typeId(T) == builtin.TypeId.Float);
+
+        const UnsignedIntType = @IntType(false, T.bit_count);
+        const f_bits = @bitCast(UnsignedIntType, f);
+
+        const exponent_bits = math.floatExponentBits(T);
+        const exponent_bias = (1 << (exponent_bits - 1)) - 1;
+        const mantissa_bits = math.floatMantissaBits(T);
+
+        const exponent_mask = (1 << exponent_bits) - 1;
+        const mantissa_mask = (1 << mantissa_bits) - 1;
+
+        var exponent = @intCast(i16, (f_bits >> mantissa_bits) & exponent_mask);
+        var mantissa = f_bits & mantissa_mask;
+
+        switch (exponent) {
+            exponent_mask => {
+                return error.NonFiniteFloat;
+            },
+            0 => {
+                // denormal
+                exponent -= exponent_bias - 1;
+            },
+            else => {
+                // normal
+                mantissa |= 1 << mantissa_bits;
+                exponent -= exponent_bias;
+            },
+        }
+
+        var shift: i16 = mantissa_bits - exponent;
+
+        // factor out powers of two early from rational
+        while (mantissa & 1 == 0 and shift > 0) {
+            mantissa >>= 1;
+            shift -= 1;
+        }
+
+        try self.p.set(mantissa);
+        self.p.setSign(f >= 0);
+
+        try self.q.set(1);
+        if (shift >= 0) {
+            try self.q.shiftLeft(self.q, @intCast(usize, shift));
+        } else {
+            try self.p.shiftLeft(self.p, @intCast(usize, -shift));
+        }
+
+        try self.reduce();
+    }
+
+    /// Return a floating-point value that is the closest value to a Rational.
+    ///
+    /// The result may not be exact if the Rational is too precise or too large for the
+    /// target type.
+    pub fn toFloat(self: Rational, comptime T: type) !T {
+        // Translated from golang.go/src/math/big/rat.go.
+        // TODO: Indicate whether the result is not exact.
+        debug.assert(@typeId(T) == builtin.TypeId.Float);
+
+        const fsize = T.bit_count;
+        const BitReprType = @IntType(false, T.bit_count);
+
+        const msize = math.floatMantissaBits(T);
+        const msize1 = msize + 1;
+        const msize2 = msize1 + 1;
+
+        const esize = math.floatExponentBits(T);
+        const ebias = (1 << (esize - 1)) - 1;
+        const emin = 1 - ebias;
+        const emax = ebias;
+
+        if (self.p.eqZero()) {
+            return 0;
+        }
+
+        // 1. left-shift a or sub so that a/b is in [1 << msize1, 1 << (msize2 + 1)]
+        var exp = @intCast(isize, self.p.bitCountTwosComp()) - @intCast(isize, self.q.bitCountTwosComp());
+
+        var a2 = try self.p.clone();
+        defer a2.deinit();
+
+        var b2 = try self.q.clone();
+        defer b2.deinit();
+
+        const shift = msize2 - exp;
+        if (shift >= 0) {
+            try a2.shiftLeft(a2, @intCast(usize, shift));
+        } else {
+            try b2.shiftLeft(b2, @intCast(usize, -shift));
+        }
+
+        // 2. compute quotient and remainder
+        var q = try Int.init(self.p.allocator.?);
+        defer q.deinit();
+
+        // unused
+        var r = try Int.init(self.p.allocator.?);
+        defer r.deinit();
+
+        try Int.divTrunc(&q, &r, a2, b2);
+
+        var mantissa = extractLowBits(q, BitReprType);
+        var have_rem = r.len() > 0;
+
+        // 3. q didn't fit in msize2 bits, redo division b2 << 1
+        if (mantissa >> msize2 == 1) {
+            if (mantissa & 1 == 1) {
+                have_rem = true;
+            }
+            mantissa >>= 1;
+            exp += 1;
+        }
+        if (mantissa >> msize1 != 1) {
+            // NOTE: This can be hit if the limb size is small (u8/16).
+            @panic("unexpected bits in result");
+        }
+
+        // 4. Rounding
+        if (emin - msize <= exp and exp <= emin) {
+            // denormal
+            const shift1 = @intCast(math.Log2Int(BitReprType), emin - (exp - 1));
+            const lost_bits = mantissa & ((@intCast(BitReprType, 1) << shift1) - 1);
+            have_rem = have_rem or lost_bits != 0;
+            mantissa >>= shift1;
+            exp = 2 - ebias;
+        }
+
+        // round q using round-half-to-even
+        var exact = !have_rem;
+        if (mantissa & 1 != 0) {
+            exact = false;
+            if (have_rem or (mantissa & 2 != 0)) {
+                mantissa += 1;
+                if (mantissa >= 1 << msize2) {
+                    // 11...1 => 100...0
+                    mantissa >>= 1;
+                    exp += 1;
+                }
+            }
+        }
+        mantissa >>= 1;
+
+        const f = math.scalbn(@intToFloat(T, mantissa), @intCast(i32, exp - msize1));
+        if (math.isInf(f)) {
+            exact = false;
+        }
+
+        return if (self.p.isPositive()) f else -f;
+    }
+
+    /// Set a rational from an integer ratio.
+    pub fn setRatio(self: *Rational, p: var, q: var) !void {
+        try self.p.set(p);
+        try self.q.set(q);
+
+        self.p.setSign(@boolToInt(self.p.isPositive()) ^ @boolToInt(self.q.isPositive()) == 0);
+        self.q.setSign(true);
+
+        try self.reduce();
+
+        if (self.q.eqZero()) {
+            @panic("cannot set rational with denominator = 0");
+        }
+    }
+
+    /// Set a Rational directly from an Int.
+    pub fn copyInt(self: *Rational, a: Int) !void {
+        try self.p.copy(a);
+        try self.q.set(1);
+    }
+
+    /// Set a Rational directly from a ratio of two Int's.
+    pub fn copyRatio(self: *Rational, a: Int, b: Int) !void {
+        try self.p.copy(a);
+        try self.q.copy(b);
+
+        self.p.setSign(@boolToInt(self.p.isPositive()) ^ @boolToInt(self.q.isPositive()) == 0);
+        self.q.setSign(true);
+
+        try self.reduce();
+    }
+
+    /// Make a Rational positive.
+    pub fn abs(r: *Rational) void {
+        r.p.abs();
+    }
+
+    /// Negate the sign of a Rational.
+    pub fn negate(r: *Rational) void {
+        r.p.negate();
+    }
+
+    /// Efficiently swap a Rational with another. This swaps the limb pointers and a full copy is not
+    /// performed. The address of the limbs field will not be the same after this function.
+    pub fn swap(r: *Rational, other: *Rational) void {
+        r.p.swap(&other.p);
+        r.q.swap(&other.q);
+    }
+
+    /// Returns -1, 0, 1 if a < b, a == b or a > b respectively.
+    pub fn cmp(a: Rational, b: Rational) !i8 {
+        return cmpInternal(a, b, true);
+    }
+
+    /// Returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.
+    pub fn cmpAbs(a: Rational, b: Rational) !i8 {
+        return cmpInternal(a, b, false);
+    }
+
+    // p/q > x/y iff p*y > x*q
+    fn cmpInternal(a: Rational, b: Rational, is_abs: bool) !i8 {
+        // TODO: Would a div compare algorithm of sorts be viable and quicker? Can we avoid
+        // the memory allocations here?
+        var q = try Int.init(a.p.allocator.?);
+        defer q.deinit();
+
+        var p = try Int.init(b.p.allocator.?);
+        defer p.deinit();
+
+        try q.mul(a.p, b.q);
+        try p.mul(b.p, a.q);
+
+        return if (is_abs) q.cmpAbs(p) else q.cmp(p);
+    }
+
+    /// rma = a + b.
+    ///
+    /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
+    ///
+    /// Returns an error if memory could not be allocated.
+    pub fn add(rma: *Rational, a: Rational, b: Rational) !void {
+        var r = rma;
+        var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;
+
+        var sr: Rational = undefined;
+        if (aliased) {
+            sr = try Rational.init(rma.p.allocator.?);
+            r = &sr;
+            aliased = true;
+        }
+        defer if (aliased) {
+            rma.swap(r);
+            r.deinit();
+        };
+
+        try r.p.mul(a.p, b.q);
+        try r.q.mul(b.p, a.q);
+        try r.p.add(r.p, r.q);
+
+        try r.q.mul(a.q, b.q);
+        try r.reduce();
+    }
+
+    /// rma = a - b.
+    ///
+    /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
+    ///
+    /// Returns an error if memory could not be allocated.
+    pub fn sub(rma: *Rational, a: Rational, b: Rational) !void {
+        var r = rma;
+        var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;
+
+        var sr: Rational = undefined;
+        if (aliased) {
+            sr = try Rational.init(rma.p.allocator.?);
+            r = &sr;
+            aliased = true;
+        }
+        defer if (aliased) {
+            rma.swap(r);
+            r.deinit();
+        };
+
+        try r.p.mul(a.p, b.q);
+        try r.q.mul(b.p, a.q);
+        try r.p.sub(r.p, r.q);
+
+        try r.q.mul(a.q, b.q);
+        try r.reduce();
+    }
+
+    /// rma = a * b.
+    ///
+    /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
+    ///
+    /// Returns an error if memory could not be allocated.
+    pub fn mul(r: *Rational, a: Rational, b: Rational) !void {
+        try r.p.mul(a.p, b.p);
+        try r.q.mul(a.q, b.q);
+        try r.reduce();
+    }
+
+    /// rma = a / b.
+    ///
+    /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
+    ///
+    /// Returns an error if memory could not be allocated.
+    pub fn div(r: *Rational, a: Rational, b: Rational) !void {
+        if (b.p.eqZero()) {
+            @panic("division by zero");
+        }
+
+        try r.p.mul(a.p, b.q);
+        try r.q.mul(b.p, a.q);
+        try r.reduce();
+    }
+
+    /// Invert the numerator and denominator fields of a Rational. p/q => q/p.
+    pub fn invert(r: *Rational) void {
+        Int.swap(&r.p, &r.q);
+    }
+
+    // reduce r/q such that gcd(r, q) = 1
+    fn reduce(r: *Rational) !void {
+        var a = try Int.init(r.p.allocator.?);
+        defer a.deinit();
+
+        const sign = r.p.isPositive();
+        r.p.abs();
+        try gcd(&a, r.p, r.q);
+        r.p.setSign(sign);
+
+        const one = Int.initFixed(([]Limb{1})[0..]);
+        if (a.cmp(one) != 0) {
+            var unused = try Int.init(r.p.allocator.?);
+            defer unused.deinit();
+
+            // TODO: divexact would be useful here
+            // TODO: don't copy r.q for div
+            try Int.divTrunc(&r.p, &unused, r.p, a);
+            try Int.divTrunc(&r.q, &unused, r.q, a);
+        }
+    }
+};
+
+const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count);
+
+fn gcd(rma: *Int, x: Int, y: Int) !void {
+    rma.assertWritable();
+    var r = rma;
+    var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr;
+
+    var sr: Int = undefined;
+    if (aliased) {
+        sr = try Int.initCapacity(rma.allocator.?, math.max(x.len(), y.len()));
+        r = &sr;
+        aliased = true;
+    }
+    defer if (aliased) {
+        rma.swap(r);
+        r.deinit();
+    };
+
+    try gcdLehmer(r, x, y);
+}
+
+// Storage must live for the lifetime of the returned value
+fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {
+    std.debug.assert(storage.len >= 2);
+
+    var A_is_positive = A >= 0;
+    const Au = @intCast(DoubleLimb, if (A < 0) -A else A);
+    storage[0] = @truncate(Limb, Au);
+    storage[1] = @truncate(Limb, Au >> Limb.bit_count);
+    var Ap = Int.initFixed(storage[0..2]);
+    Ap.setSign(A_is_positive);
+    return Ap;
+}
+
+fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
+    var x = try xa.clone();
+    x.abs();
+    defer x.deinit();
+
+    var y = try ya.clone();
+    y.abs();
+    defer y.deinit();
+
+    if (x.cmp(y) < 0) {
+        x.swap(&y);
+    }
+
+    var T = try Int.init(r.allocator.?);
+    defer T.deinit();
+
+    while (y.len() > 1) {
+        debug.assert(x.isPositive() and y.isPositive());
+        debug.assert(x.len() >= y.len());
+
+        var xh: SignedDoubleLimb = x.limbs[x.len() - 1];
+        var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1];
+
+        var A: SignedDoubleLimb = 1;
+        var B: SignedDoubleLimb = 0;
+        var C: SignedDoubleLimb = 0;
+        var D: SignedDoubleLimb = 1;
+
+        while (yh + C != 0 and yh + D != 0) {
+            const q = @divFloor(xh + A, yh + C);
+            const qp = @divFloor(xh + B, yh + D);
+            if (q != qp) {
+                break;
+            }
+
+            var t = A - q * C;
+            A = C;
+            C = t;
+            t = B - q * D;
+            B = D;
+            D = t;
+
+            t = xh - q * yh;
+            xh = yh;
+            yh = t;
+        }
+
+        if (B == 0) {
+            // T = x % y, r is unused
+            try Int.divTrunc(r, &T, x, y);
+            debug.assert(T.isPositive());
+
+            x.swap(&y);
+            y.swap(&T);
+        } else {
+            var storage: [8]Limb = undefined;
+            const Ap = FixedIntFromSignedDoubleLimb(A, storage[0..2]);
+            const Bp = FixedIntFromSignedDoubleLimb(B, storage[2..4]);
+            const Cp = FixedIntFromSignedDoubleLimb(C, storage[4..6]);
+            const Dp = FixedIntFromSignedDoubleLimb(D, storage[6..8]);
+
+            // T = Ax + By
+            try r.mul(x, Ap);
+            try T.mul(y, Bp);
+            try T.add(r.*, T);
+
+            // u = Cx + Dy, r as u
+            try x.mul(x, Cp);
+            try r.mul(y, Dp);
+            try r.add(x, r.*);
+
+            x.swap(&T);
+            y.swap(r);
+        }
+    }
+
+    // euclidean algorithm
+    debug.assert(x.cmp(y) >= 0);
+
+    while (!y.eqZero()) {
+        try Int.divTrunc(&T, r, x, y);
+        x.swap(&y);
+        y.swap(r);
+    }
+
+    r.swap(&x);
+}
+
+var buffer: [64 * 8192]u8 = undefined;
+var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]);
+var al = &fixed.allocator;
+
+test "big.rational gcd non-one small" {
+    var a = try Int.initSet(al, 17);
+    var b = try Int.initSet(al, 97);
+    var r = try Int.init(al);
+
+    try gcd(&r, a, b);
+
+    testing.expect((try r.to(u32)) == 1);
+}
+
+test "big.rational gcd non-one small" {
+    var a = try Int.initSet(al, 4864);
+    var b = try Int.initSet(al, 3458);
+    var r = try Int.init(al);
+
+    try gcd(&r, a, b);
+
+    testing.expect((try r.to(u32)) == 38);
+}
+
+test "big.rational gcd non-one large" {
+    var a = try Int.initSet(al, 0xffffffffffffffff);
+    var b = try Int.initSet(al, 0xffffffffffffffff7777);
+    var r = try Int.init(al);
+
+    try gcd(&r, a, b);
+
+    testing.expect((try r.to(u32)) == 4369);
+}
+
+test "big.rational gcd large multi-limb result" {
+    var a = try Int.initSet(al, 0x12345678123456781234567812345678123456781234567812345678);
+    var b = try Int.initSet(al, 0x12345671234567123456712345671234567123456712345671234567);
+    var r = try Int.init(al);
+
+    try gcd(&r, a, b);
+
+    testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);
+}
+
+test "big.rational gcd one large" {
+    var a = try Int.initSet(al, 1897056385327307);
+    var b = try Int.initSet(al, 2251799813685248);
+    var r = try Int.init(al);
+
+    try gcd(&r, a, b);
+
+    testing.expect((try r.to(u64)) == 1);
+}
+
+fn extractLowBits(a: Int, comptime T: type) T {
+    testing.expect(@typeId(T) == builtin.TypeId.Int);
+
+    if (T.bit_count <= Limb.bit_count) {
+        return @truncate(T, a.limbs[0]);
+    } else {
+        var r: T = 0;
+        comptime var i: usize = 0;
+
+        // Remainder is always 0 since if T.bit_count >= Limb.bit_count -> Limb | T and both
+        // are powers of two.
+        inline while (i < T.bit_count / Limb.bit_count) : (i += 1) {
+            r |= math.shl(T, a.limbs[i], i * Limb.bit_count);
+        }
+
+        return r;
+    }
+}
+
+test "big.rational extractLowBits" {
+    var a = try Int.initSet(al, 0x11112222333344441234567887654321);
+
+    const a1 = extractLowBits(a, u8);
+    testing.expect(a1 == 0x21);
+
+    const a2 = extractLowBits(a, u16);
+    testing.expect(a2 == 0x4321);
+
+    const a3 = extractLowBits(a, u32);
+    testing.expect(a3 == 0x87654321);
+
+    const a4 = extractLowBits(a, u64);
+    testing.expect(a4 == 0x1234567887654321);
+
+    const a5 = extractLowBits(a, u128);
+    testing.expect(a5 == 0x11112222333344441234567887654321);
+}
+
+test "big.rational set" {
+    var a = try Rational.init(al);
+
+    try a.setInt(5);
+    testing.expect((try a.p.to(u32)) == 5);
+    testing.expect((try a.q.to(u32)) == 1);
+
+    try a.setRatio(7, 3);
+    testing.expect((try a.p.to(u32)) == 7);
+    testing.expect((try a.q.to(u32)) == 3);
+
+    try a.setRatio(9, 3);
+    testing.expect((try a.p.to(i32)) == 3);
+    testing.expect((try a.q.to(i32)) == 1);
+
+    try a.setRatio(-9, 3);
+    testing.expect((try a.p.to(i32)) == -3);
+    testing.expect((try a.q.to(i32)) == 1);
+
+    try a.setRatio(9, -3);
+    testing.expect((try a.p.to(i32)) == -3);
+    testing.expect((try a.q.to(i32)) == 1);
+
+    try a.setRatio(-9, -3);
+    testing.expect((try a.p.to(i32)) == 3);
+    testing.expect((try a.q.to(i32)) == 1);
+}
+
+test "big.rational setFloat" {
+    var a = try Rational.init(al);
+
+    try a.setFloat(f64, 2.5);
+    testing.expect((try a.p.to(i32)) == 5);
+    testing.expect((try a.q.to(i32)) == 2);
+
+    try a.setFloat(f32, -2.5);
+    testing.expect((try a.p.to(i32)) == -5);
+    testing.expect((try a.q.to(i32)) == 2);
+
+    try a.setFloat(f32, 3.141593);
+
+    //                = 3.14159297943115234375
+    testing.expect((try a.p.to(u32)) == 3294199);
+    testing.expect((try a.q.to(u32)) == 1048576);
+
+    try a.setFloat(f64, 72.141593120712409172417410926841290461290467124);
+
+    //                = 72.1415931207124145885245525278151035308837890625
+    testing.expect((try a.p.to(u128)) == 5076513310880537);
+    testing.expect((try a.q.to(u128)) == 70368744177664);
+}
+
+test "big.rational setFloatString" {
+    var a = try Rational.init(al);
+
+    try a.setFloatString("72.14159312071241458852455252781510353");
+
+    //                  = 72.1415931207124145885245525278151035308837890625
+    testing.expect((try a.p.to(u128)) == 7214159312071241458852455252781510353);
+    testing.expect((try a.q.to(u128)) == 100000000000000000000000000000000000);
+}
+
+test "big.rational toFloat" {
+    var a = try Rational.init(al);
+
+    // = 3.14159297943115234375
+    try a.setRatio(3294199, 1048576);
+    testing.expect((try a.toFloat(f64)) == 3.14159297943115234375);
+
+    // = 72.1415931207124145885245525278151035308837890625
+    try a.setRatio(5076513310880537, 70368744177664);
+    testing.expect((try a.toFloat(f64)) == 72.141593120712409172417410926841290461290467124);
+}
+
+test "big.rational set/to Float round-trip" {
+    var a = try Rational.init(al);
+    var prng = std.rand.DefaultPrng.init(0x5EED);
+    var i: usize = 0;
+    while (i < 512) : (i += 1) {
+        const r = prng.random.float(f64);
+        try a.setFloat(f64, r);
+        testing.expect((try a.toFloat(f64)) == r);
+    }
+}
+
+test "big.rational copy" {
+    var a = try Rational.init(al);
+
+    const b = try Int.initSet(al, 5);
+
+    try a.copyInt(b);
+    testing.expect((try a.p.to(u32)) == 5);
+    testing.expect((try a.q.to(u32)) == 1);
+
+    const c = try Int.initSet(al, 7);
+    const d = try Int.initSet(al, 3);
+
+    try a.copyRatio(c, d);
+    testing.expect((try a.p.to(u32)) == 7);
+    testing.expect((try a.q.to(u32)) == 3);
+
+    const e = try Int.initSet(al, 9);
+    const f = try Int.initSet(al, 3);
+
+    try a.copyRatio(e, f);
+    testing.expect((try a.p.to(u32)) == 3);
+    testing.expect((try a.q.to(u32)) == 1);
+}
+
+test "big.rational negate" {
+    var a = try Rational.init(al);
+
+    try a.setInt(-50);
+    testing.expect((try a.p.to(i32)) == -50);
+    testing.expect((try a.q.to(i32)) == 1);
+
+    a.negate();
+    testing.expect((try a.p.to(i32)) == 50);
+    testing.expect((try a.q.to(i32)) == 1);
+
+    a.negate();
+    testing.expect((try a.p.to(i32)) == -50);
+    testing.expect((try a.q.to(i32)) == 1);
+}
+
+test "big.rational abs" {
+    var a = try Rational.init(al);
+
+    try a.setInt(-50);
+    testing.expect((try a.p.to(i32)) == -50);
+    testing.expect((try a.q.to(i32)) == 1);
+
+    a.abs();
+    testing.expect((try a.p.to(i32)) == 50);
+    testing.expect((try a.q.to(i32)) == 1);
+
+    a.abs();
+    testing.expect((try a.p.to(i32)) == 50);
+    testing.expect((try a.q.to(i32)) == 1);
+}
+
+test "big.rational swap" {
+    var a = try Rational.init(al);
+    var b = try Rational.init(al);
+
+    try a.setRatio(50, 23);
+    try b.setRatio(17, 3);
+
+    testing.expect((try a.p.to(u32)) == 50);
+    testing.expect((try a.q.to(u32)) == 23);
+
+    testing.expect((try b.p.to(u32)) == 17);
+    testing.expect((try b.q.to(u32)) == 3);
+
+    a.swap(&b);
+
+    testing.expect((try a.p.to(u32)) == 17);
+    testing.expect((try a.q.to(u32)) == 3);
+
+    testing.expect((try b.p.to(u32)) == 50);
+    testing.expect((try b.q.to(u32)) == 23);
+}
+
+test "big.rational cmp" {
+    var a = try Rational.init(al);
+    var b = try Rational.init(al);
+
+    try a.setRatio(500, 231);
+    try b.setRatio(18903, 8584);
+    testing.expect((try a.cmp(b)) < 0);
+
+    try a.setRatio(890, 10);
+    try b.setRatio(89, 1);
+    testing.expect((try a.cmp(b)) == 0);
+}
+
+test "big.rational add single-limb" {
+    var a = try Rational.init(al);
+    var b = try Rational.init(al);
+
+    try a.setRatio(500, 231);
+    try b.setRatio(18903, 8584);
+    testing.expect((try a.cmp(b)) < 0);
+
+    try a.setRatio(890, 10);
+    try b.setRatio(89, 1);
+    testing.expect((try a.cmp(b)) == 0);
+}
+
+test "big.rational add" {
+    var a = try Rational.init(al);
+    var b = try Rational.init(al);
+    var r = try Rational.init(al);
+
+    try a.setRatio(78923, 23341);
+    try b.setRatio(123097, 12441414);
+    try a.add(a, b);
+
+    try r.setRatio(984786924199, 290395044174);
+    testing.expect((try a.cmp(r)) == 0);
+}
+
+test "big.rational sub" {
+    var a = try Rational.init(al);
+    var b = try Rational.init(al);
+    var r = try Rational.init(al);
+
+    try a.setRatio(78923, 23341);
+    try b.setRatio(123097, 12441414);
+    try a.sub(a, b);
+
+    try r.setRatio(979040510045, 290395044174);
+    testing.expect((try a.cmp(r)) == 0);
+}
+
+test "big.rational mul" {
+    var a = try Rational.init(al);
+    var b = try Rational.init(al);
+    var r = try Rational.init(al);
+
+    try a.setRatio(78923, 23341);
+    try b.setRatio(123097, 12441414);
+    try a.mul(a, b);
+
+    try r.setRatio(571481443, 17082061422);
+    testing.expect((try a.cmp(r)) == 0);
+}
+
+test "big.rational div" {
+    var a = try Rational.init(al);
+    var b = try Rational.init(al);
+    var r = try Rational.init(al);
+
+    try a.setRatio(78923, 23341);
+    try b.setRatio(123097, 12441414);
+    try a.div(a, b);
+
+    try r.setRatio(75531824394, 221015929);
+    testing.expect((try a.cmp(r)) == 0);
+}
+
+test "big.rational div" {
+    var a = try Rational.init(al);
+    var r = try Rational.init(al);
+
+    try a.setRatio(78923, 23341);
+    a.invert();
+
+    try r.setRatio(23341, 78923);
+    testing.expect((try a.cmp(r)) == 0);
+
+    try a.setRatio(-78923, 23341);
+    a.invert();
+
+    try r.setRatio(-23341, 78923);
+    testing.expect((try a.cmp(r)) == 0);
+}
diff --git a/std/math/cbrt.zig b/std/math/cbrt.zig
index 9497255269..5241e31323 100644
--- a/std/math/cbrt.zig
+++ b/std/math/cbrt.zig
@@ -1,13 +1,19 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - cbrt(+-0)   = +-0
-// - cbrt(+-inf) = +-inf
-// - cbrt(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/cbrtf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/cbrt.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the cube root of x.
+///
+/// Special Cases:
+///  - cbrt(+-0)   = +-0
+///  - cbrt(+-inf) = +-inf
+///  - cbrt(nan)   = nan
 pub fn cbrt(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/ceil.zig b/std/math/ceil.zig
index da83e0ec59..5f86093a6d 100644
--- a/std/math/ceil.zig
+++ b/std/math/ceil.zig
@@ -1,14 +1,20 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - ceil(+-0)   = +-0
-// - ceil(+-inf) = +-inf
-// - ceil(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/ceilf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/ceil.c
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the least integer value greater than of equal to x.
+///
+/// Special Cases:
+///  - ceil(+-0)   = +-0
+///  - ceil(+-inf) = +-inf
+///  - ceil(nan)   = nan
 pub fn ceil(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/complex.zig b/std/math/complex.zig
index cc0573b227..e5574f9cee 100644
--- a/std/math/complex.zig
+++ b/std/math/complex.zig
@@ -23,13 +23,18 @@ pub const sqrt = @import("complex/sqrt.zig").sqrt;
 pub const tanh = @import("complex/tanh.zig").tanh;
 pub const tan = @import("complex/tan.zig").tan;
 
+/// A complex number consisting of a real an imaginary part. T must be a floating-point value.
 pub fn Complex(comptime T: type) type {
     return struct {
         const Self = @This();
 
+        /// Real part.
         re: T,
+
+        /// Imaginary part.
         im: T,
 
+        /// Create a new Complex number from the given real and imaginary parts.
         pub fn new(re: T, im: T) Self {
             return Self{
                 .re = re,
@@ -37,6 +42,7 @@ pub fn Complex(comptime T: type) type {
             };
         }
 
+        /// Returns the sum of two complex numbers.
         pub fn add(self: Self, other: Self) Self {
             return Self{
                 .re = self.re + other.re,
@@ -44,6 +50,7 @@ pub fn Complex(comptime T: type) type {
             };
         }
 
+        /// Returns the subtraction of two complex numbers.
         pub fn sub(self: Self, other: Self) Self {
             return Self{
                 .re = self.re - other.re,
@@ -51,6 +58,7 @@ pub fn Complex(comptime T: type) type {
             };
         }
 
+        /// Returns the product of two complex numbers.
         pub fn mul(self: Self, other: Self) Self {
             return Self{
                 .re = self.re * other.re - self.im * other.im,
@@ -58,6 +66,7 @@ pub fn Complex(comptime T: type) type {
             };
         }
 
+        /// Returns the quotient of two complex numbers.
         pub fn div(self: Self, other: Self) Self {
             const re_num = self.re * other.re + self.im * other.im;
             const im_num = self.im * other.re - self.re * other.im;
@@ -69,6 +78,7 @@ pub fn Complex(comptime T: type) type {
             };
         }
 
+        /// Returns the complex conjugate of a number.
         pub fn conjugate(self: Self) Self {
             return Self{
                 .re = self.re,
@@ -76,6 +86,7 @@ pub fn Complex(comptime T: type) type {
             };
         }
 
+        /// Returns the reciprocal of a complex number.
         pub fn reciprocal(self: Self) Self {
             const m = self.re * self.re + self.im * self.im;
             return Self{
@@ -84,6 +95,7 @@ pub fn Complex(comptime T: type) type {
             };
         }
 
+        /// Returns the magnitude of a complex number.
         pub fn magnitude(self: Self) T {
             return math.sqrt(self.re * self.re + self.im * self.im);
         }
diff --git a/std/math/complex/abs.zig b/std/math/complex/abs.zig
index e1368d6ef6..8105f57218 100644
--- a/std/math/complex/abs.zig
+++ b/std/math/complex/abs.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the absolute value (modulus) of z.
 pub fn abs(z: var) @typeOf(z.re) {
     const T = @typeOf(z.re);
     return math.hypot(T, z.re, z.im);
diff --git a/std/math/complex/acos.zig b/std/math/complex/acos.zig
index 8aed26a71b..f3526cc9ff 100644
--- a/std/math/complex/acos.zig
+++ b/std/math/complex/acos.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the arc-cosine of z.
 pub fn acos(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const q = cmath.asin(z);
diff --git a/std/math/complex/acosh.zig b/std/math/complex/acosh.zig
index e72bf431fe..6f0fd2e36c 100644
--- a/std/math/complex/acosh.zig
+++ b/std/math/complex/acosh.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the hyperbolic arc-cosine of z.
 pub fn acosh(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const q = cmath.acos(z);
diff --git a/std/math/complex/arg.zig b/std/math/complex/arg.zig
index 0a2441d1fd..d0c9588b8d 100644
--- a/std/math/complex/arg.zig
+++ b/std/math/complex/arg.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the angular component (in radians) of z.
 pub fn arg(z: var) @typeOf(z.re) {
     const T = @typeOf(z.re);
     return math.atan2(T, z.im, z.re);
diff --git a/std/math/complex/asin.zig b/std/math/complex/asin.zig
index 6be775d748..76f94a286c 100644
--- a/std/math/complex/asin.zig
+++ b/std/math/complex/asin.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+// Returns the arc-sine of z.
 pub fn asin(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const x = z.re;
diff --git a/std/math/complex/asinh.zig b/std/math/complex/asinh.zig
index 8e09036750..da065aad01 100644
--- a/std/math/complex/asinh.zig
+++ b/std/math/complex/asinh.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the hyperbolic arc-sine of z.
 pub fn asinh(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const q = Complex(T).new(-z.im, z.re);
diff --git a/std/math/complex/atan.zig b/std/math/complex/atan.zig
index 6b83adbd97..89bc8dfaf0 100644
--- a/std/math/complex/atan.zig
+++ b/std/math/complex/atan.zig
@@ -1,9 +1,16 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/catanf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/catan.c
+
 const std = @import("../../std.zig");
 const testing = std.testing;
 const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the arc-tangent of z.
 pub fn atan(z: var) @typeOf(z) {
     const T = @typeOf(z.re);
     return switch (T) {
diff --git a/std/math/complex/atanh.zig b/std/math/complex/atanh.zig
index 8edfb6e78e..225e7c61de 100644
--- a/std/math/complex/atanh.zig
+++ b/std/math/complex/atanh.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the hyperbolic arc-tangent of z.
 pub fn atanh(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const q = Complex(T).new(-z.im, z.re);
diff --git a/std/math/complex/conj.zig b/std/math/complex/conj.zig
index 7a42d365ea..bd71ca3c06 100644
--- a/std/math/complex/conj.zig
+++ b/std/math/complex/conj.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the complex conjugate of z.
 pub fn conj(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     return Complex(T).new(z.re, -z.im);
diff --git a/std/math/complex/cos.zig b/std/math/complex/cos.zig
index 71f9603c75..332009ffe5 100644
--- a/std/math/complex/cos.zig
+++ b/std/math/complex/cos.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the cosine of z.
 pub fn cos(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const p = Complex(T).new(-z.im, z.re);
diff --git a/std/math/complex/cosh.zig b/std/math/complex/cosh.zig
index 9806e41170..be7bfde963 100644
--- a/std/math/complex/cosh.zig
+++ b/std/math/complex/cosh.zig
@@ -1,3 +1,9 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/ccoshf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/ccosh.c
+
 const std = @import("../../std.zig");
 const testing = std.testing;
 const math = std.math;
@@ -6,6 +12,7 @@ const Complex = cmath.Complex;
 
 const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
 
+/// Returns the hyperbolic arc-cosine of z.
 pub fn cosh(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     return switch (T) {
diff --git a/std/math/complex/exp.zig b/std/math/complex/exp.zig
index c74ac2fc08..9b686bebc3 100644
--- a/std/math/complex/exp.zig
+++ b/std/math/complex/exp.zig
@@ -1,3 +1,9 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/cexpf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/cexp.c
+
 const std = @import("../../std.zig");
 const testing = std.testing;
 const math = std.math;
@@ -6,6 +12,7 @@ const Complex = cmath.Complex;
 
 const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
 
+/// Returns e raised to the power of z (e^z).
 pub fn exp(z: var) @typeOf(z) {
     const T = @typeOf(z.re);
 
diff --git a/std/math/complex/ldexp.zig b/std/math/complex/ldexp.zig
index 6b4306bf77..d6f810793f 100644
--- a/std/math/complex/ldexp.zig
+++ b/std/math/complex/ldexp.zig
@@ -1,9 +1,16 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/__cexpf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/__cexp.c
+
 const std = @import("../../std.zig");
 const debug = std.debug;
 const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns exp(z) scaled to avoid overflow.
 pub fn ldexp_cexp(z: var, expt: i32) @typeOf(z) {
     const T = @typeOf(z.re);
 
diff --git a/std/math/complex/log.zig b/std/math/complex/log.zig
index 2b43a6970f..762b4fde9a 100644
--- a/std/math/complex/log.zig
+++ b/std/math/complex/log.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the natural logarithm of z.
 pub fn log(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const r = cmath.abs(z);
diff --git a/std/math/complex/pow.zig b/std/math/complex/pow.zig
index 9174bb3626..a2480453fc 100644
--- a/std/math/complex/pow.zig
+++ b/std/math/complex/pow.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns z raised to the complex power of c.
 pub fn pow(comptime T: type, z: T, c: T) T {
     const p = cmath.log(z);
     const q = c.mul(p);
diff --git a/std/math/complex/proj.zig b/std/math/complex/proj.zig
index aadcff6ff6..c8f2d9fc6d 100644
--- a/std/math/complex/proj.zig
+++ b/std/math/complex/proj.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the projection of z onto the riemann sphere.
 pub fn proj(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
 
diff --git a/std/math/complex/sin.zig b/std/math/complex/sin.zig
index 049631f31e..9ddc3a7a80 100644
--- a/std/math/complex/sin.zig
+++ b/std/math/complex/sin.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the sine of z.
 pub fn sin(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const p = Complex(T).new(-z.im, z.re);
diff --git a/std/math/complex/sinh.zig b/std/math/complex/sinh.zig
index 0b656e5354..6286d8447f 100644
--- a/std/math/complex/sinh.zig
+++ b/std/math/complex/sinh.zig
@@ -1,3 +1,9 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/csinhf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/csinh.c
+
 const std = @import("../../std.zig");
 const testing = std.testing;
 const math = std.math;
@@ -6,6 +12,7 @@ const Complex = cmath.Complex;
 
 const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
 
+/// Returns the hyperbolic sine of z.
 pub fn sinh(z: var) @typeOf(z) {
     const T = @typeOf(z.re);
     return switch (T) {
diff --git a/std/math/complex/sqrt.zig b/std/math/complex/sqrt.zig
index e935d0b238..36f4c28e29 100644
--- a/std/math/complex/sqrt.zig
+++ b/std/math/complex/sqrt.zig
@@ -1,9 +1,17 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/csqrtf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/csqrt.c
+
 const std = @import("../../std.zig");
 const testing = std.testing;
 const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the square root of z. The real and imaginary parts of the result have the same sign
+/// as the imaginary part of z.
 pub fn sqrt(z: var) @typeOf(z) {
     const T = @typeOf(z.re);
 
diff --git a/std/math/complex/tan.zig b/std/math/complex/tan.zig
index 45e2873eb6..398b8295ca 100644
--- a/std/math/complex/tan.zig
+++ b/std/math/complex/tan.zig
@@ -4,6 +4,7 @@ const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the tanget of z.
 pub fn tan(z: var) Complex(@typeOf(z.re)) {
     const T = @typeOf(z.re);
     const q = Complex(T).new(-z.im, z.re);
diff --git a/std/math/complex/tanh.zig b/std/math/complex/tanh.zig
index de905ee3f6..5c14ec66f2 100644
--- a/std/math/complex/tanh.zig
+++ b/std/math/complex/tanh.zig
@@ -1,9 +1,16 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/ctanhf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/complex/ctanh.c
+
 const std = @import("../../std.zig");
 const testing = std.testing;
 const math = std.math;
 const cmath = math.complex;
 const Complex = cmath.Complex;
 
+/// Returns the hyperbolic tangent of z.
 pub fn tanh(z: var) @typeOf(z) {
     const T = @typeOf(z.re);
     return switch (T) {
diff --git a/std/math/copysign.zig b/std/math/copysign.zig
index ff8140b3ab..e4d90c395e 100644
--- a/std/math/copysign.zig
+++ b/std/math/copysign.zig
@@ -1,8 +1,15 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/math/copysignf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/copysign.c
+
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns a value with the magnitude of x and the sign of y.
 pub fn copysign(comptime T: type, x: T, y: T) T {
     return switch (T) {
         f16 => copysign16(x, y),
diff --git a/std/math/cos.zig b/std/math/cos.zig
index 9479482894..dc4aff59d6 100644
--- a/std/math/cos.zig
+++ b/std/math/cos.zig
@@ -1,18 +1,23 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
 //
-// - cos(+-inf) = nan
-// - cos(nan)   = nan
+// https://golang.org/src/math/sin.go
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the cosine of the radian value x.
+///
+/// Special Cases:
+///  - cos(+-inf) = nan
+///  - cos(nan)   = nan
 pub fn cos(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
-        f32 => cos32(x),
-        f64 => cos64(x),
+        f32 => cos_(f32, x),
+        f64 => cos_(f64, x),
         else => @compileError("cos not implemented for " ++ @typeName(T)),
     };
 }
@@ -33,78 +38,24 @@ const C3 = 2.48015872888517045348E-5;
 const C4 = -1.38888888888730564116E-3;
 const C5 = 4.16666666666665929218E-2;
 
-// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.
-//
-// This may have slight differences on some edge cases and may need to replaced if so.
-fn cos32(x_: f32) f32 {
-    const pi4a = 7.85398125648498535156e-1;
-    const pi4b = 3.77489470793079817668E-8;
-    const pi4c = 2.69515142907905952645E-15;
-    const m4pi = 1.273239544735162542821171882678754627704620361328125;
-
-    var x = x_;
-    if (math.isNan(x) or math.isInf(x)) {
-        return math.nan(f32);
-    }
-
-    var sign = false;
-    if (x < 0) {
-        x = -x;
-    }
-
-    var y = math.floor(x * m4pi);
-    var j = @floatToInt(i64, y);
+const pi4a = 7.85398125648498535156e-1;
+const pi4b = 3.77489470793079817668E-8;
+const pi4c = 2.69515142907905952645E-15;
+const m4pi = 1.273239544735162542821171882678754627704620361328125;
 
-    if (j & 1 == 1) {
-        j += 1;
-        y += 1;
-    }
-
-    j &= 7;
-    if (j > 3) {
-        j -= 4;
-        sign = !sign;
-    }
-    if (j > 1) {
-        sign = !sign;
-    }
-
-    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
-    const w = z * z;
-
-    const r = r: {
-        if (j == 1 or j == 2) {
-            break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
-        } else {
-            break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
-        }
-    };
-
-    if (sign) {
-        return -r;
-    } else {
-        return r;
-    }
-}
-
-fn cos64(x_: f64) f64 {
-    const pi4a = 7.85398125648498535156e-1;
-    const pi4b = 3.77489470793079817668E-8;
-    const pi4c = 2.69515142907905952645E-15;
-    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+fn cos_(comptime T: type, x_: T) T {
+    const I = @IntType(true, T.bit_count);
 
     var x = x_;
     if (math.isNan(x) or math.isInf(x)) {
-        return math.nan(f64);
+        return math.nan(T);
     }
 
     var sign = false;
-    if (x < 0) {
-        x = -x;
-    }
+    x = math.fabs(x);
 
     var y = math.floor(x * m4pi);
-    var j = @floatToInt(i64, y);
+    var j = @floatToInt(I, y);
 
     if (j & 1 == 1) {
         j += 1;
@@ -123,56 +74,51 @@ fn cos64(x_: f64) f64 {
     const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
     const w = z * z;
 
-    const r = r: {
-        if (j == 1 or j == 2) {
-            break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
-        } else {
-            break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
-        }
-    };
+    const r = if (j == 1 or j == 2)
+        z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))))
+    else
+        1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
 
-    if (sign) {
-        return -r;
-    } else {
-        return r;
-    }
+    return if (sign) -r else r;
 }
 
 test "math.cos" {
-    expect(cos(f32(0.0)) == cos32(0.0));
-    expect(cos(f64(0.0)) == cos64(0.0));
+    expect(cos(f32(0.0)) == cos_(f32, 0.0));
+    expect(cos(f64(0.0)) == cos_(f64, 0.0));
 }
 
 test "math.cos32" {
     const epsilon = 0.000001;
 
-    expect(math.approxEq(f32, cos32(0.0), 1.0, epsilon));
-    expect(math.approxEq(f32, cos32(0.2), 0.980067, epsilon));
-    expect(math.approxEq(f32, cos32(0.8923), 0.627623, epsilon));
-    expect(math.approxEq(f32, cos32(1.5), 0.070737, epsilon));
-    expect(math.approxEq(f32, cos32(37.45), 0.969132, epsilon));
-    expect(math.approxEq(f32, cos32(89.123), 0.400798, epsilon));
+    expect(math.approxEq(f32, cos_(f32, 0.0), 1.0, epsilon));
+    expect(math.approxEq(f32, cos_(f32, 0.2), 0.980067, epsilon));
+    expect(math.approxEq(f32, cos_(f32, 0.8923), 0.627623, epsilon));
+    expect(math.approxEq(f32, cos_(f32, 1.5), 0.070737, epsilon));
+    expect(math.approxEq(f32, cos_(f32, -1.5), 0.070737, epsilon));
+    expect(math.approxEq(f32, cos_(f32, 37.45), 0.969132, epsilon));
+    expect(math.approxEq(f32, cos_(f32, 89.123), 0.400798, epsilon));
 }
 
 test "math.cos64" {
     const epsilon = 0.000001;
 
-    expect(math.approxEq(f64, cos64(0.0), 1.0, epsilon));
-    expect(math.approxEq(f64, cos64(0.2), 0.980067, epsilon));
-    expect(math.approxEq(f64, cos64(0.8923), 0.627623, epsilon));
-    expect(math.approxEq(f64, cos64(1.5), 0.070737, epsilon));
-    expect(math.approxEq(f64, cos64(37.45), 0.969132, epsilon));
-    expect(math.approxEq(f64, cos64(89.123), 0.40080, epsilon));
+    expect(math.approxEq(f64, cos_(f64, 0.0), 1.0, epsilon));
+    expect(math.approxEq(f64, cos_(f64, 0.2), 0.980067, epsilon));
+    expect(math.approxEq(f64, cos_(f64, 0.8923), 0.627623, epsilon));
+    expect(math.approxEq(f64, cos_(f64, 1.5), 0.070737, epsilon));
+    expect(math.approxEq(f64, cos_(f64, -1.5), 0.070737, epsilon));
+    expect(math.approxEq(f64, cos_(f64, 37.45), 0.969132, epsilon));
+    expect(math.approxEq(f64, cos_(f64, 89.123), 0.40080, epsilon));
 }
 
 test "math.cos32.special" {
-    expect(math.isNan(cos32(math.inf(f32))));
-    expect(math.isNan(cos32(-math.inf(f32))));
-    expect(math.isNan(cos32(math.nan(f32))));
+    expect(math.isNan(cos_(f32, math.inf(f32))));
+    expect(math.isNan(cos_(f32, -math.inf(f32))));
+    expect(math.isNan(cos_(f32, math.nan(f32))));
 }
 
 test "math.cos64.special" {
-    expect(math.isNan(cos64(math.inf(f64))));
-    expect(math.isNan(cos64(-math.inf(f64))));
-    expect(math.isNan(cos64(math.nan(f64))));
+    expect(math.isNan(cos_(f64, math.inf(f64))));
+    expect(math.isNan(cos_(f64, -math.inf(f64))));
+    expect(math.isNan(cos_(f64, math.nan(f64))));
 }
diff --git a/std/math/cosh.zig b/std/math/cosh.zig
index eb3082e5fa..75c5c15ec1 100644
--- a/std/math/cosh.zig
+++ b/std/math/cosh.zig
@@ -1,8 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - cosh(+-0)   = 1
-// - cosh(+-inf) = +inf
-// - cosh(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/coshf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/cosh.c
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
@@ -11,6 +11,12 @@ const expo2 = @import("expo2.zig").expo2;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns the hyperbolic cosine of x.
+///
+/// Special Cases:
+///  - cosh(+-0)   = 1
+///  - cosh(+-inf) = +inf
+///  - cosh(nan)   = nan
 pub fn cosh(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/exp.zig b/std/math/exp.zig
index aabccffd0b..ad058646a4 100644
--- a/std/math/exp.zig
+++ b/std/math/exp.zig
@@ -1,13 +1,19 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - exp(+inf) = +inf
-// - exp(nan)  = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/expf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/exp.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const assert = std.debug.assert;
 const builtin = @import("builtin");
 
+/// Returns e raised to the power of x (e^x).
+///
+/// Special Cases:
+///  - exp(+inf) = +inf
+///  - exp(nan)  = nan
 pub fn exp(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/exp2.zig b/std/math/exp2.zig
index 392d45bf68..07a39576b1 100644
--- a/std/math/exp2.zig
+++ b/std/math/exp2.zig
@@ -1,12 +1,18 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - exp2(+inf) = +inf
-// - exp2(nan)  = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/exp2f.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/exp2.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns 2 raised to the power of x (2^x).
+///
+/// Special Cases:
+///  - exp2(+inf) = +inf
+///  - exp2(nan)  = nan
 pub fn exp2(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/expm1.zig b/std/math/expm1.zig
index b83cc3e82e..5e347f86f6 100644
--- a/std/math/expm1.zig
+++ b/std/math/expm1.zig
@@ -1,14 +1,23 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - expm1(+inf) = +inf
-// - expm1(-inf) = -1
-// - expm1(nan)  = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/expmf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/expm.c
+
+// TODO: Updated recently.
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns e raised to the power of x, minus 1 (e^x - 1). This is more accurate than exp(e, x) - 1
+/// when x is near 0.
+///
+/// Special Cases:
+///  - expm1(+inf) = +inf
+///  - expm1(-inf) = -1
+///  - expm1(nan)  = nan
 pub fn expm1(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/expo2.zig b/std/math/expo2.zig
index 57e2bf34c9..c00098a5a7 100644
--- a/std/math/expo2.zig
+++ b/std/math/expo2.zig
@@ -1,5 +1,12 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/math/__expo2f.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/__expo2.c
+
 const math = @import("../math.zig");
 
+/// Returns exp(x) / 2 for x >= log(maxFloat(T)).
 pub fn expo2(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/fabs.zig b/std/math/fabs.zig
index e30f788ae7..6469f38835 100644
--- a/std/math/fabs.zig
+++ b/std/math/fabs.zig
@@ -1,13 +1,19 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - fabs(+-inf) = +inf
-// - fabs(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/fabsf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/fabs.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns the absolute value of x.
+///
+/// Special Cases:
+///  - fabs(+-inf) = +inf
+///  - fabs(nan)   = nan
 pub fn fabs(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/floor.zig b/std/math/floor.zig
index ab45a8fee7..e5ff2b1fc1 100644
--- a/std/math/floor.zig
+++ b/std/math/floor.zig
@@ -1,14 +1,20 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - floor(+-0)   = +-0
-// - floor(+-inf) = +-inf
-// - floor(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/floorf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/floor.c
 
 const builtin = @import("builtin");
 const expect = std.testing.expect;
 const std = @import("../std.zig");
 const math = std.math;
 
+/// Returns the greatest integer value less than or equal to x.
+///
+/// Special Cases:
+///  - floor(+-0)   = +-0
+///  - floor(+-inf) = +-inf
+///  - floor(nan)   = nan
 pub fn floor(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/fma.zig b/std/math/fma.zig
index a317bc96de..19c306fa2a 100644
--- a/std/math/fma.zig
+++ b/std/math/fma.zig
@@ -1,7 +1,14 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c
+
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns x * y + z with a single rounding error.
 pub fn fma(comptime T: type, x: T, y: T, z: T) T {
     return switch (T) {
         f32 => fma32(x, y, z),
@@ -16,7 +23,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {
     const u = @bitCast(u64, xy_z);
     const e = (u >> 52) & 0x7FF;
 
-    if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or xy_z - xy == z) {
+    if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or (xy_z - xy == z and xy_z - z == xy)) {
         return @floatCast(f32, xy_z);
     } else {
         // TODO: Handle inexact case with double-rounding
@@ -24,6 +31,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {
     }
 }
 
+// NOTE: Upstream fma.c has been rewritten completely to raise fp exceptions more accurately.
 fn fma64(x: f64, y: f64, z: f64) f64 {
     if (!math.isFinite(x) or !math.isFinite(y)) {
         return x * y + z;
diff --git a/std/math/frexp.zig b/std/math/frexp.zig
index 35eec315d9..2759cd6492 100644
--- a/std/math/frexp.zig
+++ b/std/math/frexp.zig
@@ -1,8 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - frexp(+-0)   = +-0, 0
-// - frexp(+-inf) = +-inf, 0
-// - frexp(nan)   = nan, undefined
+// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/frexp.c
 
 const std = @import("../std.zig");
 const math = std.math;
@@ -17,6 +17,13 @@ fn frexp_result(comptime T: type) type {
 pub const frexp32_result = frexp_result(f32);
 pub const frexp64_result = frexp_result(f64);
 
+/// Breaks x into a normalized fraction and an integral power of two.
+/// f == frac * 2^exp, with |frac| in the interval [0.5, 1).
+///
+/// Special Cases:
+///  - frexp(+-0)   = +-0, 0
+///  - frexp(+-inf) = +-inf, 0
+///  - frexp(nan)   = nan, undefined
 pub fn frexp(x: var) frexp_result(@typeOf(x)) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/hypot.zig b/std/math/hypot.zig
index fddbe7c068..c15da1495e 100644
--- a/std/math/hypot.zig
+++ b/std/math/hypot.zig
@@ -1,15 +1,21 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - hypot(+-inf, y)  = +inf
-// - hypot(x, +-inf)  = +inf
-// - hypot(nan, y)    = nan
-// - hypot(x, nan)    = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/hypotf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/hypot.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns sqrt(x * x + y * y), avoiding unncessary overflow and underflow.
+///
+/// Special Cases:
+///  - hypot(+-inf, y)  = +inf
+///  - hypot(x, +-inf)  = +inf
+///  - hypot(nan, y)    = nan
+///  - hypot(x, nan)    = nan
 pub fn hypot(comptime T: type, x: T, y: T) T {
     return switch (T) {
         f32 => hypot32(x, y),
diff --git a/std/math/ilogb.zig b/std/math/ilogb.zig
index 2dfd42a62c..fe4158a6dd 100644
--- a/std/math/ilogb.zig
+++ b/std/math/ilogb.zig
@@ -1,8 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - ilogb(+-inf) = maxInt(i32)
-// - ilogb(0)     = maxInt(i32)
-// - ilogb(nan)   = maxInt(i32)
+// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogbf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogb.c
 
 const std = @import("../std.zig");
 const math = std.math;
@@ -10,6 +10,12 @@ const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 const minInt = std.math.minInt;
 
+/// Returns the binary exponent of x as an integer.
+///
+/// Special Cases:
+///  - ilogb(+-inf) = maxInt(i32)
+///  - ilogb(0)     = maxInt(i32)
+///  - ilogb(nan)   = maxInt(i32)
 pub fn ilogb(x: var) i32 {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/inf.zig b/std/math/inf.zig
index e1bfbb197a..86ff245533 100644
--- a/std/math/inf.zig
+++ b/std/math/inf.zig
@@ -1,6 +1,7 @@
 const std = @import("../std.zig");
 const math = std.math;
 
+/// Returns value inf for the type T.
 pub fn inf(comptime T: type) T {
     return switch (T) {
         f16 => math.inf_f16,
diff --git a/std/math/isfinite.zig b/std/math/isfinite.zig
index ee8a5ff590..99eba668f9 100644
--- a/std/math/isfinite.zig
+++ b/std/math/isfinite.zig
@@ -3,6 +3,7 @@ const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns whether x is a finite value.
 pub fn isFinite(x: var) bool {
     const T = @typeOf(x);
     switch (T) {
diff --git a/std/math/isinf.zig b/std/math/isinf.zig
index 1b1759bd36..37934f4cf4 100644
--- a/std/math/isinf.zig
+++ b/std/math/isinf.zig
@@ -3,6 +3,7 @@ const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns whether x is an infinity, ignoring sign.
 pub fn isInf(x: var) bool {
     const T = @typeOf(x);
     switch (T) {
@@ -28,6 +29,7 @@ pub fn isInf(x: var) bool {
     }
 }
 
+/// Returns whether x is an infinity with a positive sign.
 pub fn isPositiveInf(x: var) bool {
     const T = @typeOf(x);
     switch (T) {
@@ -49,6 +51,7 @@ pub fn isPositiveInf(x: var) bool {
     }
 }
 
+/// Returns whether x is an infinity with a negative sign.
 pub fn isNegativeInf(x: var) bool {
     const T = @typeOf(x);
     switch (T) {
diff --git a/std/math/isnan.zig b/std/math/isnan.zig
index 9e541bf0a2..cf8cd2e1c2 100644
--- a/std/math/isnan.zig
+++ b/std/math/isnan.zig
@@ -3,13 +3,15 @@ const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns whether x is a nan.
 pub fn isNan(x: var) bool {
     return x != x;
 }
 
-/// Note: A signalling nan is identical to a standard nan right now but may have a different bit
-/// representation in the future when required.
+/// Returns whether x is a signalling nan.
 pub fn isSignalNan(x: var) bool {
+    // Note: A signalling nan is identical to a standard nan right now but may have a different bit
+    // representation in the future when required.
     return isNan(x);
 }
 
diff --git a/std/math/isnormal.zig b/std/math/isnormal.zig
index cddcada1d3..f8611ef805 100644
--- a/std/math/isnormal.zig
+++ b/std/math/isnormal.zig
@@ -3,6 +3,7 @@ const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+// Returns whether x has a normalized representation (i.e. integer part of mantissa is 1).
 pub fn isNormal(x: var) bool {
     const T = @typeOf(x);
     switch (T) {
diff --git a/std/math/ln.zig b/std/math/ln.zig
index 82b212f00f..c5d4c9ff25 100644
--- a/std/math/ln.zig
+++ b/std/math/ln.zig
@@ -1,9 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - ln(+inf)  = +inf
-// - ln(0)     = -inf
-// - ln(x)     = nan if x < 0
-// - ln(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/lnf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/ln.c
 
 const std = @import("../std.zig");
 const math = std.math;
@@ -11,6 +10,13 @@ const expect = std.testing.expect;
 const builtin = @import("builtin");
 const TypeId = builtin.TypeId;
 
+/// Returns the natural logarithm of x.
+///
+/// Special Cases:
+///  - ln(+inf)  = +inf
+///  - ln(0)     = -inf
+///  - ln(x)     = nan if x < 0
+///  - ln(nan)   = nan
 pub fn ln(x: var) @typeOf(x) {
     const T = @typeOf(x);
     switch (@typeId(T)) {
diff --git a/std/math/log.zig b/std/math/log.zig
index 100dc2fb7f..77f3639fd2 100644
--- a/std/math/log.zig
+++ b/std/math/log.zig
@@ -1,9 +1,16 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/math/logf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/log.c
+
 const std = @import("../std.zig");
 const math = std.math;
 const builtin = @import("builtin");
 const TypeId = builtin.TypeId;
 const expect = std.testing.expect;
 
+/// Returns the logarithm of x for the provided base.
 pub fn log(comptime T: type, base: T, x: T) T {
     if (base == 2) {
         return math.log2(x);
diff --git a/std/math/log10.zig b/std/math/log10.zig
index 7311778ddd..9b0bc3ac52 100644
--- a/std/math/log10.zig
+++ b/std/math/log10.zig
@@ -1,9 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - log10(+inf)  = +inf
-// - log10(0)     = -inf
-// - log10(x)     = nan if x < 0
-// - log10(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/log10f.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/log10.c
 
 const std = @import("../std.zig");
 const math = std.math;
@@ -12,6 +11,13 @@ const builtin = @import("builtin");
 const TypeId = builtin.TypeId;
 const maxInt = std.math.maxInt;
 
+/// Returns the base-10 logarithm of x.
+///
+/// Special Cases:
+///  - log10(+inf)  = +inf
+///  - log10(0)     = -inf
+///  - log10(x)     = nan if x < 0
+///  - log10(nan)   = nan
 pub fn log10(x: var) @typeOf(x) {
     const T = @typeOf(x);
     switch (@typeId(T)) {
diff --git a/std/math/log1p.zig b/std/math/log1p.zig
index c9be1132be..bae6deb536 100644
--- a/std/math/log1p.zig
+++ b/std/math/log1p.zig
@@ -1,16 +1,22 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - log1p(+inf)  = +inf
-// - log1p(+-0)   = +-0
-// - log1p(-1)    = -inf
-// - log1p(x)     = nan if x < -1
-// - log1p(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/log1pf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/log1p.c
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the natural logarithm of 1 + x with greater accuracy when x is near zero.
+///
+/// Special Cases:
+///  - log1p(+inf)  = +inf
+///  - log1p(+-0)   = +-0
+///  - log1p(-1)    = -inf
+///  - log1p(x)     = nan if x < -1
+///  - log1p(nan)   = nan
 pub fn log1p(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/log2.zig b/std/math/log2.zig
index e2dbf4105a..88450a7ffd 100644
--- a/std/math/log2.zig
+++ b/std/math/log2.zig
@@ -1,9 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - log2(+inf)  = +inf
-// - log2(0)     = -inf
-// - log2(x)     = nan if x < 0
-// - log2(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/log2f.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/log2.c
 
 const std = @import("../std.zig");
 const math = std.math;
@@ -12,6 +11,13 @@ const builtin = @import("builtin");
 const TypeId = builtin.TypeId;
 const maxInt = std.math.maxInt;
 
+/// Returns the base-2 logarithm of x.
+///
+/// Special Cases:
+///  - log2(+inf)  = +inf
+///  - log2(0)     = -inf
+///  - log2(x)     = nan if x < 0
+///  - log2(nan)   = nan
 pub fn log2(x: var) @typeOf(x) {
     const T = @typeOf(x);
     switch (@typeId(T)) {
diff --git a/std/math/modf.zig b/std/math/modf.zig
index e5a4964e63..92194d4c75 100644
--- a/std/math/modf.zig
+++ b/std/math/modf.zig
@@ -1,7 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - modf(+-inf) = +-inf, nan
-// - modf(nan)   = nan, nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/modff.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/modf.c
 
 const std = @import("../std.zig");
 const math = std.math;
@@ -17,6 +18,12 @@ fn modf_result(comptime T: type) type {
 pub const modf32_result = modf_result(f32);
 pub const modf64_result = modf_result(f64);
 
+/// Returns the integer and fractional floating-point numbers that sum to x. The sign of each
+/// result is the same as the sign of x.
+///
+/// Special Cases:
+///  - modf(+-inf) = +-inf, nan
+///  - modf(nan)   = nan, nan
 pub fn modf(x: var) modf_result(@typeOf(x)) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/nan.zig b/std/math/nan.zig
index 1a11a141d5..5a01a5b3bd 100644
--- a/std/math/nan.zig
+++ b/std/math/nan.zig
@@ -1,5 +1,6 @@
 const math = @import("../math.zig");
 
+/// Returns the nan representation for type T.
 pub fn nan(comptime T: type) T {
     return switch (T) {
         f16 => math.nan_f16,
@@ -10,9 +11,10 @@ pub fn nan(comptime T: type) T {
     };
 }
 
-// Note: A signalling nan is identical to a standard right now by may have a different bit
-// representation in the future when required.
+/// Returns the signalling nan representation for type T.
 pub fn snan(comptime T: type) T {
+    // Note: A signalling nan is identical to a standard right now by may have a different bit
+    // representation in the future when required.
     return switch (T) {
         f16 => @bitCast(f16, math.nan_u16),
         f32 => @bitCast(f32, math.nan_u32),
diff --git a/std/math/pow.zig b/std/math/pow.zig
index 81bb2d95d5..18c9f80634 100644
--- a/std/math/pow.zig
+++ b/std/math/pow.zig
@@ -1,32 +1,36 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
 //
-//  pow(x, +-0)    = 1 for any x
-//  pow(1, y)      = 1 for any y
-//  pow(x, 1)      = x for any x
-//  pow(nan, y)    = nan
-//  pow(x, nan)    = nan
-//  pow(+-0, y)    = +-inf for y an odd integer < 0
-//  pow(+-0, -inf) = +inf
-//  pow(+-0, +inf) = +0
-//  pow(+-0, y)    = +inf for finite y < 0 and not an odd integer
-//  pow(+-0, y)    = +-0 for y an odd integer > 0
-//  pow(+-0, y)    = +0 for finite y > 0 and not an odd integer
-//  pow(-1, +-inf) = 1
-//  pow(x, +inf)   = +inf for |x| > 1
-//  pow(x, -inf)   = +0 for |x| > 1
-//  pow(x, +inf)   = +0 for |x| < 1
-//  pow(x, -inf)   = +inf for |x| < 1
-//  pow(+inf, y)   = +inf for y > 0
-//  pow(+inf, y)   = +0 for y < 0
-//  pow(-inf, y)   = pow(-0, -y)
-//  pow(x, y)      = nan for finite x < 0 and finite non-integer y
+// https://golang.org/src/math/pow.go
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
-// This implementation is taken from the go stlib, musl is a bit more complex.
+/// Returns x raised to the power of y (x^y).
+///
+/// Special Cases:
+///  - pow(x, +-0)    = 1 for any x
+///  - pow(1, y)      = 1 for any y
+///  - pow(x, 1)      = x for any x
+///  - pow(nan, y)    = nan
+///  - pow(x, nan)    = nan
+///  - pow(+-0, y)    = +-inf for y an odd integer < 0
+///  - pow(+-0, -inf) = +inf
+///  - pow(+-0, +inf) = +0
+///  - pow(+-0, y)    = +inf for finite y < 0 and not an odd integer
+///  - pow(+-0, y)    = +-0 for y an odd integer > 0
+///  - pow(+-0, y)    = +0 for finite y > 0 and not an odd integer
+///  - pow(-1, +-inf) = 1
+///  - pow(x, +inf)   = +inf for |x| > 1
+///  - pow(x, -inf)   = +0 for |x| > 1
+///  - pow(x, +inf)   = +0 for |x| < 1
+///  - pow(x, -inf)   = +inf for |x| < 1
+///  - pow(+inf, y)   = +inf for y > 0
+///  - pow(+inf, y)   = +0 for y < 0
+///  - pow(-inf, y)   = pow(-0, -y)
+///  - pow(x, y)      = nan for finite x < 0 and finite non-integer y
 pub fn pow(comptime T: type, x: T, y: T) T {
     if (@typeInfo(T) == builtin.TypeId.Int) {
         return math.powi(T, x, y) catch unreachable;
@@ -53,15 +57,6 @@ pub fn pow(comptime T: type, x: T, y: T) T {
         return x;
     }
 
-    // special case sqrt
-    if (y == 0.5) {
-        return math.sqrt(x);
-    }
-
-    if (y == -0.5) {
-        return 1 / math.sqrt(x);
-    }
-
     if (x == 0) {
         if (y < 0) {
             // pow(+-0, y) = +- 0   for y an odd integer
@@ -112,14 +107,16 @@ pub fn pow(comptime T: type, x: T, y: T) T {
         }
     }
 
-    var ay = y;
-    var flip = false;
-    if (ay < 0) {
-        ay = -ay;
-        flip = true;
+    // special case sqrt
+    if (y == 0.5) {
+        return math.sqrt(x);
+    }
+
+    if (y == -0.5) {
+        return 1 / math.sqrt(x);
     }
 
-    const r1 = math.modf(ay);
+    const r1 = math.modf(math.fabs(y));
     var yi = r1.ipart;
     var yf = r1.fpart;
 
@@ -148,8 +145,18 @@ pub fn pow(comptime T: type, x: T, y: T) T {
     var xe = r2.exponent;
     var x1 = r2.significand;
 
-    var i = @floatToInt(i32, yi);
+    var i = @floatToInt(@IntType(true, T.bit_count), yi);
     while (i != 0) : (i >>= 1) {
+        const overflow_shift = math.floatExponentBits(T) + 1;
+        if (xe < -(1 << overflow_shift) or (1 << overflow_shift) < xe) {
+            // catch xe before it overflows the left shift below
+            // Since i != 0 it has at least one bit still set, so ae will accumulate xe
+            // on at least one more iteration, ae += xe is a lower bound on ae
+            // the lower bound on ae exceeds the size of a float exp
+            // so the final call to Ldexp will produce under/overflow (0/Inf)
+            ae += xe;
+            break;
+        }
         if (i & 1 == 1) {
             a1 *= x1;
             ae += xe;
@@ -163,7 +170,7 @@ pub fn pow(comptime T: type, x: T, y: T) T {
     }
 
     // a *= a1 * 2^ae
-    if (flip) {
+    if (y < 0) {
         a1 = 1 / a1;
         ae = -ae;
     }
@@ -202,6 +209,9 @@ test "math.pow.special" {
     expect(pow(f32, 45, 1.0) == 45);
     expect(pow(f32, -45, 1.0) == -45);
     expect(math.isNan(pow(f32, math.nan(f32), 5.0)));
+    expect(math.isPositiveInf(pow(f32, -math.inf(f32), 0.5)));
+    expect(math.isPositiveInf(pow(f32, -0, -0.5)));
+    expect(pow(f32, -0, 0.5) == 0);
     expect(math.isNan(pow(f32, 5.0, math.nan(f32))));
     expect(math.isPositiveInf(pow(f32, 0.0, -1.0)));
     //expect(math.isNegativeInf(pow(f32, -0.0, -3.0))); TODO is this required?
@@ -232,3 +242,11 @@ test "math.pow.special" {
     expect(math.isNan(pow(f32, -1.0, 1.2)));
     expect(math.isNan(pow(f32, -12.4, 78.5)));
 }
+
+test "math.pow.overflow" {
+    expect(math.isPositiveInf(pow(f64, 2, 1 << 32)));
+    expect(pow(f64, 2, -(1 << 32)) == 0);
+    expect(math.isNegativeInf(pow(f64, -2, (1 << 32) + 1)));
+    expect(pow(f64, 0.5, 1 << 45) == 0);
+    expect(math.isPositiveInf(pow(f64, 0.5, -(1 << 45))));
+}
diff --git a/std/math/powi.zig b/std/math/powi.zig
index b7212efcbf..d80700e5cd 100644
--- a/std/math/powi.zig
+++ b/std/math/powi.zig
@@ -1,12 +1,7 @@
-// Special Cases:
+// Based on Rust, which is licensed under the MIT license.
+// https://github.com/rust-lang/rust/blob/360432f1e8794de58cd94f34c9c17ad65871e5b5/LICENSE-MIT
 //
-//  powi(x, +-0)   = 1 for any x
-//  powi(0, y)     = 0 for any y
-//  powi(1, y)     = 1 for any y
-//  powi(-1, y)    = -1 for for y an odd integer
-//  powi(-1, y)    = 1 for for y an even integer
-//  powi(x, y)     = Overflow for for y >= @sizeOf(x) - 1 y > 0
-//  powi(x, y)     = Underflow for for y > @sizeOf(x) - 1 y < 0
+// https://github.com/rust-lang/rust/blob/360432f1e8794de58cd94f34c9c17ad65871e5b5/src/libcore/num/mod.rs#L3423
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
@@ -14,7 +9,16 @@ const math = std.math;
 const assert = std.debug.assert;
 const testing = std.testing;
 
-// This implementation is based on that from the rust stlib
+/// Returns the power of x raised by the integer y (x^y).
+///
+/// Special Cases:
+///  - powi(x, +-0)   = 1 for any x
+///  - powi(0, y)     = 0 for any y
+///  - powi(1, y)     = 1 for any y
+///  - powi(-1, y)    = -1 for y an odd integer
+///  - powi(-1, y)    = 1 for y an even integer
+///  - powi(x, y)     = Overflow for y >= @sizeOf(x) - 1 or y > 0
+///  - powi(x, y)     = Underflow for y > @sizeOf(x) - 1 or y < 0
 pub fn powi(comptime T: type, x: T, y: T) (error{
     Overflow,
     Underflow,
diff --git a/std/math/round.zig b/std/math/round.zig
index 39ff56ca79..0b80a46ce5 100644
--- a/std/math/round.zig
+++ b/std/math/round.zig
@@ -1,14 +1,20 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - round(+-0)   = +-0
-// - round(+-inf) = +-inf
-// - round(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/roundf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/round.c
 
 const builtin = @import("builtin");
 const expect = std.testing.expect;
 const std = @import("../std.zig");
 const math = std.math;
 
+/// Returns x rounded to the nearest integer, rounding half away from zero.
+///
+/// Special Cases:
+///  - round(+-0)   = +-0
+///  - round(+-inf) = +-inf
+///  - round(nan)   = nan
 pub fn round(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/scalbn.zig b/std/math/scalbn.zig
index d1338f5acb..d5716d621c 100644
--- a/std/math/scalbn.zig
+++ b/std/math/scalbn.zig
@@ -1,7 +1,14 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbnf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbn.c
+
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns x * 2^n.
 pub fn scalbn(x: var, n: i32) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/signbit.zig b/std/math/signbit.zig
index 9727152b07..e5c5909292 100644
--- a/std/math/signbit.zig
+++ b/std/math/signbit.zig
@@ -2,6 +2,7 @@ const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns whether x is negative or negative 0.
 pub fn signbit(x: var) bool {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/sin.zig b/std/math/sin.zig
index e25b8a292b..f21db4054e 100644
--- a/std/math/sin.zig
+++ b/std/math/sin.zig
@@ -1,19 +1,24 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
 //
-// - sin(+-0)   = +-0
-// - sin(+-inf) = nan
-// - sin(nan)   = nan
+// https://golang.org/src/math/sin.go
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the sine of the radian value x.
+///
+/// Special Cases:
+///  - sin(+-0)   = +-0
+///  - sin(+-inf) = nan
+///  - sin(nan)   = nan
 pub fn sin(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
-        f32 => sin32(x),
-        f64 => sin64(x),
+        f32 => sin_(T, x),
+        f64 => sin_(T, x),
         else => @compileError("sin not implemented for " ++ @typeName(T)),
     };
 }
@@ -34,83 +39,27 @@ const C3 = 2.48015872888517045348E-5;
 const C4 = -1.38888888888730564116E-3;
 const C5 = 4.16666666666665929218E-2;
 
-// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.
-//
-// This may have slight differences on some edge cases and may need to replaced if so.
-fn sin32(x_: f32) f32 {
-    const pi4a = 7.85398125648498535156e-1;
-    const pi4b = 3.77489470793079817668E-8;
-    const pi4c = 2.69515142907905952645E-15;
-    const m4pi = 1.273239544735162542821171882678754627704620361328125;
-
-    var x = x_;
-    if (x == 0 or math.isNan(x)) {
-        return x;
-    }
-    if (math.isInf(x)) {
-        return math.nan(f32);
-    }
-
-    var sign = false;
-    if (x < 0) {
-        x = -x;
-        sign = true;
-    }
-
-    var y = math.floor(x * m4pi);
-    var j = @floatToInt(i64, y);
-
-    if (j & 1 == 1) {
-        j += 1;
-        y += 1;
-    }
+const pi4a = 7.85398125648498535156e-1;
+const pi4b = 3.77489470793079817668E-8;
+const pi4c = 2.69515142907905952645E-15;
+const m4pi = 1.273239544735162542821171882678754627704620361328125;
 
-    j &= 7;
-    if (j > 3) {
-        j -= 4;
-        sign = !sign;
-    }
-
-    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
-    const w = z * z;
-
-    const r = r: {
-        if (j == 1 or j == 2) {
-            break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
-        } else {
-            break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
-        }
-    };
-
-    if (sign) {
-        return -r;
-    } else {
-        return r;
-    }
-}
-
-fn sin64(x_: f64) f64 {
-    const pi4a = 7.85398125648498535156e-1;
-    const pi4b = 3.77489470793079817668E-8;
-    const pi4c = 2.69515142907905952645E-15;
-    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+fn sin_(comptime T: type, x_: T) T {
+    const I = @IntType(true, T.bit_count);
 
     var x = x_;
     if (x == 0 or math.isNan(x)) {
         return x;
     }
     if (math.isInf(x)) {
-        return math.nan(f64);
+        return math.nan(T);
     }
 
-    var sign = false;
-    if (x < 0) {
-        x = -x;
-        sign = true;
-    }
+    var sign = x < 0;
+    x = math.fabs(x);
 
     var y = math.floor(x * m4pi);
-    var j = @floatToInt(i64, y);
+    var j = @floatToInt(I, y);
 
     if (j & 1 == 1) {
         j += 1;
@@ -126,61 +75,56 @@ fn sin64(x_: f64) f64 {
     const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
     const w = z * z;
 
-    const r = r: {
-        if (j == 1 or j == 2) {
-            break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
-        } else {
-            break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
-        }
-    };
+    const r = if (j == 1 or j == 2)
+        1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))))
+    else
+        z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
 
-    if (sign) {
-        return -r;
-    } else {
-        return r;
-    }
+    return if (sign) -r else r;
 }
 
 test "math.sin" {
-    expect(sin(f32(0.0)) == sin32(0.0));
-    expect(sin(f64(0.0)) == sin64(0.0));
+    expect(sin(f32(0.0)) == sin_(f32, 0.0));
+    expect(sin(f64(0.0)) == sin_(f64, 0.0));
     expect(comptime (math.sin(f64(2))) == math.sin(f64(2)));
 }
 
 test "math.sin32" {
     const epsilon = 0.000001;
 
-    expect(math.approxEq(f32, sin32(0.0), 0.0, epsilon));
-    expect(math.approxEq(f32, sin32(0.2), 0.198669, epsilon));
-    expect(math.approxEq(f32, sin32(0.8923), 0.778517, epsilon));
-    expect(math.approxEq(f32, sin32(1.5), 0.997495, epsilon));
-    expect(math.approxEq(f32, sin32(37.45), -0.246544, epsilon));
-    expect(math.approxEq(f32, sin32(89.123), 0.916166, epsilon));
+    expect(math.approxEq(f32, sin_(f32, 0.0), 0.0, epsilon));
+    expect(math.approxEq(f32, sin_(f32, 0.2), 0.198669, epsilon));
+    expect(math.approxEq(f32, sin_(f32, 0.8923), 0.778517, epsilon));
+    expect(math.approxEq(f32, sin_(f32, 1.5), 0.997495, epsilon));
+    expect(math.approxEq(f32, sin_(f32, -1.5), -0.997495, epsilon));
+    expect(math.approxEq(f32, sin_(f32, 37.45), -0.246544, epsilon));
+    expect(math.approxEq(f32, sin_(f32, 89.123), 0.916166, epsilon));
 }
 
 test "math.sin64" {
     const epsilon = 0.000001;
 
-    expect(math.approxEq(f64, sin64(0.0), 0.0, epsilon));
-    expect(math.approxEq(f64, sin64(0.2), 0.198669, epsilon));
-    expect(math.approxEq(f64, sin64(0.8923), 0.778517, epsilon));
-    expect(math.approxEq(f64, sin64(1.5), 0.997495, epsilon));
-    expect(math.approxEq(f64, sin64(37.45), -0.246543, epsilon));
-    expect(math.approxEq(f64, sin64(89.123), 0.916166, epsilon));
+    expect(math.approxEq(f64, sin_(f64, 0.0), 0.0, epsilon));
+    expect(math.approxEq(f64, sin_(f64, 0.2), 0.198669, epsilon));
+    expect(math.approxEq(f64, sin_(f64, 0.8923), 0.778517, epsilon));
+    expect(math.approxEq(f64, sin_(f64, 1.5), 0.997495, epsilon));
+    expect(math.approxEq(f64, sin_(f64, -1.5), -0.997495, epsilon));
+    expect(math.approxEq(f64, sin_(f64, 37.45), -0.246543, epsilon));
+    expect(math.approxEq(f64, sin_(f64, 89.123), 0.916166, epsilon));
 }
 
 test "math.sin32.special" {
-    expect(sin32(0.0) == 0.0);
-    expect(sin32(-0.0) == -0.0);
-    expect(math.isNan(sin32(math.inf(f32))));
-    expect(math.isNan(sin32(-math.inf(f32))));
-    expect(math.isNan(sin32(math.nan(f32))));
+    expect(sin_(f32, 0.0) == 0.0);
+    expect(sin_(f32, -0.0) == -0.0);
+    expect(math.isNan(sin_(f32, math.inf(f32))));
+    expect(math.isNan(sin_(f32, -math.inf(f32))));
+    expect(math.isNan(sin_(f32, math.nan(f32))));
 }
 
 test "math.sin64.special" {
-    expect(sin64(0.0) == 0.0);
-    expect(sin64(-0.0) == -0.0);
-    expect(math.isNan(sin64(math.inf(f64))));
-    expect(math.isNan(sin64(-math.inf(f64))));
-    expect(math.isNan(sin64(math.nan(f64))));
+    expect(sin_(f64, 0.0) == 0.0);
+    expect(sin_(f64, -0.0) == -0.0);
+    expect(math.isNan(sin_(f64, math.inf(f64))));
+    expect(math.isNan(sin_(f64, -math.inf(f64))));
+    expect(math.isNan(sin_(f64, math.nan(f64))));
 }
diff --git a/std/math/sinh.zig b/std/math/sinh.zig
index cf4363c4a9..73ee65ea6f 100644
--- a/std/math/sinh.zig
+++ b/std/math/sinh.zig
@@ -1,8 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - sinh(+-0)   = +-0
-// - sinh(+-inf) = +-inf
-// - sinh(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/sinhf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/sinh.c
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
@@ -11,6 +11,12 @@ const expect = std.testing.expect;
 const expo2 = @import("expo2.zig").expo2;
 const maxInt = std.math.maxInt;
 
+/// Returns the hyperbolic sine of x.
+///
+/// Special Cases:
+///  - sinh(+-0)   = +-0
+///  - sinh(+-inf) = +-inf
+///  - sinh(nan)   = nan
 pub fn sinh(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/sqrt.zig b/std/math/sqrt.zig
index 9062f598a1..30af5915d4 100644
--- a/std/math/sqrt.zig
+++ b/std/math/sqrt.zig
@@ -1,10 +1,3 @@
-// Special Cases:
-//
-// - sqrt(+inf)  = +inf
-// - sqrt(+-0)   = +-0
-// - sqrt(x)     = nan if x < 0
-// - sqrt(nan)   = nan
-
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
@@ -12,6 +5,13 @@ const builtin = @import("builtin");
 const TypeId = builtin.TypeId;
 const maxInt = std.math.maxInt;
 
+/// Returns the square root of x.
+///
+/// Special Cases:
+///  - sqrt(+inf)  = +inf
+///  - sqrt(+-0)   = +-0
+///  - sqrt(x)     = nan if x < 0
+///  - sqrt(nan)   = nan
 pub fn sqrt(x: var) (if (@typeId(@typeOf(x)) == TypeId.Int) @IntType(false, @typeOf(x).bit_count / 2) else @typeOf(x)) {
     const T = @typeOf(x);
     switch (@typeId(T)) {
diff --git a/std/math/tan.zig b/std/math/tan.zig
index fc11ebdef7..e8259ee7ad 100644
--- a/std/math/tan.zig
+++ b/std/math/tan.zig
@@ -1,19 +1,24 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
 //
-// - tan(+-0)   = +-0
-// - tan(+-inf) = nan
-// - tan(nan)   = nan
+// https://golang.org/src/math/tan.go
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 
+/// Returns the tangent of the radian value x.
+///
+/// Special Cases:
+///  - tan(+-0)   = +-0
+///  - tan(+-inf) = nan
+///  - tan(nan)   = nan
 pub fn tan(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
-        f32 => tan32(x),
-        f64 => tan64(x),
+        f32 => tan_(f32, x),
+        f64 => tan_(f64, x),
         else => @compileError("tan not implemented for " ++ @typeName(T)),
     };
 }
@@ -27,80 +32,27 @@ const Tq2 = -1.32089234440210967447E6;
 const Tq3 = 2.50083801823357915839E7;
 const Tq4 = -5.38695755929454629881E7;
 
-// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.
-//
-// This may have slight differences on some edge cases and may need to replaced if so.
-fn tan32(x_: f32) f32 {
-    const pi4a = 7.85398125648498535156e-1;
-    const pi4b = 3.77489470793079817668E-8;
-    const pi4c = 2.69515142907905952645E-15;
-    const m4pi = 1.273239544735162542821171882678754627704620361328125;
-
-    var x = x_;
-    if (x == 0 or math.isNan(x)) {
-        return x;
-    }
-    if (math.isInf(x)) {
-        return math.nan(f32);
-    }
-
-    var sign = false;
-    if (x < 0) {
-        x = -x;
-        sign = true;
-    }
-
-    var y = math.floor(x * m4pi);
-    var j = @floatToInt(i64, y);
-
-    if (j & 1 == 1) {
-        j += 1;
-        y += 1;
-    }
-
-    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
-    const w = z * z;
-
-    var r = r: {
-        if (w > 1e-14) {
-            break :r z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4));
-        } else {
-            break :r z;
-        }
-    };
-
-    if (j & 2 == 2) {
-        r = -1 / r;
-    }
-    if (sign) {
-        r = -r;
-    }
-
-    return r;
-}
+const pi4a = 7.85398125648498535156e-1;
+const pi4b = 3.77489470793079817668E-8;
+const pi4c = 2.69515142907905952645E-15;
+const m4pi = 1.273239544735162542821171882678754627704620361328125;
 
-fn tan64(x_: f64) f64 {
-    const pi4a = 7.85398125648498535156e-1;
-    const pi4b = 3.77489470793079817668E-8;
-    const pi4c = 2.69515142907905952645E-15;
-    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+fn tan_(comptime T: type, x_: T) T {
+    const I = @IntType(true, T.bit_count);
 
     var x = x_;
     if (x == 0 or math.isNan(x)) {
         return x;
     }
     if (math.isInf(x)) {
-        return math.nan(f64);
+        return math.nan(T);
     }
 
-    var sign = false;
-    if (x < 0) {
-        x = -x;
-        sign = true;
-    }
+    var sign = x < 0;
+    x = math.fabs(x);
 
     var y = math.floor(x * m4pi);
-    var j = @floatToInt(i64, y);
+    var j = @floatToInt(I, y);
 
     if (j & 1 == 1) {
         j += 1;
@@ -110,63 +62,57 @@ fn tan64(x_: f64) f64 {
     const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
     const w = z * z;
 
-    var r = r: {
-        if (w > 1e-14) {
-            break :r z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4));
-        } else {
-            break :r z;
-        }
-    };
+    var r = if (w > 1e-14)
+        z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4))
+    else
+        z;
 
     if (j & 2 == 2) {
         r = -1 / r;
     }
-    if (sign) {
-        r = -r;
-    }
 
-    return r;
+    return if (sign) -r else r;
 }
 
 test "math.tan" {
-    expect(tan(f32(0.0)) == tan32(0.0));
-    expect(tan(f64(0.0)) == tan64(0.0));
+    expect(tan(f32(0.0)) == tan_(f32, 0.0));
+    expect(tan(f64(0.0)) == tan_(f64, 0.0));
 }
 
 test "math.tan32" {
     const epsilon = 0.000001;
 
-    expect(math.approxEq(f32, tan32(0.0), 0.0, epsilon));
-    expect(math.approxEq(f32, tan32(0.2), 0.202710, epsilon));
-    expect(math.approxEq(f32, tan32(0.8923), 1.240422, epsilon));
-    expect(math.approxEq(f32, tan32(1.5), 14.101420, epsilon));
-    expect(math.approxEq(f32, tan32(37.45), -0.254397, epsilon));
-    expect(math.approxEq(f32, tan32(89.123), 2.285852, epsilon));
+    expect(math.approxEq(f32, tan_(f32, 0.0), 0.0, epsilon));
+    expect(math.approxEq(f32, tan_(f32, 0.2), 0.202710, epsilon));
+    expect(math.approxEq(f32, tan_(f32, 0.8923), 1.240422, epsilon));
+    expect(math.approxEq(f32, tan_(f32, 1.5), 14.101420, epsilon));
+    expect(math.approxEq(f32, tan_(f32, 37.45), -0.254397, epsilon));
+    expect(math.approxEq(f32, tan_(f32, 89.123), 2.285852, epsilon));
 }
 
 test "math.tan64" {
     const epsilon = 0.000001;
 
-    expect(math.approxEq(f64, tan64(0.0), 0.0, epsilon));
-    expect(math.approxEq(f64, tan64(0.2), 0.202710, epsilon));
-    expect(math.approxEq(f64, tan64(0.8923), 1.240422, epsilon));
-    expect(math.approxEq(f64, tan64(1.5), 14.101420, epsilon));
-    expect(math.approxEq(f64, tan64(37.45), -0.254397, epsilon));
-    expect(math.approxEq(f64, tan64(89.123), 2.2858376, epsilon));
+    expect(math.approxEq(f64, tan_(f64, 0.0), 0.0, epsilon));
+    expect(math.approxEq(f64, tan_(f64, 0.2), 0.202710, epsilon));
+    expect(math.approxEq(f64, tan_(f64, 0.8923), 1.240422, epsilon));
+    expect(math.approxEq(f64, tan_(f64, 1.5), 14.101420, epsilon));
+    expect(math.approxEq(f64, tan_(f64, 37.45), -0.254397, epsilon));
+    expect(math.approxEq(f64, tan_(f64, 89.123), 2.2858376, epsilon));
 }
 
 test "math.tan32.special" {
-    expect(tan32(0.0) == 0.0);
-    expect(tan32(-0.0) == -0.0);
-    expect(math.isNan(tan32(math.inf(f32))));
-    expect(math.isNan(tan32(-math.inf(f32))));
-    expect(math.isNan(tan32(math.nan(f32))));
+    expect(tan_(f32, 0.0) == 0.0);
+    expect(tan_(f32, -0.0) == -0.0);
+    expect(math.isNan(tan_(f32, math.inf(f32))));
+    expect(math.isNan(tan_(f32, -math.inf(f32))));
+    expect(math.isNan(tan_(f32, math.nan(f32))));
 }
 
 test "math.tan64.special" {
-    expect(tan64(0.0) == 0.0);
-    expect(tan64(-0.0) == -0.0);
-    expect(math.isNan(tan64(math.inf(f64))));
-    expect(math.isNan(tan64(-math.inf(f64))));
-    expect(math.isNan(tan64(math.nan(f64))));
+    expect(tan_(f64, 0.0) == 0.0);
+    expect(tan_(f64, -0.0) == -0.0);
+    expect(math.isNan(tan_(f64, math.inf(f64))));
+    expect(math.isNan(tan_(f64, -math.inf(f64))));
+    expect(math.isNan(tan_(f64, math.nan(f64))));
 }
diff --git a/std/math/tanh.zig b/std/math/tanh.zig
index f88ecfdc9c..48d26d091e 100644
--- a/std/math/tanh.zig
+++ b/std/math/tanh.zig
@@ -1,8 +1,8 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - sinh(+-0)   = +-0
-// - sinh(+-inf) = +-1
-// - sinh(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/tanhf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/tanh.c
 
 const builtin = @import("builtin");
 const std = @import("../std.zig");
@@ -11,6 +11,12 @@ const expect = std.testing.expect;
 const expo2 = @import("expo2.zig").expo2;
 const maxInt = std.math.maxInt;
 
+/// Returns the hyperbolic tangent of x.
+///
+/// Special Cases:
+///  - sinh(+-0)   = +-0
+///  - sinh(+-inf) = +-1
+///  - sinh(nan)   = nan
 pub fn tanh(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {
diff --git a/std/math/trunc.zig b/std/math/trunc.zig
index 9449d0e7eb..219bcd4914 100644
--- a/std/math/trunc.zig
+++ b/std/math/trunc.zig
@@ -1,14 +1,20 @@
-// Special Cases:
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
 //
-// - trunc(+-0)   = +-0
-// - trunc(+-inf) = +-inf
-// - trunc(nan)   = nan
+// https://git.musl-libc.org/cgit/musl/tree/src/math/truncf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/trunc.c
 
 const std = @import("../std.zig");
 const math = std.math;
 const expect = std.testing.expect;
 const maxInt = std.math.maxInt;
 
+/// Returns the integer value of x.
+///
+/// Special Cases:
+///  - trunc(+-0)   = +-0
+///  - trunc(+-inf) = +-inf
+///  - trunc(nan)   = nan
 pub fn trunc(x: var) @typeOf(x) {
     const T = @typeOf(x);
     return switch (T) {