# hare

The Hare programming language
git clone https://git.torresjrjr.com/hare.git

```commit e29ec221cd43655821a89857c8706fc41cb94421
parent 3c362bec22fb22ac81f7929365770cc612e0e29f
Date:   Tue, 17 Aug 2021 17:37:53 +0200

Add math.ha and clean up floats.ha

Diffstat:
Mmath/floats.ha | 457++++++++++++++++++++++++++++++++++++++++++++-----------------------------------
Amath/ints.ha | 138+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Amath/math.ha | 690+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Mscripts/gen-stdlib | 2+-
Mstdlib.mk | 12++++++++----
```
```5 files changed, 1090 insertions(+), 209 deletions(-)
diff --git a/math/floats.ha b/math/floats.ha
@@ -1,144 +1,5 @@
use types;

-// The floating point value representing Not a Number, i.e. an undefined or
-// unrepresentable value. You cannot test if a number is NaN by comparing to
-// this value; see [[isnan]] instead.
-export def NAN: f32 = 0.0 / 0.0;
-
-// The floating point value representing positive infinity. Use -[[INF]] for
-// negative infinity.
-export def INF: f32 = 1.0 / 0.0;
-
-// Returns true if the given f64 is NaN.
-export fn isnanf64(n: f64) bool = n != n;
-
-// Returns true if the given f32 is NaN.
-export fn isnanf32(n: f32) bool = n != n;
-
-// Returns true if the given floating-point number is NaN.
-export fn isnan(n: types::floating) bool = {
-	return match (n) {
-		f: f64 => isnanf64(f),
-		f: f32 => isnanf32(f),
-	};
-};
-
-@test fn isnan() void = {
-	assert(isnan(NAN));
-	assert(isnan(-NAN));
-	assert(isnan(f64frombits(0xfffabcdef1234567)));
-	assert(!isnan(INF));
-	assert(!isnan(1.23f32));
-};
-
-// Returns true if the given floating-point number is infinite.
-export fn isinf(n: f64) bool = {
-	const bits = f64bits(n);
-	const mant = bits & F64_MANTISSA_MASK;
-	const exp = bits >> F64_MANTISSA_BITS & F64_EXPONENT_MASK;
-	return exp == F64_EXPONENT_MASK && mant == 0;
-};
-
-@test fn isinf() void = {
-	assert(isinf(INF));
-	assert(isinf(-INF));
-	assert(!isinf(NAN));
-	assert(!isinf(1.23));
-	assert(!isinf(-1.23f32));
-};
-
-// Returns true if the given floating-point number is normal.
-export fn isnormal(n: types::floating) bool = {
-	return match (n) {
-		n: f32 => isnormalf32(n),
-		n: f64 => isnormalf64(n),
-	};
-};
-
-// Returns true if the given f32 is normal.
-export fn isnormalf32(n: f32) bool = {
-	const bits = f32bits(n);
-	const mant = bits & F32_MANTISSA_MASK;
-	const exp = bits >> F32_MANTISSA_BITS & F32_EXPONENT_MASK;
-	return exp != F32_EXPONENT_MASK && (exp > 0 || mant == 0);
-};
-
-// Returns true if the given f64 is normal.
-export fn isnormalf64(n: f64) bool = {
-	const bits = f64bits(n);
-	const mant = bits & F64_MANTISSA_MASK;
-	const exp = bits >> F64_MANTISSA_BITS & F64_EXPONENT_MASK;
-	return exp != F64_EXPONENT_MASK && (exp > 0 || mant == 0);
-};
-
-// Returns true if the given floating-point number is subnormal.
-export fn issubnormal(n: types::floating) bool = {
-	return match (n) {
-		n: f32 => issubnormalf32(n),
-		n: f64 => issubnormalf64(n),
-	};
-};
-
-// Returns true if the given f32 is subnormal.
-export fn issubnormalf32(n: f32) bool = {
-	const bits = f32bits(n);
-	const mant = bits & F32_MANTISSA_MASK;
-	const exp = bits >> F32_MANTISSA_BITS & F32_EXPONENT_MASK;
-	return exp == 0 && mant != 0;
-};
-
-// Returns true if the given f64 is subnormal.
-export fn issubnormalf64(n: f64) bool = {
-	const bits = f64bits(n);
-	const mant = bits & F64_MANTISSA_MASK;
-	const exp = bits >> F64_MANTISSA_BITS & F64_EXPONENT_MASK;
-	return exp == 0 && mant != 0;
-};
-
-@test fn float_normality() void = {
-	assert(isnormal(0.0));
-	assert(isnormal(1.0));
-	assert(!isnormal(NAN));
-	assert(!isnormal(INF));
-	assert(!isnormal(1.0e-310));
-	assert(!isnormal(1.0e-40f32));
-
-	assert(isnormalf32(1.0));
-	assert(isnormalf32(0.0));
-	assert(!isnormalf32(NAN));
-	assert(!isnormalf32(INF));
-	assert(!isnormalf32(-1.0e-40));
-	assert(isnormalf32(-1.0e-50));
-
-	assert(isnormalf64(1.0));
-	assert(isnormalf64(0.0));
-	assert(!isnormalf64(NAN));
-	assert(!isnormalf64(INF));
-	assert(!isnormalf64(-1.0e-320));
-	assert(isnormalf64(-1.0e-330));
-
-	assert(issubnormal(1.0e-320));
-	assert(issubnormal(1.0e-42f32));
-	assert(!issubnormal(NAN));
-	assert(!issubnormal(INF));
-	assert(!issubnormal(1.0));
-	assert(!issubnormal(0.0));
-
-	assert(issubnormalf32(1.0e-45));
-	assert(issubnormalf32(-1.0e-39));
-	assert(!issubnormalf32(-NAN));
-	assert(!issubnormalf32(-INF));
-	assert(!issubnormalf32(0.0));
-	assert(!issubnormalf32(-1.0e-49));
-
-	assert(issubnormalf64(5.0e-324));
-	assert(issubnormalf64(-2.0e-310));
-	assert(!issubnormalf64(-NAN));
-	assert(!issubnormalf64(-INF));
-	assert(!issubnormalf64(-1.0e-400));
-	assert(!issubnormalf64(0.0));
-};
-
// Returns the binary representation of the given f64.
export fn f64bits(n: f64) u64 = *(&n: *u64);

@@ -235,7 +96,8 @@ def F32_SIGN_MASK: u32 = 1u32 << 31;

// The f32 that contains only an exponent that evaluates to zero.
-def F32_EXP_ZERO: u32 = ((F32_EXPONENT_BIAS: u32) - 1) << (F32_MANTISSA_BITS: u32);
+def F32_EXP_ZERO: u32 =
+	((F32_EXPONENT_BIAS: u32) - 1) << (F32_MANTISSA_BITS: u32);

// Contains information about the structure of a specific floating point number
// type.
@@ -270,36 +132,179 @@ export const f32info: floatinfo = floatinfo {
expmask = (1 << 8) - 1,
};

-// Returns the absolute value of n.
+// The floating point value representing Not a Number, i.e. an undefined or
+// unrepresentable value. You cannot test if a number is NaN by comparing to
+// this value; see [[isnan]] instead.
+export def NAN: f32 = 0.0 / 0.0;
+
+// The floating point value representing positive infinity. Use -[[INF]] for
+// negative infinity.
+export def INF: f32 = 1.0 / 0.0;
+
+// Returns true if the given floating-point number is NaN.
+export fn isnan(n: f64) bool = n != n;
+
+@test fn isnan() void = {
+	assert(isnan(NAN));
+	assert(isnan(-NAN));
+	assert(isnan(f64frombits(0xfffabcdef1234567)));
+	assert(!isnan(INF));
+	assert(!isnan(1.23f32));
+};
+
+// Returns true if the given floating-point number is infinite.
+export fn isinf(n: f64) bool = {
+	const bits = f64bits(n);
+	const mant = bits & F64_MANTISSA_MASK;
+	const exp = bits >> F64_MANTISSA_BITS & F64_EXPONENT_MASK;
+	return exp == F64_EXPONENT_MASK && mant == 0;
+};
+
+// Returns true if the given floating-point number is positive infinity.
+export fn isposinf(n: f64) bool = {
+	return isinf(n) && signf64(n) == 1i64;
+};
+
+// Returns true if the given floating-point number is negative infinity.
+export fn isneginf(n: f64) bool = {
+	return isinf(n) && signf64(n) == -1i64;
+};
+
+@test fn isinf() void = {
+	assert(isinf(INF));
+	assert(isinf(-INF));
+	assert(!isinf(NAN));
+	assert(!isinf(1.23));
+	assert(!isinf(-1.23f32));
+	assert(isposinf(math::INF));
+	assert(!isposinf(-math::INF));
+	assert(isneginf(-math::INF));
+	assert(!isneginf(math::INF));
+};
+
+// Returns true if the given floating-point number is normal.
+export fn isnormal(n: types::floating) bool = {
+	return match (n) {
+		n: f32 => isnormalf32(n),
+		n: f64 => isnormalf64(n),
+	};
+};
+
+// Returns true if the given f64 is normal.
+export fn isnormalf64(n: f64) bool = {
+	const bits = f64bits(n);
+	const mant = bits & F64_MANTISSA_MASK;
+	const exp = bits >> F64_MANTISSA_BITS & F64_EXPONENT_MASK;
+	return exp != F64_EXPONENT_MASK && (exp > 0 || mant == 0);
+};
+
+// Returns true if the given f32 is normal.
+export fn isnormalf32(n: f32) bool = {
+	const bits = f32bits(n);
+	const mant = bits & F32_MANTISSA_MASK;
+	const exp = bits >> F32_MANTISSA_BITS & F32_EXPONENT_MASK;
+	return exp != F32_EXPONENT_MASK && (exp > 0 || mant == 0);
+};
+
+// Returns true if the given floating-point number is subnormal.
+export fn issubnormal(n: types::floating) bool = {
+	return match (n) {
+		n: f32 => issubnormalf32(n),
+		n: f64 => issubnormalf64(n),
+	};
+};
+
+// Returns true if the given f64 is subnormal.
+export fn issubnormalf64(n: f64) bool = {
+	const bits = f64bits(n);
+	const mant = bits & F64_MANTISSA_MASK;
+	const exp = bits >> F64_MANTISSA_BITS & F64_EXPONENT_MASK;
+	return exp == 0 && mant != 0;
+};
+
+// Returns true if the given f32 is subnormal.
+export fn issubnormalf32(n: f32) bool = {
+	const bits = f32bits(n);
+	const mant = bits & F32_MANTISSA_MASK;
+	const exp = bits >> F32_MANTISSA_BITS & F32_EXPONENT_MASK;
+	return exp == 0 && mant != 0;
+};
+
+@test fn float_normality() void = {
+	assert(isnormal(0.0));
+	assert(isnormal(1.0));
+	assert(!isnormal(NAN));
+	assert(!isnormal(INF));
+	assert(!isnormal(1.0e-310));
+	assert(!isnormal(1.0e-40f32));
+
+	assert(isnormalf32(1.0));
+	assert(isnormalf32(0.0));
+	assert(!isnormalf32(NAN));
+	assert(!isnormalf32(INF));
+	assert(!isnormalf32(-1.0e-40));
+	assert(isnormalf32(-1.0e-50));
+
+	assert(isnormalf64(1.0));
+	assert(isnormalf64(0.0));
+	assert(!isnormalf64(NAN));
+	assert(!isnormalf64(INF));
+	assert(!isnormalf64(-1.0e-320));
+	assert(isnormalf64(-1.0e-330));
+
+	assert(issubnormal(1.0e-320));
+	assert(issubnormal(1.0e-42f32));
+	assert(!issubnormal(NAN));
+	assert(!issubnormal(INF));
+	assert(!issubnormal(1.0));
+	assert(!issubnormal(0.0));
+
+	assert(issubnormalf32(1.0e-45));
+	assert(issubnormalf32(-1.0e-39));
+	assert(!issubnormalf32(-NAN));
+	assert(!issubnormalf32(-INF));
+	assert(!issubnormalf32(0.0));
+	assert(!issubnormalf32(-1.0e-49));
+
+	assert(issubnormalf64(5.0e-324));
+	assert(issubnormalf64(-2.0e-310));
+	assert(!issubnormalf64(-NAN));
+	assert(!issubnormalf64(-INF));
+	assert(!issubnormalf64(-1.0e-400));
+	assert(!issubnormalf64(0.0));
+};
+
+// Returns the absolute value of f64 n.
export fn absf64(n: f64) f64 = {
};

-// Returns the absolute value of n.
+// Returns the absolute value of f32 n.
export fn absf32(n: f32) f32 = {
};

-// Returns the absolute value of n.
-export fn abs(n: types::floating) f64 = {
+// Returns the absolute value of floating-point number n.
+export fn absf(n: types::floating) f64 = {
return match (n) {
n: f64 => absf64(n),
n: f32 => (absf32(n): f64),
};
};

-@test fn test_abs() void = {
+@test fn absf() void = {
assert(absf64(2.0f64) == 2.0f64);
assert(absf32(2.0f32) == 2.0f32);
-	assert(abs(2.0f64) == 2.0f64);
-	assert(abs(2.0f32) == 2.0f64);
-	assert(abs(-2.0f64) == 2.0f64);
-	assert(abs(-2.0f32) == 2.0f32);
-	assert(abs(0f32) == 0f32);
-	assert(abs(0f64) == 0f64);
+	assert(absf(2.0f64) == 2.0f64);
+	assert(absf(2.0f32) == 2.0f64);
+	assert(absf(-2.0f64) == 2.0f64);
+	assert(absf(-2.0f32) == 2.0f32);
+	assert(absf(0f32) == 0f32);
+	assert(absf(0f64) == 0f64);
};

-// Returns 1 if x is positive and -1 if x is negative. Note that zero is also signed.
+// Returns 1 if x is positive and -1 if x is negative. Note that zero is also
+// signed.
export fn signf64(x: f64) i64 = {
if (f64bits(x) & F64_SIGN_MASK == 0) {
return 1i64;
@@ -308,7 +313,8 @@ export fn signf64(x: f64) i64 = {
};
};

-// Returns 1 if x is positive and -1 if x is negative. Note that zero is also signed.
+// Returns 1 if x is positive and -1 if x is negative. Note that zero is also
+// signed.
export fn signf32(x: f32) i64 = {
if (f32bits(x) & F32_SIGN_MASK == 0) {
return 1i64;
@@ -317,21 +323,54 @@ export fn signf32(x: f32) i64 = {
};
};

-// Returns 1 if x is positive and -1 if x is negative. Note that zero is also signed.
-export fn sign(x: types::floating) i64 = {
+// Returns 1 if x is positive and -1 if x is negative. Note that zero is also
+// signed.
+export fn signf(x: types::floating) i64 = {
return match (x) {
n: f64 => signf64(n),
n: f32 => signf32(n),
};
};

-@test fn test_sign() void = {
-	assert(sign(0.0f64) == 1i64);
-	assert(sign(-0.0f64) == -1i64);
-	assert(sign(0.0f32) == 1i64);
-	assert(sign(-0.0f32) == -1i64);
-	assert(sign(1.5f64) == 1i64);
-	assert(sign(-1.5f64) == -1i64);
+@test fn signf() void = {
+	assert(signf(0.0f64) == 1i64);
+	assert(signf(-0.0f64) == -1i64);
+	assert(signf(0.0f32) == 1i64);
+	assert(signf(-0.0f32) == -1i64);
+	assert(signf(1.5f64) == 1i64);
+	assert(signf(-1.5f64) == -1i64);
+	assert(ispositive(1.0f64));
+	assert(!ispositive(-1.0f64));
+	assert(isnegative(-1.0f64));
+	assert(!isnegative(1.0f64));
+};
+
+// Returns whether or not x is positive.
+export fn ispositivef64(x: f64) bool = signf64(x) == 1i64;
+
+// Returns whether or not x is positive.
+export fn ispositivef32(x: f32) bool = signf32(x) == 1i32;
+
+// Returns whether or not x is positive.
+export fn ispositive(x: types::floating) bool = {
+	return match (x) {
+		n: f64 => ispositivef64(n),
+		n: f32 => ispositivef32(n),
+	};
+};
+
+// Returns whether or not x is negative.
+export fn isnegativef64(x: f64) bool = signf64(x) == -1i64;
+
+// Returns whether or not x is negative.
+export fn isnegativef32(x: f32) bool = signf32(x) == -1i32;
+
+// Returns whether or not x is negative.
+export fn isnegative(x: types::floating) bool = {
+	return match (x) {
+		n: f64 => isnegativef64(n),
+		n: f32 => isnegativef32(n),
+	};
};

// Returns x, but with the sign of y.
@@ -354,7 +393,7 @@ export fn copysign(x: types::floating, y: types::floating) f64 = {
};
};

-@test fn test_copysign() void = {
+@test fn copysign() void = {
assert(copysign(100.0f64, 1.0f64) == 100.0f64);
assert(copysign(100.0f64, -1.0f64) == -100.0f64);
assert(copysign(100.0f32, 1.0f32) == 100.0f32);
@@ -362,9 +401,9 @@ export fn copysign(x: types::floating, y: types::floating) f64 = {
assert(copysign(100.0f64, 0.0f64) == 100.0f64);
assert(copysign(100.0f64, -0.0f64) == -100.0f64);
assert(copysign(0.0f64, 100.0f64) == 0.0f64);
-	assert(sign(copysign(0.0f64, 100.0f64)) > 0);
+	assert(signf(copysign(0.0f64, 100.0f64)) > 0);
assert(copysign(0.0f64, -100.0f64) == 0.0f64);
-	assert(sign(copysign(0.0f64, -100.0f64)) < 0);
+	assert(signf(copysign(0.0f64, -100.0f64)) < 0);
};

// Takes a potentially subnormal f64 n and returns a normal f64 normal_float
@@ -389,7 +428,7 @@ export fn normalizef32(n: f32) (f64, i64) = {
return (n, 0);
};

-@test fn test_normalize() void = {
+@test fn normalize() void = {
let res = normalizef64(5.0e-320);
assert(res.0 > F64_MIN_NORMAL);
assert(res.1 < 0i64);
@@ -398,8 +437,8 @@ export fn normalizef32(n: f32) (f64, i64) = {
assert(res.1 == 0i64);
};

-// Breaks a f64 down into its mantissa and exponent. The mantissa will be between 0.5 and
-// 1.
+// Breaks a f64 down into its mantissa and exponent. The mantissa will be
+// between 0.5 and 1.
export fn frexpf64(n: f64) (f64, i64) = {
if (isnan(n) || isinf(n) || n == 0f64) {
return (n, 0);
@@ -411,12 +450,13 @@ export fn frexpf64(n: f64) (f64, i64) = {
const raw_exp: u64 = (bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK;
const exp: i64 = normalization_exp +
(raw_exp: i64) - (F64_EXPONENT_BIAS: i64) + 1;
-	const mantissa: f64 = f64frombits((bits & F64_EXP_REMOVAL_MASK) | F64_EXP_ZERO);
+	const mantissa: f64 =
+		f64frombits((bits & F64_EXP_REMOVAL_MASK) | F64_EXP_ZERO);
return (mantissa, exp);
};

-// Breaks a f32 down into its mantissa and exponent. The mantissa will be between 0.5 and
-// 1.
+// Breaks a f32 down into its mantissa and exponent. The mantissa will be
+// between 0.5 and 1.
export fn frexpf32(n: f32) (f64, i64) = {
if (isnan(n) || isinf(n) || n == 0f32) {
return (n, 0);
@@ -429,12 +469,13 @@ export fn frexpf32(n: f32) (f64, i64) = {
const exp: i64 = normalization_exp +
(raw_exp: i64) - (F32_EXPONENT_BIAS: i64) + 1;
-	const mantissa: f32 = f32frombits((bits & F32_EXP_REMOVAL_MASK) | F32_EXP_ZERO);
+	const mantissa: f32 =
+		f32frombits((bits & F32_EXP_REMOVAL_MASK) | F32_EXP_ZERO);
return (mantissa, exp);
};

-// Breaks a float down into its mantissa and exponent. The mantissa will be between 0.5
-// and 1.
+// Breaks a float down into its mantissa and exponent. The mantissa will be
+// between 0.5 and 1.
export fn frexp(n: types::floating) (f64, i64) = {
return match (n) {
n: f64 => frexpf64(n),
@@ -442,7 +483,7 @@ export fn frexp(n: types::floating) (f64, i64) = {
};
};

-@test fn test_frexp() void = {
+@test fn frexp() void = {
let res = frexp(3.0f64);
assert(res.0 == 0.75f64);
assert(res.1 == 2i64);
@@ -464,7 +505,8 @@ export fn ldexpf64(mantissa: f64, exp: i64) f64 = {
const normal_float = normalized.0;
const normalization_exp = normalized.1;
const bits = f64bits(normal_float);
-	const mantissa_exp = (((bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK): i64) -
+	const mantissa_exp =
+		(((bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK): i64) -
(F64_EXPONENT_BIAS: i64);
let res_exp = exp + normalization_exp + mantissa_exp;
// Underflow
@@ -503,7 +545,8 @@ export fn ldexpf32(mantissa: f32, exp: i64) f32 = {
const normal_float = normalized.0;
const normalization_exp = normalized.1;
const bits = f32bits(normal_float);
-	const mantissa_exp = (((bits >> F32_MANTISSA_BITS) & F32_EXPONENT_MASK): i32) -
+	const mantissa_exp =
+		(((bits >> F32_MANTISSA_BITS) & F32_EXPONENT_MASK): i32) -
(F32_EXPONENT_BIAS: i32);
let res_exp = exp + normalization_exp + mantissa_exp;
// Underflow
@@ -533,7 +576,7 @@ export fn ldexpf32(mantissa: f32, exp: i64) f32 = {
return subnormal_factor * f32frombits(res);
};

-@test fn test_frexp_ldexp() void = {
+@test fn frexp_ldexp() void = {
const tests64: [_]f64 = [INF, -INF,
0.0, 1.0, -1.0, 2.42, 123456789.0,
F64_MIN_NORMAL, F64_MAX_NORMAL,
@@ -545,7 +588,7 @@ export fn ldexpf32(mantissa: f32, exp: i64) f32 = {
};
assert(ldexpf64(1.0f64, -1076i64) == 0.0f64);
assert(ldexpf64(-1.0f64, -1076i64) == -0.0f64);
-	assert(sign(ldexpf64(-1.0f64, -1076i64)) < 0);
+	assert(signf(ldexpf64(-1.0f64, -1076i64)) < 0);
assert(ldexpf64(2.0f64, 1024i64) == INF);
assert(ldexpf64(-2.0f64, 1024i64) == -INF);

@@ -560,22 +603,22 @@ export fn ldexpf32(mantissa: f32, exp: i64) f32 = {
};
assert(ldexpf32(1.0f32, -1076i32) == 0.0f32);
assert(ldexpf32(-1.0f32, -1076i32) == -0.0f32);
-	assert(sign(ldexpf32(-1.0f32, -1076i32)) < 0);
+	assert(signf(ldexpf32(-1.0f32, -1076i32)) < 0);
assert(ldexpf32(2.0f32, 1024i32) == INF);
assert(ldexpf32(-2.0f32, 1024i32) == -INF);
};

// Returns the integer and fractional parts of an f64.
-export fn modf64(n: f64) (f64, f64) = {
+export fn modf64(n: f64) (i64, f64) = {
if (n < 1.0f64) {
if (n < 0.0f64) {
let positive_parts = modf64(-n);
return (-positive_parts.0, -positive_parts.1);
};
if (n == 0.0f64) {
-			return (n, n);
+			return ((n: i64), n);
};
-		return (0.0f64, n);
+		return (0i64, n);
};
let bits = f64bits(n);
const exp = (((bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK): i64) -
@@ -589,20 +632,20 @@ export fn modf64(n: f64) (f64, f64) = {
};
const int_part = f64frombits(bits);
const frac_part = n - int_part;
-	return (int_part, frac_part);
+	return ((int_part: i64), frac_part);
};

// Returns the integer and fractional parts of an f32.
-export fn modf32(n: f32) (f32, f32) = {
+export fn modf32(n: f32) (i32, f32) = {
if (n < 1.0f32) {
if (n < 0.0f32) {
let positive_parts = modf32(-n);
return (-positive_parts.0, -positive_parts.1);
};
if (n == 0.0f32) {
-			return (n, n);
+			return ((n: i32), n);
};
-		return (0.0f32, n);
+		return (0i32, n);
};
let bits = f32bits(n);
const exp = (((bits >> F32_MANTISSA_BITS) & F32_EXPONENT_MASK): i32) -
@@ -616,45 +659,51 @@ export fn modf32(n: f32) (f32, f32) = {
};
const int_part = f32frombits(bits);
const frac_part = n - int_part;
-	return (int_part, frac_part);
+	return ((int_part: i32), frac_part);
};

-@test fn test_modf() void = {
+@test fn modf() void = {
// 64
let res = modf64(1.75f64);
-	assert(res.0 == 1.0f64);
+	assert(res.0 == 1i64);
assert(res.1 == 0.75f64);
res = modf64(0.75f64);
-	assert(res.0 == 0.0f64);
+	assert(res.0 == 0i64);
assert(res.1 == 0.75f64);
res = modf64(-0.75f64);
-	assert(res.0 == -0.0f64);
+	assert(res.0 == -0i64);
assert(res.1 == -0.75f64);
res = modf64(0.0f64);
-	assert(res.0 == 0.0f64);
+	assert(res.0 == 0i64);
assert(res.1 == 0.0f64);
-	assert(sign(res.1) > 0);
+	assert(signf(res.1) > 0);
res = modf64(-0.0f64);
-	assert(res.0 == -0.0f64);
+	assert(res.0 == -0i64);
assert(res.1 == -0.0f64);
-	assert(sign(res.1) < 0);
+	assert(signf(res.1) < 0);
+	res = modf64(23.50f64);
+	assert(res.0 == 23i64);
+	assert(res.1 == 0.50f64);

// 32
let res = modf32(1.75f32);
-	assert(res.0 == 1.0f32);
+	assert(res.0 == 1i32);
assert(res.1 == 0.75f32);
res = modf32(0.75f32);
-	assert(res.0 == 0.0f32);
+	assert(res.0 == 0i32);
assert(res.1 == 0.75f32);
res = modf32(-0.75f32);
-	assert(res.0 == -0.0f32);
+	assert(res.0 == -0i32);
assert(res.1 == -0.75f32);
res = modf32(0.0f32);
-	assert(res.0 == 0.0f32);
+	assert(res.0 == 0i32);
assert(res.1 == 0.0f32);
-	assert(sign(res.1) > 0);
+	assert(signf(res.1) > 0);
res = modf32(-0.0f32);
-	assert(res.0 == -0.0f32);
-	assert(res.0 == -0.0f32);
-	assert(sign(res.1) < 0);
+	assert(res.0 == -0i32);
+	assert(res.1 == -0.0f32);
+	assert(signf(res.1) < 0);
+	res = modf32(23.50f32);
+	assert(res.0 == 23i32);
+	assert(res.1 == 0.50f32);
};
diff --git a/math/ints.ha b/math/ints.ha
@@ -0,0 +1,138 @@
+use types;
+
+// Returns the absolute value of signed integer n.
+export fn absi8(n: i8) i8 = {
+	if (n < 0i8) {
+		return -n;
+	} else {
+		return n;
+	};
+};
+
+// Returns the absolute value of signed integer n.
+export fn absi16(n: i16) i16 = {
+	if (n < 0i16) {
+		return -n;
+	} else {
+		return n;
+	};
+};
+
+// Returns the absolute value of signed integer n.
+export fn absi32(n: i32) i32 = {
+	if (n < 0i32) {
+		return -n;
+	} else {
+		return n;
+	};
+};
+
+// Returns the absolute value of signed integer n.
+export fn absi64(n: i64) i64 = {
+	if (n < 0i64) {
+		return -n;
+	} else {
+		return n;
+	};
+};
+
+// Returns the absolute value of signed integer n.
+export fn absi(n: types::integer) i64 = {
+	return match (n) {
+		n: i8 => (absi8(n): i64),
+		n: i16 => (absi16(n): i64),
+		n: i32 => (absi32(n): i64),
+		n: i64 => (absi64(n): i64),
+	};
+};
+
+@test fn absi() void = {
+	assert(absi8(2i8) == 2i8);
+	assert(absi8(-2i8) == 2i8);
+	assert(absi16(2i16) == 2i16);
+	assert(absi16(-2i16) == 2i16);
+	assert(absi32(2i32) == 2i32);
+	assert(absi32(-2i32) == 2i32);
+	assert(absi64(2i64) == 2i64);
+	assert(absi64(-2i64) == 2i64);
+	assert(absi(2i8) == 2i64);
+	assert(absi(-2i8) == 2i64);
+	assert(absi(2i16) == 2i64);
+	assert(absi(-2i16) == 2i64);
+	assert(absi(2i32) == 2i64);
+	assert(absi(-2i32) == 2i64);
+	assert(absi(2i64) == 2i64);
+	assert(absi(-2i64) == 2i64);
+};
+
+// Return 1 if n is positive, -1 if it's negative and 0 if it's 0.
+export fn signi8(n: i8) i8 = {
+	if (n > 0i8) {
+		return 1i8;
+	};
+	if (n < 0i8) {
+		return -1i8;
+	};
+	return 0i8;
+};
+
+// Return 1 if n is positive, -1 if it's negative and 0 if it's 0.
+export fn signi16(n: i16) i16 = {
+	if (n > 0i16) {
+		return 1i16;
+	};
+	if (n < 0i16) {
+		return -1i16;
+	};
+	return 0i16;
+};
+
+// Return 1 if n is positive, -1 if it's negative and 0 if it's 0.
+export fn signi32(n: i32) i32 = {
+	if (n > 0i32) {
+		return 1i32;
+	};
+	if (n < 0i32) {
+		return -1i32;
+	};
+	return 0i32;
+};
+
+// Return 1 if n is positive, -1 if it's negative and 0 if it's 0.
+export fn signi64(n: i64) i64 = {
+	if (n > 0i64) {
+		return 1i64;
+	};
+	if (n < 0i64) {
+		return -1i64;
+	};
+	return 0i64;
+};
+
+// Return 1 if n is positive, -1 if it's negative and 0 if it's 0.
+export fn signi(n: types::integer) i64 = {
+	return match (n) {
+		n: i8 => (signi8(n): i64),
+		n: i16 => (signi16(n): i64),
+		n: i32 => (signi32(n): i64),
+		n: i64 => (signi64(n): i64),
+	};
+};
+
+@test fn signi() void = {
+	assert(signi8(2i8) == 1i8);
+	assert(signi8(-2i8) == -1i8);
+	assert(signi8(0i8) == 0i8);
+	assert(signi16(2i16) == 1i16);
+	assert(signi16(-2i16) == -1i16);
+	assert(signi16(0i16) == 0i16);
+	assert(signi32(2i32) == 1i32);
+	assert(signi32(-2i32) == -1i32);
+	assert(signi32(0i32) == 0i32);
+	assert(signi64(2i64) == 1i64);
+	assert(signi64(-2i64) == -1i64);
+	assert(signi64(0i64) == 0i64);
+	assert(signi(2i16) == 1i64);
+	assert(signi(-2i16) == -1i64);
+	assert(signi(0i16) == 0i64);
+};
diff --git a/math/math.ha b/math/math.ha
@@ -0,0 +1,690 @@
+// Sections of the code below, in particular log() and exp(), are based on Go's
+// implementation, which is, in turn, based on FreeBSD's. The original C code,
+// as well as the respective comments and constants are from
+// /usr/src/lib/msun/src/e_log.c.
+//
+// ====================================================
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// ====================================================
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+//    * Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//    * Redistributions in binary form must reproduce the above
+// copyright notice, this list of conditions and the following disclaimer
+// in the documentation and/or other materials provided with the
+// distribution.
+//    * Neither the name of Google Inc. nor the names of its
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+// ====================================================
+
+use types;
+
+// Returns whether x and y are within tol of each other.
+export fn eqwithinf64(x: f64, y: f64, tol: f64) bool = {
+	return absf64(x - y) < tol;
+};
+
+// Returns whether x and y are within tol of each other.
+export fn eqwithinf32(x: f32, y: f32, tol: f32) bool = {
+	return absf32(x - y) < tol;
+};
+
+// Returns whether x and y are within tol of each other.
+export fn eqwithin(x: types::floating, y: types::floating,
+		tol: types::floating) bool = {
+	return match (x) {
+		n: f64 => eqwithinf64(n, y as f64, tol as f64),
+		n: f32 => eqwithinf32(n, y as f32, tol as f32),
+	};
+};
+
+@test fn equalwithin() void = {
+	assert(eqwithin(1.0f64, 2.0f64, 2.0f64));
+	assert(eqwithin(1.0f32, 2.0f32, 2.0f32));
+	assert(!eqwithin(1.0005f32, 1.0004f32, 0.00001f32));
+};
+
+// e - https://oeis.org/A001113
+export def E: f64 = 2.71828182845904523536028747135266249775724709369995957496696763;
+// pi - https://oeis.org/A000796
+export def PI: f64 = 3.14159265358979323846264338327950288419716939937510582097494459;
+// phi - https://oeis.org/A001622
+export def PHI: f64 = 1.61803398874989484820458683436563811772030917980576286213544862;
+// sqrt(2) - https://oeis.org/A002193
+def SQRT_2: f64 = 1.41421356237309504880168872420969807856967187537694807317667974;
+// sqrt(e) - https://oeis.org/A019774
+def SQRT_E: f64 = 1.64872127070012814684865078781416357165377610071014801157507931;
+// sqrt(pi) - https://oeis.org/A002161
+def SQRT_PI: f64 = 1.77245385090551602729816748334114518279754945612238712821380779;
+// sqrt(phi) - https://oeis.org/A139339
+def SQRT_PHI: f64 = 1.27201964951406896425242246173749149171560804184009624861664038;
+// ln(2) - https://oeis.org/A002162
+def LN_2: f64 = 0.693147180559945309417232121458176568075500134360255254120680009;
+def LN2_HI: f64 = 6.93147180369123816490e-01;
+def LN2_LO: f64 = 1.90821492927058770002e-10;
+// log_{2}(e)
+def LOG2_E: f64 = 1.0f64 / LN_2;
+// ln(10) - https://oeis.org/A002392
+def LN_10: f64 = 2.30258509299404568401799145468436420760110148862877297603332790;
+// log_{10}(e)
+def LOG10_E: f64 = 1.0f64 / LN_10;
+
+// __ieee754_log(x)
+// Return the logarithm of x
+//
+// Method :
+//   1. Argument Reduction: find k and f such that
+//			x = 2**k * (1+f),
+//	   where  sqrt(2)/2 < 1+f < sqrt(2) .
+//
+//   2. Approximation of log(1+f).
+//	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+//		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+//	     	 = 2s + s*R
+//      We use a special Reme algorithm on [0,0.1716] to generate
+//	a polynomial of degree 14 to approximate R.  The maximum error
+//	of this polynomial approximation is bounded by 2**-58.45. In
+//	other words,
+//		        2      4      6      8      10      12      14
+//	    R(z) ~ L1*s +L2*s +L3*s +L4*s +L5*s  +L6*s  +L7*s
+//	(the values of L1 to L7 are listed in the program) and
+//	    |      2          14          |     -58.45
+//	    | L1*s +...+L7*s    -  R(z) | <= 2
+//	    |                             |
+//	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+//	In order to guarantee error in log below 1ulp, we compute log by
+//		log(1+f) = f - s*(f - R)		(if f is not too large)
+//		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
+//
+//	3. Finally,  log(x) = k*Ln2 + log(1+f).
+//			    = k*Ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*Ln2_lo)))
+//	   Here Ln2 is split into two floating point number:
+//			Ln2_hi + Ln2_lo,
+//	   where n*Ln2_hi is always exact for |n| < 2000.
+//
+// Special cases:
+//	log(x) is NaN with signal if x < 0 (including -INF) ;
+//	log(+INF) is +INF; log(0) is -INF with signal;
+//	log(NaN) is that NaN with no signal.
+//
+// Accuracy:
+//	according to an error analysis, the error is always less than
+//	1 ulp (unit in the last place).
+//
+// Constants:
+// The hexadecimal values are the intended ones for the following
+// constants. The decimal values may be used, provided that the
+// compiler will convert from decimal to binary accurately enough
+// to produce the hexadecimal values shown.
+
+// Returns the natural logarithm of x.
+export fn logf64(x: f64) f64 = {
+	const L1 = 6.666666666666735130e-01; // 3FE55555 55555593
+	const L2 = 3.999999999940941908e-01; // 3FD99999 9997FA04
+	const L3 = 2.857142874366239149e-01; // 3FD24924 94229359
+	const L4 = 2.222219843214978396e-01; // 3FCC71C5 1D8E78AF
+	const L5 = 1.818357216161805012e-01; // 3FC74664 96CB03DE
+	const L6 = 1.531383769920937332e-01; // 3FC39A09 D078C69F
+	const L7 = 1.479819860511658591e-01; // 3FC2F112 DF3E5244
+
+	// special cases
+	if (isnan(x) || isposinf(x)) {
+		return x;
+	} else if (x < 0.0f64) {
+		return NAN;
+	} else if (x == 0.0f64) {
+		return -INF;
+	};
+
+	// Reduce
+	const f1_and_ki = frexp(x);
+	let f1 = f1_and_ki.0;
+	let ki = f1_and_ki.1;
+	if (f1 < (SQRT_2 / 2.0f64)) {
+		f1 *= 2.0f64;
+		ki -= 1i64;
+	};
+	let f = f1 - 1.0f64;
+	let k = (ki: f64);
+
+	// Compute
+	const s = f / (2.0f64 + f);
+	const s2 = s * s;
+	const s4 = s2 * s2;
+	const t1 = s2 * (L1 + s4 * (L3 + s4 * (L5 + s4 * L7)));
+	const t2 = s4 * (L2 + s4 * (L4 + s4 * L6));
+	const R = t1 + t2;
+	const hfsq = 0.5f64 * f * f;
+	return k * LN2_HI - ((hfsq - (s * (hfsq + R) + k * LN2_LO)) - f);
+};
+
+@test fn logf64() void = {
+	assert(logf64(E) == 1.0f64);
+	assert(logf64(54.598150033144239078110261202860878402790f64) == 4.0f64);
+	assert(isnan(logf64(-1.0f64)));
+	assert(isposinf(logf64(INF)));
+	assert(isneginf(logf64(0.0f64)));
+	assert(isnan(logf64(NAN)));
+};
+
+// exp(x)
+// Returns the exponential of x.
+//
+// Method
+//   1. Argument reduction:
+//      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+//      Given x, find r and integer k such that
+//
+//               x = k*ln2 + r,  |r| <= 0.5*ln2.
+//
+//      Here r will be represented as r = hi-lo for better
+//      accuracy.
+//
+//   2. Approximation of exp(r) by a special rational function on
+//      the interval [0,0.34658]:
+//      Write
+//          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+//      We use a special Remez algorithm on [0,0.34658] to generate
+//      a polynomial of degree 5 to approximate R. The maximum error
+//      of this polynomial approximation is bounded by 2**-59. In
+//      other words,
+//          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+//      (where z=r*r, and the values of P1 to P5 are listed below)
+//      and
+//          |                  5          |     -59
+//          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
+//          |                             |
+//      The computation of exp(r) thus becomes
+//                             2*r
+//              exp(r) = 1 + -------
+//                            R - r
+//                                 r*R1(r)
+//                     = 1 + r + ----------- (for better accuracy)
+//                                2 - R1(r)
+//      where
+//                               2       4             10
+//              R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
+//
+//   3. Scale back to obtain exp(x):
+//      From step 1, we have
+//         exp(x) = 2**k * exp(r)
+//
+// Special cases:
+//      exp(INF) is INF, exp(NaN) is NaN;
+//      exp(-INF) is 0, and
+//      for finite argument, only exp(0)=1 is exact.
+//
+// Accuracy:
+//      according to an error analysis, the error is always less than
+//      1 ulp (unit in the last place).
+//
+// Misc. info.
+//      For IEEE double
+//          if x >  7.09782712893383973096e+02 then exp(x) overflow
+//          if x < -7.45133219101941108420e+02 then exp(x) underflow
+//
+// Constants:
+// The hexadecimal values are the intended ones for the following
+// constants. The decimal values may be used, provided that the
+// compiler will convert from decimal to binary accurately enough
+// to produce the hexadecimal values shown.
+
+// Returns e^r * 2^k where r = hi - lo and |r| <= (ln(2) / 2).
+export fn expmultif64(hi: f64, lo: f64, k: i64) f64 = {
+	const P1 = 1.66666666666666657415e-01; // 0x3FC55555; 0x55555555
+	const P2 = -2.77777777770155933842e-03; // 0xBF66C16C; 0x16BEBD93
+	const P3 = 6.61375632143793436117e-05; // 0x3F11566A; 0xAF25DE2C
+	const P4 = -1.65339022054652515390e-06; // 0xBEBBBD41; 0xC5D26BF1
+	const P5 = 4.13813679705723846039e-08; // 0x3E663769; 0x72BEA4D0
+
+	let r = hi - lo;
+	let t = r * r;
+	let c = r - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
+	let y = 1.0f64 - ((lo - (r * c) / (2.0f64 - c)) - hi);
+	return ldexpf64(y, k);
+};
+
+// Returns e^x.
+export fn expf64(x: f64) f64 = {
+	const overflow = 7.09782712893383973096e+02;
+	const underflow = -7.45133219101941108420e+02;
+	const near_zero = 1.0f64 / ((1i64 << 28i64): f64);
+
+	// Special cases
+	if (isnan(x) || isposinf(x)) {
+		return x;
+	} else if (isneginf(x)) {
+		return 0.0f64;
+	} else if (x > overflow) {
+		return INF;
+	} else if (x < underflow) {
+		return 0.0f64;
+	} else if (-near_zero < x && x < near_zero) {
+		return 1.0f64 + x;
+	};
+
+	// Reduce; computed as r = hi - lo for extra precision.
+	let k = 0i64;
+	if (x < 0.0f64) {
+		k = (((LOG2_E * x) - 0.5): i64);
+	} else if (x > 0.0f64) {
+		k = (((LOG2_E * x) + 0.5): i64);
+	};
+	let hi = x - ((k: f64) * LN2_HI);
+	let lo = (k: f64) * LN2_LO;
+
+	// Compute
+	return expmultif64(hi, lo, k);
+};
+
+// Returns 2^x.
+export fn exp2f64(x: f64) f64 = {
+	const overflow = 1.0239999999999999e+03;
+	const underflow = -1.0740e+03;
+
+	// Special cases
+	if (isnan(x) || isposinf(x)) {
+		return x;
+	} else if (isneginf(x)) {
+		return 0.0f64;
+	} else if (x > overflow) {
+		return INF;
+	} else if (x < underflow) {
+		return 0.0f64;
+	};
+
+	// Argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2.
+	// Computed as r = hi - lo for extra precision.
+	let k = 0i64;
+	if (x > 0.0f64) {
+		k = ((x + 0.5): i64);
+	} else if (x < 0.0f64) {
+		k = ((x - 0.5): i64);
+	};
+	let t = x - (k: f64);
+	let hi = t * LN2_HI;
+	let lo = -t * LN2_LO;
+
+	// Compute
+	return expmultif64(hi, lo, k);
+};
+
+@test fn expf64() void = {
+	assert(expf64(1.0f64) == E);
+	assert(isnan(expf64(NAN)));
+	assert(isinf(expf64(INF)));
+	assert(expf64(-INF) == 0.0f64);
+	assert(isinf(expf64(99999.0f64)));
+	assert(expf64(-99999.0f64) == 0.0f64);
+	assert(expf64(0.5e-20) == 1.0f64);
+};
+
+@test fn exp2f64() void = {
+	assert(exp2f64(0.0f64) == 1.0f64);
+	assert(exp2f64(3.0f64) == 8.0f64);
+	assert(exp2f64(-2.0f64) == 0.25f64);
+	assert(!isinf(exp2f64(256.0f64)));
+	assert(isinf(exp2f64(99999.0f64)));
+	assert(exp2f64(-99999.0f64) == 0.0f64);
+	assert(isnan(exp2f64(NAN)));
+	assert(isinf(exp2f64(INF)));
+	assert(exp2f64(-INF) == 0.0f64);
+};
+
+// __ieee754_sqrt(x)
+// Return correctly rounded sqrt.
+//           -----------------------------------------
+//           | Use the hardware sqrt if you have one |
+//           -----------------------------------------
+// Method:
+//   Bit by bit method using integer arithmetic. (Slow, but portable)
+//   1. Normalization
+//      Scale x to y in [1,4) with even powers of 2:
+//      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then
+//              sqrt(x) = 2**k * sqrt(y)
+//   2. Bit by bit computation
+//      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
+//           i                                                   0
+//                                     i+1         2
+//          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)
+//           i      i            i                 i
+//
+//      To compute q    from q , one checks whether
+//                  i+1       i
+//
+//                            -(i+1) 2
+//                      (q + 2      )  <= y.                     (2)
+//                        i
+//                                                            -(i+1)
+//      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
+//                             i+1   i             i+1   i
+//
+//      With some algebraic manipulation, it is not difficult to see
+//      that (2) is equivalent to
+//                             -(i+1)
+//                      s  +  2       <= y                       (3)
+//                       i                i
+//
+//      The advantage of (3) is that s  and y  can be computed by
+//                                    i      i
+//      the following recurrence formula:
+//          if (3) is false
+//
+//          s     =  s  ,       y    = y   ;                     (4)
+//           i+1      i          i+1    i
+//
+//      otherwise,
+//                         -i                      -(i+1)
+//          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)
+//           i+1      i          i+1    i     i
+//
+//      One may easily use induction to prove (4) and (5).
+//      Note. Since the left hand side of (3) contain only i+2 bits,
+//            it is not necessary to do a full (53-bit) comparison
+//            in (3).
+//   3. Final rounding
+//      After generating the 53 bits result, we compute one more bit.
+//      Together with the remainder, we can decide whether the
+//      result is exact, bigger than 1/2ulp, or less than 1/2ulp
+//      (it will never equal to 1/2ulp).
+//      The rounding mode can be detected by checking whether
+//      huge + tiny is equal to huge, and whether huge - tiny is
+//      equal to huge for some floating point number "huge" and "tiny".
+
+// Returns the square root of x.
+export fn sqrtf64(x: f64) f64 = {
+	// Special cases
+	if (x == 0.0f64) {
+		return x;
+	} else if (isnan(x) || isposinf(x)) {
+		return x;
+	} else if (x < 0.0f64) {
+		return NAN;
+	};
+
+	let bits = f64bits(x);
+
+	// Normalize x
+	let exp = (((bits >> F64_MANTISSA_BITS) & F64_EXPONENT_MASK): i64);
+	if (exp == 0i64) {
+		// Subnormal x
+		for (bits & (1 << F64_MANTISSA_BITS) == 0) {
+			bits <<= 1;
+			exp -= 1;
+		};
+		exp += 1;
+	};
+	// Unbias exponent
+	exp -= (F64_EXPONENT_BIAS: i64);
+	bits = bits & ~(F64_EXPONENT_MASK << F64_MANTISSA_BITS);
+	bits = bits | (1u64 << (F64_MANTISSA_BITS: u64));
+	// Odd exp, double x to make it even
+	if (exp & 1i64 == 1i64) {
+		bits <<= 1;
+	};
+	// exp = exp/2, exponent of square root
+	exp >>= 1;
+	// Generate sqrt(x) bit by bit
+	bits <<= 1;
+	// q = sqrt(x)
+	let q = 0u64;
+	let s = 0u64;
+	// r = moving bit from MSB to LSB
+	let r = ((1u64 << (F64_MANTISSA_BITS + 1u64)): u64);
+	for (r != 0) {
+		let t = s + r;
+		if (t <= bits) {
+			s = t + r;
+			bits -= t;
+			q += r;
+		};
+		bits <<= 1u64;
+		r >>= 1u64;
+	};
+	// Final rounding
+	if (bits != 0) {
+		// Remainder, result not exact
+		// Round according to extra bit
+		q += q & 1;
+	};
+	// significand + biased exponent
+	bits = (q >> 1) + (
+		((exp - 1i64 + (F64_EXPONENT_BIAS: i64)): u64) <<
+		F64_MANTISSA_BITS);
+	return f64frombits(bits);
+};
+
+@test fn sqrt() void = {
+	assert(sqrtf64(2.0f64) == SQRT_2);
+	assert(sqrtf64(4.0f64) == 2.0f64);
+	assert(sqrtf64(16.0f64) == 4.0f64);
+	assert(sqrtf64(65536.0f64) == 256.0f64);
+	assert(sqrtf64(powf64(123.0f64, 2.0f64)) == 123.0f64);
+	assert(sqrtf64(0.0f64) == 0.0f64);
+	assert(isnan(sqrtf64(NAN)));
+	assert(isposinf(sqrtf64(INF)));
+	assert(isnan(sqrtf64(-2.0f64)));
+};
+
+fn is_f64_odd_int(x: f64) bool = {
+	let x_int_frac = modf64(x);
+	let x_int = x_int_frac.0;
+	let x_frac = x_int_frac.1;
+	let has_no_frac = (x_frac == 0.0f64);
+	let is_odd = ((x_int & 1i64) == 1i64);
+	return has_no_frac && is_odd;
+};
+
+// Returns x^p.
+export fn powf64(x: f64, p: f64) f64 = {
+	if (x == 1.0f64 || p == 0.0f64) {
+		return 1.0f64;
+	} else if (p == 1.0f64) {
+		return x;
+	} else if (isnan(x)) {
+		return NAN;
+	} else if (isnan(p)) {
+		return NAN;
+	} else if (x == 0.0f64) {
+		if (p < 0.0f64) {
+			if (is_f64_odd_int(p)) {
+				return copysignf64(INF, x);
+			} else {
+				return INF;
+			};
+		} else if (p > 0.0f64) {
+			if (is_f64_odd_int(p)) {
+				return x;
+			} else {
+				return 0.0f64;
+			};
+		};
+	} else if (isinf(p)) {
+		if (x == -1.0f64) {
+			return 1.0f64;
+		} else if ((absf64(x) < 1.0f64) == isposinf(p)) {
+			return 0.0f64;
+		};
+		return INF;
+	} else if (isinf(x)) {
+		if (isneginf(x)) {
+			return powf64(-0.0f64, -p);
+		} else if (p < 0.0f64) {
+			return 0.0f64;
+		} else if (p > 0.0f64) {
+			return INF;
+		};
+	} else if (p == 0.5f64) {
+		return sqrtf64(x);
+	} else if (p == -0.5f64) {
+		return 1.0f64 / sqrtf64(x);
+	};
+
+	let x_parts = frexp(x);
+	let x_mantissa = x_parts.0;
+	let x_exp = x_parts.1;
+
+	let p_int_frac = modf64(absf64(p));
+	let p_int = p_int_frac.0;
+	let p_frac = p_int_frac.1;
+
+	let res_mantissa = 1.0f64;
+	let res_exp = 0i64;
+
+	// The method used later in this function doesn't apply to fractional
+	// powers, so we have to handle these separately with
+	// x^p = e^{p * ln(x)}
+	if (p_frac != 0.0f64) {
+		if (p_frac > 0.5f64) {
+			p_frac -= 1.0f64;
+			p_int += 1i64;
+		};
+		res_mantissa = expf64(p_frac * logf64(x));
+	};
+
+	// Repeatedly square our number x, for each bit in our power p.
+	// If the current bit is 1 in p, add the respective power of x to our
+	// result.
+	for (let i = p_int; i != 0; i >>= 1) {
+		// Check for over/underflow.
+		if (x_exp <= -1i64 << (F64_EXPONENT_BITS: i64)) {
+			return 0.0f64;
+	};
+		if (x_exp >= 1i64 << (F64_EXPONENT_BITS: i64)) {
+			return INF;
+		};
+		// Perform squaring.
+		if (i & 1i64 == 1i64) {
+			res_mantissa *= x_mantissa;
+			res_exp += x_exp;
+		};
+		x_mantissa *= x_mantissa;
+		x_exp <<= 1;
+		// Correct mantisa to be in [0.5, 1).
+		if (x_mantissa < 0.5f64) {
+			x_mantissa += x_mantissa;
+			x_exp -= 1;
+		};
+	};
+
+	if (p < 0.0f64) {
+		res_mantissa = 1.0f64 / res_mantissa;
+		res_exp = -res_exp;
+	};
+
+	let res = ldexpf64(res_mantissa, res_exp);
+	return res;
+};
+
+@test fn powf64() void = {
+	// Positive integer
+	assert(powf64(2.0f64, 2.0f64) == 4.0f64);
+	assert(powf64(3.0f64, 3.0f64) == 27.0f64);
+	// Very large positive integer
+	assert(!isinf(powf64(2.0f64, 1020.0f64)));
+	assert(isinf(powf64(2.0f64, 1050.0f64)));
+	// Negative integer
+	assert(powf64(2.0f64, -1.0f64) == 0.5f64);
+	assert(powf64(2.0f64, -2.0f64) == 0.25f64);
+	// Very small negative integer
+	assert(powf64(2.0f64, -1020.0f64) > 0.0f64);
+	assert(powf64(2.0f64, -1080.0f64) == 0.0f64);
+	// Positive fractional powers
+	assert(eqwithin(powf64(2.0f64, 1.5f64),
+		2.8284271247461900976033774f64,
+		1e-10f64));
+	assert(eqwithin(powf64(2.0f64, 5.5f64),
+		45.254833995939041561654039f64,
+		1e-10f64));
+	// Negative fractional powers
+	assert(eqwithin(powf64(2.0f64, -1.5f64),
+		0.3535533905932737622004221f64,
+		1e-10f64));
+	assert(eqwithin(powf64(2.0f64, -5.5f64),
+		0.0220970869120796101375263f64,
+		1e-10f64));
+
+	// Special cases
+	// pow(x, ±0) = 1 for any x
+	assert(powf64(123.0f64, 0.0f64) == 1.0f64);
+	// pow(1, y) = 1 for any y
+	assert(powf64(1.0f64, 123.0f64) == 1.0f64);
+	// pow(x, 1) = x for any x
+	assert(powf64(123.0f64, 1.0f64) == 123.0f64);
+	// pow(NaN, y) = NaN
+	assert(isnan(powf64(NAN, 123.0f64)));
+	// pow(x, NaN) = NaN
+	assert(isnan(powf64(123.0f64, NAN)));
+	// pow(±0, y) = ±Inf for y an odd integer < 0
+	assert(isposinf(powf64(0.0f64, -3.0f64)));
+	assert(isneginf(powf64(-0.0f64, -3.0f64)));
+	// pow(±0, -Inf) = +Inf
+	assert(isposinf(powf64(0.0f64, -INF)));
+	assert(isposinf(powf64(-0.0f64, -INF)));
+	// pow(±0, +Inf) = +0
+	assert(powf64(0.0f64, INF) == 0.0f64);
+	assert(powf64(-0.0f64, INF) == 0.0f64);
+	// pow(±0, y) = +Inf for finite y < 0 and not an odd integer
+	assert(isposinf(powf64(0.0f64, -2.0f64)));
+	assert(isposinf(powf64(-0.0f64, -2.0f64)));
+	//pow(±0, y) = ±0 for y an odd integer > 0
+	assert(powf64(0.0f64, 123.0f64) == 0.0f64);
+	let neg_zero = powf64(-0.0f64, 123.0f64);
+	assert(neg_zero == 0.0f64);
+	assert(isnegative(neg_zero));
+	// pow(±0, y) = +0 for finite y > 0 and not an odd integer
+	assert(powf64(0.0f64, 8.0f64) == 0.0f64);
+	// pow(-1, ±Inf) = 1
+	assert(powf64(-1.0f64, INF) == 1.0f64);
+	assert(powf64(-1.0f64, -INF) == 1.0f64);
+	// pow(x, +Inf) = +Inf for |x| > 1
+	assert(isposinf(powf64(123.0f64, INF)));
+	// pow(x, -Inf) = +0 for |x| > 1
+	assert(powf64(123.0f64, -INF) == 0.0f64);
+	// pow(x, +Inf) = +0 for |x| < 1
+	assert(powf64(0.5f64, INF) == 0.0f64);
+	assert(powf64(-0.5f64, INF) == 0.0f64);
+	// pow(x, -Inf) = +Inf for |x| < 1
+	assert(isposinf(powf64(0.5f64, -INF)));
+	assert(isposinf(powf64(-0.5f64, -INF)));
+	// pow(+Inf, y) = +Inf for y > 0
+	assert(isposinf(powf64(INF, 123.0f64)));
+	// pow(+Inf, y) = +0 for y < 0
+	assert(powf64(INF, -1.0f64) == 0.0f64);
+	// pow(-Inf, y) = pow(-0, -y)
+	assert(powf64(-INF, 123.0f64) == powf64(-0.0f64, -123.0f64));
+	// pow(x, y) = NaN for finite x < 0 and finite non-integer y
+	assert(isnan(powf64(-2.0f64, 1.23f64)));
+	// sqrt
+	assert(powf64(4.0f64, 0.5f64) == sqrtf64(4.0f64));
+	assert(powf64(4.0f64, 0.5f64) == 2.0f64);
+	assert(powf64(4.0f64, -0.5f64) == (1.0f64 / sqrtf64(4.0f64)));
+	assert(powf64(4.0f64, -0.5f64) == (1.0f64 / 2.0f64));
+};
diff --git a/scripts/gen-stdlib b/scripts/gen-stdlib
@@ -550,7 +550,7 @@ linux_vdso() {
math() {
printf '# math\n'
gen_srcs math \
-		floats.ha
+		math.ha floats.ha ints.ha
gen_ssa math types
}

diff --git a/stdlib.mk b/stdlib.mk
@@ -844,9 +844,11 @@ \$(HARECACHE)/linux/vdso/linux_vdso.ssa: \$(stdlib_linux_vdso_srcs) \$(stdlib_rt) \$
# math
# math
stdlib_math_srcs= \
-	\$(STDLIB)/math/floats.ha
+	\$(STDLIB)/math/math.ha \
+	\$(STDLIB)/math/floats.ha \
+	\$(STDLIB)/math/ints.ha

-\$(HARECACHE)/math/math.ssa: \$(stdlib_math_srcs) \$(stdlib_rt)
+\$(HARECACHE)/math/math.ssa: \$(stdlib_math_srcs) \$(stdlib_rt) \$(stdlib_types)
@printf 'HAREC \t\$@\n'
@mkdir -p \$(HARECACHE)/math
@HARECACHE=\$(HARECACHE) \$(HAREC) \$(HAREFLAGS) -o \$@ -Nmath \
@@ -2061,9 +2063,11 @@ \$(TESTCACHE)/linux/vdso/linux_vdso.ssa: \$(testlib_linux_vdso_srcs) \$(testlib_rt)
# math
# math
testlib_math_srcs= \
-	\$(STDLIB)/math/floats.ha
+	\$(STDLIB)/math/math.ha \
+	\$(STDLIB)/math/floats.ha \
+	\$(STDLIB)/math/ints.ha

-\$(TESTCACHE)/math/math.ssa: \$(testlib_math_srcs) \$(testlib_rt)
+\$(TESTCACHE)/math/math.ssa: \$(testlib_math_srcs) \$(testlib_rt) \$(testlib_types)
@printf 'HAREC \t\$@\n'
@mkdir -p \$(TESTCACHE)/math
@HARECACHE=\$(TESTCACHE) \$(HAREC) \$(TESTHAREFLAGS) -o \$@ -Nmath \
```