From 1cfc4198f520117d0528c66621d952d19468b4bc Mon Sep 17 00:00:00 2001
From: Ulises Jeremias Cornejo Fandos <ulisesjcf_24@outlook.com>
Date: Sun, 22 Aug 2021 18:35:28 -0300
Subject: [PATCH] vlib/math: Add a pure V backend for vlib/math (#11267)

---
 vlib/math/ROADMAP.md                   |   4 +
 vlib/math/abs.c.v                      |   8 +
 vlib/math/abs.js.v                     |   9 +
 vlib/math/abs.v                        |  18 +
 vlib/math/bits.c.v                     |  12 -
 vlib/math/bits.v                       |  49 +-
 vlib/math/cbrt.c.v                     |   9 +
 vlib/math/cbrt.js.v                    |   9 +
 vlib/math/cbrt.v                       |  52 ++
 vlib/math/complex/complex.v            |   2 +-
 vlib/math/const.v                      |   2 +
 vlib/math/div.c.v                      |   9 +
 vlib/math/div.v                        |  87 +++
 vlib/math/erf.c.v                      |  17 +
 vlib/math/erf.v                        | 938 ++++++++-----------------
 vlib/math/erf_test.v                   |  60 +-
 vlib/math/evaluate.v                   |  18 -
 vlib/math/exp.c.v                      |  17 +
 vlib/math/exp.js.v                     |  12 +
 vlib/math/exp.v                        | 214 ++++++
 vlib/math/{factorial => }/factorial.v  |  43 +-
 vlib/math/factorial/factorial_tables.v | 375 ----------
 vlib/math/factorial/factorial_test.v   |  14 -
 vlib/math/factorial_tables.v           | 711 +++++++++++++++++++
 vlib/math/factorial_test.v             |  13 +
 vlib/math/floor.c.v                    |  34 +
 vlib/math/floor.js.v                   |  34 +
 vlib/math/floor.v                      | 105 +++
 vlib/math/gamma.c.v                    |  17 +
 vlib/math/gamma.v                      | 335 +++++++++
 vlib/math/gamma_tables.v               | 163 +++++
 vlib/math/hypot.c.v                    |   9 +
 vlib/math/hypot.v                      |  24 +
 vlib/math/internal/machine.v           |  58 ++
 vlib/math/invhyp.v                     |  51 ++
 vlib/math/invtrig.c.v                  |  33 +
 vlib/math/invtrig.js.v                 |  33 +
 vlib/math/invtrig.v                    | 219 ++++++
 vlib/math/log.c.v                      |  25 +
 vlib/math/log.js.v                     |   9 +
 vlib/math/log.v                        |  76 ++
 vlib/math/math.c.v                     | 284 --------
 vlib/math/math.js.v                    | 281 --------
 vlib/math/math.v                       |  76 +-
 vlib/math/math_test.v                  | 247 ++++++-
 vlib/math/modf.v                       |  29 +
 vlib/math/nextafter.v                  |  45 ++
 vlib/math/poly.v                       |  65 ++
 vlib/math/pow.c.v                      |  17 +
 vlib/math/pow.js.v                     |   3 +
 vlib/math/pow.v                        |  37 +
 vlib/math/q_rsqrt.v                    |  12 +
 vlib/math/sin.c.v                      |  33 +
 vlib/math/sin.js.v                     |  17 +
 vlib/math/sin.v                        | 179 +++++
 vlib/math/sinh.c.v                     |  17 +
 vlib/math/sinh.js.v                    |  17 +
 vlib/math/sinh.v                       |  49 ++
 vlib/math/sqrt.c.v                     |  17 +
 vlib/math/sqrt.v                       |  37 +
 vlib/math/tan.c.v                      |  17 +
 vlib/math/tan.js.v                     |   9 +
 vlib/math/tan.v                        | 113 +++
 vlib/math/tanh.c.v                     |   9 +
 vlib/math/tanh.js.v                    |   9 +
 vlib/math/tanh.v                       |  45 ++
 vlib/math/{unsafe.c.v => unsafe.v}     |   0
 67 files changed, 3799 insertions(+), 1792 deletions(-)
 create mode 100644 vlib/math/ROADMAP.md
 create mode 100644 vlib/math/abs.c.v
 create mode 100644 vlib/math/abs.js.v
 create mode 100644 vlib/math/abs.v
 delete mode 100644 vlib/math/bits.c.v
 create mode 100644 vlib/math/cbrt.c.v
 create mode 100644 vlib/math/cbrt.js.v
 create mode 100644 vlib/math/cbrt.v
 create mode 100644 vlib/math/div.c.v
 create mode 100644 vlib/math/div.v
 create mode 100644 vlib/math/erf.c.v
 delete mode 100644 vlib/math/evaluate.v
 create mode 100644 vlib/math/exp.c.v
 create mode 100644 vlib/math/exp.js.v
 create mode 100644 vlib/math/exp.v
 rename vlib/math/{factorial => }/factorial.v (60%)
 delete mode 100644 vlib/math/factorial/factorial_tables.v
 delete mode 100644 vlib/math/factorial/factorial_test.v
 create mode 100644 vlib/math/factorial_tables.v
 create mode 100644 vlib/math/factorial_test.v
 create mode 100644 vlib/math/floor.c.v
 create mode 100644 vlib/math/floor.js.v
 create mode 100644 vlib/math/floor.v
 create mode 100644 vlib/math/gamma.c.v
 create mode 100644 vlib/math/gamma.v
 create mode 100644 vlib/math/gamma_tables.v
 create mode 100644 vlib/math/hypot.c.v
 create mode 100644 vlib/math/hypot.v
 create mode 100644 vlib/math/internal/machine.v
 create mode 100644 vlib/math/invhyp.v
 create mode 100644 vlib/math/invtrig.c.v
 create mode 100644 vlib/math/invtrig.js.v
 create mode 100644 vlib/math/invtrig.v
 create mode 100644 vlib/math/log.c.v
 create mode 100644 vlib/math/log.js.v
 create mode 100644 vlib/math/log.v
 delete mode 100644 vlib/math/math.js.v
 create mode 100644 vlib/math/modf.v
 create mode 100644 vlib/math/nextafter.v
 create mode 100644 vlib/math/poly.v
 create mode 100644 vlib/math/pow.c.v
 create mode 100644 vlib/math/pow.js.v
 create mode 100644 vlib/math/pow.v
 create mode 100644 vlib/math/q_rsqrt.v
 create mode 100644 vlib/math/sin.c.v
 create mode 100644 vlib/math/sin.js.v
 create mode 100644 vlib/math/sin.v
 create mode 100644 vlib/math/sinh.c.v
 create mode 100644 vlib/math/sinh.js.v
 create mode 100644 vlib/math/sinh.v
 create mode 100644 vlib/math/sqrt.c.v
 create mode 100644 vlib/math/sqrt.v
 create mode 100644 vlib/math/tan.c.v
 create mode 100644 vlib/math/tan.js.v
 create mode 100644 vlib/math/tan.v
 create mode 100644 vlib/math/tanh.c.v
 create mode 100644 vlib/math/tanh.js.v
 create mode 100644 vlib/math/tanh.v
 rename vlib/math/{unsafe.c.v => unsafe.v} (100%)

diff --git a/vlib/math/ROADMAP.md b/vlib/math/ROADMAP.md
new file mode 100644
index 000000000..856424417
--- /dev/null
+++ b/vlib/math/ROADMAP.md
@@ -0,0 +1,4 @@
+- [x] Move `vsl/vmath` to `vlib/math` as default backend
+- [ ] Implement `log` in pure V
+- [ ] Implement `pow` in pure V
+- [ ] Define functions for initial release of hardware implementations
diff --git a/vlib/math/abs.c.v b/vlib/math/abs.c.v
new file mode 100644
index 000000000..983ddcbb9
--- /dev/null
+++ b/vlib/math/abs.c.v
@@ -0,0 +1,8 @@
+module math
+
+fn C.fabs(x f64) f64
+
+[inline]
+pub fn abs(a f64) f64 {
+	return C.fabs(a)
+}
diff --git a/vlib/math/abs.js.v b/vlib/math/abs.js.v
new file mode 100644
index 000000000..abf90b085
--- /dev/null
+++ b/vlib/math/abs.js.v
@@ -0,0 +1,9 @@
+module math
+
+fn JS.Math.abs(x f64) f64
+
+// Returns the absolute value.
+[inline]
+pub fn abs(a f64) f64 {
+	return JS.Math.abs(a)
+}
diff --git a/vlib/math/abs.v b/vlib/math/abs.v
new file mode 100644
index 000000000..55bcff2b7
--- /dev/null
+++ b/vlib/math/abs.v
@@ -0,0 +1,18 @@
+module math
+
+// Returns the absolute value.
+[inline]
+pub fn abs(x f64) f64 {
+	if x > 0.0 {
+		return x
+	}
+	return -x
+}
+
+[inline]
+pub fn fabs(x f64) f64 {
+	if x > 0.0 {
+		return x
+	}
+	return -x
+}
diff --git a/vlib/math/bits.c.v b/vlib/math/bits.c.v
deleted file mode 100644
index 27aab5eeb..000000000
--- a/vlib/math/bits.c.v
+++ /dev/null
@@ -1,12 +0,0 @@
-module math
-
-// inf returns positive infinity if sign >= 0, negative infinity if sign < 0.
-pub fn inf(sign int) f64 {
-	v := if sign >= 0 { uvinf } else { uvneginf }
-	return f64_from_bits(v)
-}
-
-// nan returns an IEEE 754 ``not-a-number'' value.
-pub fn nan() f64 {
-	return f64_from_bits(uvnan)
-}
diff --git a/vlib/math/bits.v b/vlib/math/bits.v
index dbf2237f8..deaf96256 100644
--- a/vlib/math/bits.v
+++ b/vlib/math/bits.v
@@ -4,17 +4,29 @@
 module math
 
 const (
-	uvnan     = u64(0x7FF8000000000001)
-	uvinf     = u64(0x7FF0000000000000)
-	uvneginf  = u64(0xFFF0000000000000)
-	uvone     = u64(0x3FF0000000000000)
-	mask      = 0x7FF
-	shift     = 64 - 11 - 1
-	bias      = 1023
-	sign_mask = (u64(1) << 63)
-	frac_mask = ((u64(1) << u64(shift)) - u64(1))
+	uvnan                   = u64(0x7FF8000000000001)
+	uvinf                   = u64(0x7FF0000000000000)
+	uvneginf                = u64(0xFFF0000000000000)
+	uvone                   = u64(0x3FF0000000000000)
+	mask                    = 0x7FF
+	shift                   = 64 - 11 - 1
+	bias                    = 1023
+	normalize_smallest_mask = (u64(1) << 52)
+	sign_mask               = (u64(1) << 63)
+	frac_mask               = ((u64(1) << u64(shift)) - u64(1))
 )
 
+// inf returns positive infinity if sign >= 0, negative infinity if sign < 0.
+pub fn inf(sign int) f64 {
+	v := if sign >= 0 { math.uvinf } else { math.uvneginf }
+	return f64_from_bits(v)
+}
+
+// nan returns an IEEE 754 ``not-a-number'' value.
+pub fn nan() f64 {
+	return f64_from_bits(math.uvnan)
+}
+
 // is_nan reports whether f is an IEEE 754 ``not-a-number'' value.
 pub fn is_nan(f f64) bool {
 	// IEEE 754 says that only NaNs satisfy f != f.
@@ -36,13 +48,16 @@ pub fn is_inf(f f64, sign int) bool {
 	return (sign >= 0 && f > max_f64) || (sign <= 0 && f < -max_f64)
 }
 
-// NOTE: (joe-c) exponent notation is borked
+pub fn is_finite(f f64) bool {
+	return !is_nan(f) && !is_inf(f, 0)
+}
+
 // normalize returns a normal number y and exponent exp
 // satisfying x == y × 2**exp. It assumes x is finite and non-zero.
-// pub fn normalize(x f64) (f64, int) {
-// smallest_normal := 2.2250738585072014e-308 // 2**-1022
-// if abs(x) < smallest_normal {
-// return x * (1 << 52), -52
-// }
-// return x, 0
-// }
+pub fn normalize(x f64) (f64, int) {
+	smallest_normal := 2.2250738585072014e-308 // 2**-1022
+	if abs(x) < smallest_normal {
+		return x * math.normalize_smallest_mask, -52
+	}
+	return x, 0
+}
diff --git a/vlib/math/cbrt.c.v b/vlib/math/cbrt.c.v
new file mode 100644
index 000000000..892075a5e
--- /dev/null
+++ b/vlib/math/cbrt.c.v
@@ -0,0 +1,9 @@
+module math
+
+fn C.cbrt(x f64) f64
+
+// cbrt calculates cubic root.
+[inline]
+pub fn cbrt(a f64) f64 {
+	return C.cbrt(a)
+}
diff --git a/vlib/math/cbrt.js.v b/vlib/math/cbrt.js.v
new file mode 100644
index 000000000..306bba224
--- /dev/null
+++ b/vlib/math/cbrt.js.v
@@ -0,0 +1,9 @@
+module math
+
+fn JS.Math.cbrt(x f64) f64
+
+// cbrt calculates cubic root.
+[inline]
+pub fn cbrt(a f64) f64 {
+	return JS.Math.cbrt(a)
+}
diff --git a/vlib/math/cbrt.v b/vlib/math/cbrt.v
new file mode 100644
index 000000000..2c34ef2c3
--- /dev/null
+++ b/vlib/math/cbrt.v
@@ -0,0 +1,52 @@
+module math
+
+// cbrt returns the cube root of a.
+//
+// special cases are:
+// cbrt(±0) = ±0
+// cbrt(±inf) = ±inf
+// cbrt(nan) = nan
+pub fn cbrt(a f64) f64 {
+	mut x := a
+	b1 := 715094163 // (682-0.03306235651)*2**20
+	b2 := 696219795 // (664-0.03306235651)*2**20
+	c := 5.42857142857142815906e-01 // 19/35     = 0x3FE15F15F15F15F1
+	d := -7.05306122448979611050e-01 // -864/1225 = 0xBFE691DE2532C834
+	e_ := 1.41428571428571436819e+00 // 99/70     = 0x3FF6A0EA0EA0EA0F
+	f := 1.60714285714285720630e+00 // 45/28     = 0x3FF9B6DB6DB6DB6E
+	g := 3.57142857142857150787e-01 // 5/14      = 0x3FD6DB6DB6DB6DB7
+	smallest_normal := 2.22507385850720138309e-308 // 2**-1022  = 0x0010000000000000
+	if x == 0.0 || is_nan(x) || is_inf(x, 0) {
+		return x
+	}
+	mut sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+	// rough cbrt to 5 bits
+	mut t := f64_from_bits(f64_bits(x) / u64(3 + (u64(b1) << 32)))
+	if x < smallest_normal {
+		// subnormal number
+		t = f64(u64(1) << 54) // set t= 2**54
+		t *= x
+		t = f64_from_bits(f64_bits(t) / u64(3 + (u64(b2) << 32)))
+	}
+	// new cbrt to 23 bits
+	mut r := t * t / x
+	mut s := c + r * t
+	t *= g + f / (s + e_ + d / s)
+	// chop to 22 bits, make larger than cbrt(x)
+	t = f64_from_bits(f64_bits(t) & (u64(0xffffffffc) << 28) + (u64(1) << 30))
+	// one step newton iteration to 53 bits with error less than 0.667ulps
+	s = t * t // t*t is exact
+	r = x / s
+	w := t + t
+	r = (r - t) / (w + r) // r-s is exact
+	t = t + t * r
+	// restore the sign bit
+	if sign {
+		t = -t
+	}
+	return t
+}
diff --git a/vlib/math/complex/complex.v b/vlib/math/complex/complex.v
index 9fd7cc8cb..b7ec6aa74 100644
--- a/vlib/math/complex/complex.v
+++ b/vlib/math/complex/complex.v
@@ -27,7 +27,7 @@ pub fn (c Complex) str() string {
 // Complex Modulus value
 // mod() and abs() return the same
 pub fn (c Complex) abs() f64 {
-	return C.hypot(c.re, c.im)
+	return math.hypot(c.re, c.im)
 }
 
 pub fn (c Complex) mod() f64 {
diff --git a/vlib/math/const.v b/vlib/math/const.v
index 5a831a97c..7c473f79f 100644
--- a/vlib/math/const.v
+++ b/vlib/math/const.v
@@ -6,6 +6,8 @@ module math
 pub const (
 	e        = 2.71828182845904523536028747135266249775724709369995957496696763
 	pi       = 3.14159265358979323846264338327950288419716939937510582097494459
+	pi_2     = pi / 2.0
+	pi_4     = pi / 4.0
 	phi      = 1.61803398874989484820458683436563811772030917980576286213544862
 	tau      = 6.28318530717958647692528676655900576839433879875021164194988918
 	sqrt2    = 1.41421356237309504880168872420969807856967187537694807317667974
diff --git a/vlib/math/div.c.v b/vlib/math/div.c.v
new file mode 100644
index 000000000..90cee1af9
--- /dev/null
+++ b/vlib/math/div.c.v
@@ -0,0 +1,9 @@
+module math
+
+fn C.fmod(x f64, y f64) f64
+
+// fmod returns the floating-point remainder of number / denom (rounded towards zero):
+[inline]
+pub fn fmod(x f64, y f64) f64 {
+	return C.fmod(x, y)
+}
diff --git a/vlib/math/div.v b/vlib/math/div.v
new file mode 100644
index 000000000..6ad3c0d81
--- /dev/null
+++ b/vlib/math/div.v
@@ -0,0 +1,87 @@
+module math
+
+// Floating-point mod function.
+// mod returns the floating-point remainder of x/y.
+// The magnitude of the result is less than y and its
+// sign agrees with that of x.
+//
+// special cases are:
+// mod(±inf, y) = nan
+// mod(nan, y) = nan
+// mod(x, 0) = nan
+// mod(x, ±inf) = x
+// mod(x, nan) = nan
+pub fn mod(x f64, y f64) f64 {
+	return fmod(x, y)
+}
+
+// fmod returns the floating-point remainder of number / denom (rounded towards zero)
+pub fn fmod(x f64, y f64) f64 {
+	if y == 0 || is_inf(x, 0) || is_nan(x) || is_nan(y) {
+		return nan()
+	}
+	abs_y := abs(y)
+	abs_y_fr, abs_y_exp := frexp(abs_y)
+	mut r := x
+	if x < 0 {
+		r = -x
+	}
+	for r >= abs_y {
+		rfr, mut rexp := frexp(r)
+		if rfr < abs_y_fr {
+			rexp = rexp - 1
+		}
+		r = r - ldexp(abs_y, rexp - abs_y_exp)
+	}
+	if x < 0 {
+		r = -r
+	}
+	return r
+}
+
+// gcd calculates greatest common (positive) divisor (or zero if a and b are both zero).
+pub fn gcd(a_ i64, b_ i64) i64 {
+	mut a := a_
+	mut b := b_
+	if a < 0 {
+		a = -a
+	}
+	if b < 0 {
+		b = -b
+	}
+	for b != 0 {
+		a %= b
+		if a == 0 {
+			return b
+		}
+		b %= a
+	}
+	return a
+}
+
+// egcd returns (gcd(a, b), x, y) such that |a*x + b*y| = gcd(a, b)
+pub fn egcd(a i64, b i64) (i64, i64, i64) {
+	mut old_r, mut r := a, b
+	mut old_s, mut s := i64(1), i64(0)
+	mut old_t, mut t := i64(0), i64(1)
+
+	for r != 0 {
+		quot := old_r / r
+		old_r, r = r, old_r % r
+		old_s, s = s, old_s - quot * s
+		old_t, t = t, old_t - quot * t
+	}
+	return if old_r < 0 { -old_r } else { old_r }, old_s, old_t
+}
+
+// lcm calculates least common (non-negative) multiple.
+pub fn lcm(a i64, b i64) i64 {
+	if a == 0 {
+		return a
+	}
+	res := a * (b / gcd(b, a))
+	if res < 0 {
+		return -res
+	}
+	return res
+}
diff --git a/vlib/math/erf.c.v b/vlib/math/erf.c.v
new file mode 100644
index 000000000..160615210
--- /dev/null
+++ b/vlib/math/erf.c.v
@@ -0,0 +1,17 @@
+module math
+
+fn C.erf(x f64) f64
+
+fn C.erfc(x f64) f64
+
+// erf computes the error function value
+[inline]
+pub fn erf(a f64) f64 {
+	return C.erf(a)
+}
+
+// erfc computes the complementary error function value
+[inline]
+pub fn erfc(a f64) f64 {
+	return C.erfc(a)
+}
diff --git a/vlib/math/erf.v b/vlib/math/erf.v
index 043da317e..23757891b 100644
--- a/vlib/math/erf.v
+++ b/vlib/math/erf.v
@@ -1,659 +1,327 @@
-// Provides the [error](https://en.wikipedia.org/wiki/Error_function) and related functions
-// based on https://github.com/unovor/frame/blob/master/statrs-0.10.0/src/function/erf.rs
-//
-// NOTE: This impl does not have the same precision as glibc impl of erf,erfc and others, we should fix this
-// in the future.
 module math
 
-// Coefficients for erf_impl polynominal
+/*
+*                           x
+ *                    2      |\
+ *     erf(x)  =  ---------  | exp(-t*t)dt
+ *                 sqrt(pi) \|
+ *                           0
+ *
+ *     erfc(x) =  1-erf(x)
+ *  Note that
+ *              erf(-x) = -erf(x)
+ *              erfc(-x) = 2 - erfc(x)
+ *
+ * Method:
+ *      1. For |x| in [0, 0.84375]
+ *          erf(x)  = x + x*R(x**2)
+ *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
+ *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
+ *         where R = P/Q where P is an odd poly of degree 8 and
+ *         Q is an odd poly of degree 10.
+ *                                               -57.90
+ *                      | R - (erf(x)-x)/x | <= 2
+ *
+ *
+ *         Remark. The formula is derived by noting
+ *          erf(x) = (2/sqrt(pi))*(x - x**3/3 + x**5/10 - x**7/42 + ....)
+ *         and that
+ *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
+ *         is close to one. The interval is chosen because the fix
+ *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
+ *         near 0.6174), and by some experiment, 0.84375 is chosen to
+ *         guarantee the error is less than one ulp for erf.
+ *
+ *      2. For |x| in [0.84375,1.25], let s_ = |x| - 1, and
+ *         c = 0.84506291151 rounded to single (24 bits)
+ *              erf(x)  = sign(x) * (c  + P1(s_)/Q1(s_))
+ *              erfc(x) = (1-c)  - P1(s_)/Q1(s_) if x > 0
+ *                        1+(c+P1(s_)/Q1(s_))    if x < 0
+ *              |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
+ *         Remark: here we use the taylor series expansion at x=1.
+ *              erf(1+s_) = erf(1) + s_*Poly(s_)
+ *                       = 0.845.. + P1(s_)/Q1(s_)
+ *         That is, we use rational approximation to approximate
+ *                      erf(1+s_) - (c = (single)0.84506291151)
+ *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
+ *         where
+ *              P1(s_) = degree 6 poly in s_
+ *              Q1(s_) = degree 6 poly in s_
+ *
+ *      3. For x in [1.25,1/0.35(~2.857143)],
+ *              erfc(x) = (1/x)*exp(-x*x-0.5625+R1/s1)
+ *              erf(x)  = 1 - erfc(x)
+ *         where
+ *              R1(z) = degree 7 poly in z, (z=1/x**2)
+ *              s1(z) = degree 8 poly in z
+ *
+ *      4. For x in [1/0.35,28]
+ *              erfc(x) = (1/x)*exp(-x*x-0.5625+R2/s2) if x > 0
+ *                      = 2.0 - (1/x)*exp(-x*x-0.5625+R2/s2) if -6<x<0
+ *                      = 2.0 - tiny            (if x <= -6)
+ *              erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
+ *              erf(x)  = sign(x)*(1.0 - tiny)
+ *         where
+ *              R2(z) = degree 6 poly in z, (z=1/x**2)
+ *              s2(z) = degree 7 poly in z
+ *
+ *      Note1:
+ *         To compute exp(-x*x-0.5625+R/s), let s_ be a single
+ *         precision number and s_ := x; then
+ *              -x*x = -s_*s_ + (s_-x)*(s_+x)
+ *              exp(-x*x-0.5626+R/s) =
+ *                      exp(-s_*s_-0.5625)*exp((s_-x)*(s_+x)+R/s);
+ *      Note2:
+ *         Here 4 and 5 make use of the asymptotic series
+ *                        exp(-x*x)
+ *              erfc(x) ~ ---------- * ( 1 + Poly(1/x**2) )
+ *                        x*sqrt(pi)
+ *         We use rational approximation to approximate
+ *              g(s_)=f(1/x**2) = log(erfc(x)*x) - x*x + 0.5625
+ *         Here is the error bound for R1/s1 and R2/s2
+ *              |R1/s1 - f(x)|  < 2**(-62.57)
+ *              |R2/s2 - f(x)|  < 2**(-61.52)
+ *
+ *      5. For inf > x >= 28
+ *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
+ *              erfc(x) = tiny*tiny (raise underflow) if x > 0
+ *                      = 2 - tiny if x<0
+ *
+ *      7. special case:
+ *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
+ *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
+ *              erfc/erf(nan) is nan
+*/
 const (
-	// Polynomial coefficients for a numerator of `erf_impl`
-	// in the interval [1e-10, 0.5].
-	erf_impl_an     = [0.00337916709551257388990745, -0.00073695653048167948530905,
-		-0.374732337392919607868241, 0.0817442448733587196071743, -0.0421089319936548595203468,
-		0.0070165709512095756344528, -0.00495091255982435110337458, 0.000871646599037922480317225]
-	// Polynomial coefficients for a denominator of `erf_impl`
-	// in the interval [1e-10, 0.5]
-	erf_impl_ad     = [1.0, -0.218088218087924645390535, 0.412542972725442099083918,
-		-0.0841891147873106755410271, 0.0655338856400241519690695, -0.0120019604454941768171266,
-		0.00408165558926174048329689, -0.000615900721557769691924509]
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [0.5, 0.75].
-	erf_impl_bn     = [-0.0361790390718262471360258, 0.292251883444882683221149,
-		0.281447041797604512774415, 0.125610208862766947294894, 0.0274135028268930549240776,
-		0.00250839672168065762786937,
-	]
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [0.5, 0.75].
-	erf_impl_bd     = [1.0, 1.8545005897903486499845, 1.43575803037831418074962,
-		0.582827658753036572454135, 0.124810476932949746447682, 0.0113724176546353285778481]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [0.75, 1.25].
-	erf_impl_cn     = [
-		-0.0397876892611136856954425,
-		0.153165212467878293257683,
-		0.191260295600936245503129,
-		0.10276327061989304213645,
-		0.029637090615738836726027,
-		0.0046093486780275489468812,
-		0.000307607820348680180548455,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [0.75, 1.25].
-	erf_impl_cd     = [
-		1.0,
-		1.95520072987627704987886,
-		1.64762317199384860109595,
-		0.768238607022126250082483,
-		0.209793185936509782784315,
-		0.0319569316899913392596356,
-		0.00213363160895785378615014,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [1.25, 2.25].
-	erf_impl_dn     = [
-		-0.0300838560557949717328341,
-		0.0538578829844454508530552,
-		0.0726211541651914182692959,
-		0.0367628469888049348429018,
-		0.00964629015572527529605267,
-		0.00133453480075291076745275,
-		0.778087599782504251917881e-4,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [1.25, 2.25].
-	erf_impl_dd     = [
-		1.0,
-		1.75967098147167528287343,
-		1.32883571437961120556307,
-		0.552528596508757581287907,
-		0.133793056941332861912279,
-		0.0179509645176280768640766,
-		0.00104712440019937356634038,
-		-0.106640381820357337177643e-7,
-	]
-
-	//  Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [2.25, 3.5].
-	erf_impl_en     = [
-		-0.0117907570137227847827732,
-		0.014262132090538809896674,
-		0.0202234435902960820020765,
-		0.00930668299990432009042239,
-		0.00213357802422065994322516,
-		0.00025022987386460102395382,
-		0.120534912219588189822126e-4,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [2.25, 3.5].
-	erf_impl_ed     = [
-		1.0,
-		1.50376225203620482047419,
-		0.965397786204462896346934,
-		0.339265230476796681555511,
-		0.0689740649541569716897427,
-		0.00771060262491768307365526,
-		0.000371421101531069302990367,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [3.5, 5.25].
-	erf_impl_fn     = [
-		-0.00546954795538729307482955,
-		0.00404190278731707110245394,
-		0.0054963369553161170521356,
-		0.00212616472603945399437862,
-		0.000394984014495083900689956,
-		0.365565477064442377259271e-4,
-		0.135485897109932323253786e-5,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [3.5, 5.25].
-	erf_impl_fd     = [
-		1.0,
-		1.21019697773630784832251,
-		0.620914668221143886601045,
-		0.173038430661142762569515,
-		0.0276550813773432047594539,
-		0.00240625974424309709745382,
-		0.891811817251336577241006e-4,
-		-0.465528836283382684461025e-11,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [5.25, 8].
-	erf_impl_gn     = [
-		-0.00270722535905778347999196,
-		0.0013187563425029400461378,
-		0.00119925933261002333923989,
-		0.00027849619811344664248235,
-		0.267822988218331849989363e-4,
-		0.923043672315028197865066e-6,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [5.25, 8].
-	erf_impl_gd     = [
-		1.0,
-		0.814632808543141591118279,
-		0.268901665856299542168425,
-		0.0449877216103041118694989,
-		0.00381759663320248459168994,
-		0.000131571897888596914350697,
-		0.404815359675764138445257e-11,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [8, 11.5].
-	erf_impl_hn     = [
-		-0.00109946720691742196814323,
-		0.000406425442750422675169153,
-		0.000274499489416900707787024,
-		0.465293770646659383436343e-4,
-		0.320955425395767463401993e-5,
-		0.778286018145020892261936e-7,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [8, 11.5].
-	erf_impl_hd     = [
-		1.0,
-		0.588173710611846046373373,
-		0.139363331289409746077541,
-		0.0166329340417083678763028,
-		0.00100023921310234908642639,
-		0.24254837521587225125068e-4,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [11.5, 17].
-	erf_impl_in     = [
-		-0.00056907993601094962855594,
-		0.000169498540373762264416984,
-		0.518472354581100890120501e-4,
-		0.382819312231928859704678e-5,
-		0.824989931281894431781794e-7,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [11.5, 17].
-	erf_impl_id     = [
-		1.0,
-		0.339637250051139347430323,
-		0.043472647870310663055044,
-		0.00248549335224637114641629,
-		0.535633305337152900549536e-4,
-		-0.117490944405459578783846e-12,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [17, 24].
-	erf_impl_jn     = [
-		-0.000241313599483991337479091,
-		0.574224975202501512365975e-4,
-		0.115998962927383778460557e-4,
-		0.581762134402593739370875e-6,
-		0.853971555085673614607418e-8,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [17, 24].
-	erf_impl_jd     = [
-		1.0,
-		0.233044138299687841018015,
-		0.0204186940546440312625597,
-		0.000797185647564398289151125,
-		0.117019281670172327758019e-4,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [24, 38].
-	erf_impl_kn     = [
-		-0.000146674699277760365803642,
-		0.162666552112280519955647e-4,
-		0.269116248509165239294897e-5,
-		0.979584479468091935086972e-7,
-		0.101994647625723465722285e-8,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [24, 38].
-	erf_impl_kd     = [
-		1.0,
-		0.165907812944847226546036,
-		0.0103361716191505884359634,
-		0.000286593026373868366935721,
-		0.298401570840900340874568e-5,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [38, 60].
-	erf_impl_ln     = [
-		-0.583905797629771786720406e-4,
-		0.412510325105496173512992e-5,
-		0.431790922420250949096906e-6,
-		0.993365155590013193345569e-8,
-		0.653480510020104699270084e-10,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [38, 60].
-	erf_impl_ld     = [
-		1.0,
-		0.105077086072039915406159,
-		0.00414278428675475620830226,
-		0.726338754644523769144108e-4,
-		0.477818471047398785369849e-6,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [60, 85].
-	erf_impl_mn     = [
-		-0.196457797609229579459841e-4,
-		0.157243887666800692441195e-5,
-		0.543902511192700878690335e-7,
-		0.317472492369117710852685e-9,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [60, 85].
-	erf_impl_md     = [
-		1.0,
-		0.052803989240957632204885,
-		0.000926876069151753290378112,
-		0.541011723226630257077328e-5,
-		0.535093845803642394908747e-15,
-	]
-
-	// Polynomial coefficients for a numerator in `erf_impl`
-	// in the interval [85, 110].
-	erf_impl_nn     = [
-		-0.789224703978722689089794e-5,
-		0.622088451660986955124162e-6,
-		0.145728445676882396797184e-7,
-		0.603715505542715364529243e-10,
-	]
-
-	// Polynomial coefficients for a denominator in `erf_impl`
-	// in the interval [85, 110].
-	erf_impl_nd     = [
-		1.0,
-		0.0375328846356293715248719,
-		0.000467919535974625308126054,
-		0.193847039275845656900547e-5,
-	]
-
-	// **********************************************************
-	// ********** Coefficients for erf_inv_impl polynomial ******
-	// **********************************************************
-	// Polynomial coefficients for a numerator of `erf_inv_impl`
-	// in the interval [0, 0.5].
-	erf_inv_impl_an = [
-		-0.000508781949658280665617,
-		-0.00836874819741736770379,
-		0.0334806625409744615033,
-		-0.0126926147662974029034,
-		-0.0365637971411762664006,
-		0.0219878681111168899165,
-		0.00822687874676915743155,
-		-0.00538772965071242932965,
-	]
-
-	// Polynomial coefficients for a denominator of `erf_inv_impl`
-	// in the interval [0, 0.5].
-	erf_inv_impl_ad = [
-		1.0,
-		-0.970005043303290640362,
-		-1.56574558234175846809,
-		1.56221558398423026363,
-		0.662328840472002992063,
-		-0.71228902341542847553,
-		-0.0527396382340099713954,
-		0.0795283687341571680018,
-		-0.00233393759374190016776,
-		0.000886216390456424707504,
-	]
-
-	// Polynomial coefficients for a numerator of `erf_inv_impl`
-	// in the interval [0.5, 0.75].
-	erf_inv_impl_bn = [
-		-0.202433508355938759655,
-		0.105264680699391713268,
-		8.37050328343119927838,
-		17.6447298408374015486,
-		-18.8510648058714251895,
-		-44.6382324441786960818,
-		17.445385985570866523,
-		21.1294655448340526258,
-		-3.67192254707729348546,
-	]
-
-	// Polynomial coefficients for a denominator of `erf_inv_impl`
-	// in the interval [0.5, 0.75].
-	erf_inv_impl_bd = [
-		1.0,
-		6.24264124854247537712,
-		3.9713437953343869095,
-		-28.6608180499800029974,
-		-20.1432634680485188801,
-		48.5609213108739935468,
-		10.8268667355460159008,
-		-22.6436933413139721736,
-		1.72114765761200282724,
-	]
-
-	// Polynomial coefficients for a numerator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x less than 3.
-	erf_inv_impl_cn = [
-		-0.131102781679951906451,
-		-0.163794047193317060787,
-		0.117030156341995252019,
-		0.387079738972604337464,
-		0.337785538912035898924,
-		0.142869534408157156766,
-		0.0290157910005329060432,
-		0.00214558995388805277169,
-		-0.679465575181126350155e-6,
-		0.285225331782217055858e-7,
-		-0.681149956853776992068e-9,
-	]
-
-	// Polynomial coefficients for a denominator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x less than 3.
-	erf_inv_impl_cd = [
-		1.0,
-		3.46625407242567245975,
-		5.38168345707006855425,
-		4.77846592945843778382,
-		2.59301921623620271374,
-		0.848854343457902036425,
-		0.152264338295331783612,
-		0.01105924229346489121,
-	]
-
-	// Polynomial coefficients for a numerator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x between 3 and 6.
-	erf_inv_impl_dn = [
-		-0.0350353787183177984712,
-		-0.00222426529213447927281,
-		0.0185573306514231072324,
-		0.00950804701325919603619,
-		0.00187123492819559223345,
-		0.000157544617424960554631,
-		0.460469890584317994083e-5,
-		-0.230404776911882601748e-9,
-		0.266339227425782031962e-11,
-	]
-
-	// Polynomial coefficients for a denominator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x between 3 and 6.
-	erf_inv_impl_dd = [
-		1.0,
-		1.3653349817554063097,
-		0.762059164553623404043,
-		0.220091105764131249824,
-		0.0341589143670947727934,
-		0.00263861676657015992959,
-		0.764675292302794483503e-4,
-	]
-
-	// Polynomial coefficients for a numerator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x between 6 and 18.
-	erf_inv_impl_en = [
-		-0.0167431005076633737133,
-		-0.00112951438745580278863,
-		0.00105628862152492910091,
-		0.000209386317487588078668,
-		0.149624783758342370182e-4,
-		0.449696789927706453732e-6,
-		0.462596163522878599135e-8,
-		-0.281128735628831791805e-13,
-		0.99055709973310326855e-16,
-	]
-
-	// Polynomial coefficients for a denominator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x between 6 and 18.
-	erf_inv_impl_ed = [
-		1.0,
-		0.591429344886417493481,
-		0.138151865749083321638,
-		0.0160746087093676504695,
-		0.000964011807005165528527,
-		0.275335474764726041141e-4,
-		0.282243172016108031869e-6,
-	]
-
-	// Polynomial coefficients for a numerator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x between 18 and 44.
-	erf_inv_impl_fn = [
-		-0.0024978212791898131227,
-		-0.779190719229053954292e-5,
-		0.254723037413027451751e-4,
-		0.162397777342510920873e-5,
-		0.396341011304801168516e-7,
-		0.411632831190944208473e-9,
-		0.145596286718675035587e-11,
-		-0.116765012397184275695e-17,
-	]
-
-	// Polynomial coefficients for a denominator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x between 18 and 44.
-	erf_inv_impl_fd = [
-		1.0,
-		0.207123112214422517181,
-		0.0169410838120975906478,
-		0.000690538265622684595676,
-		0.145007359818232637924e-4,
-		0.144437756628144157666e-6,
-		0.509761276599778486139e-9,
-	]
-
-	// Polynomial coefficients for a numerator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x greater than 44.
-	erf_inv_impl_gn = [
-		-0.000539042911019078575891,
-		-0.28398759004727721098e-6,
-		0.899465114892291446442e-6,
-		0.229345859265920864296e-7,
-		0.225561444863500149219e-9,
-		0.947846627503022684216e-12,
-		0.135880130108924861008e-14,
-		-0.348890393399948882918e-21,
-	]
-
-	// Polynomial coefficients for a denominator of `erf_inv_impl`
-	// in the interval [0.75, 1] with x greater than 44.
-	erf_inv_impl_gd = [
-		1.0,
-		0.0845746234001899436914,
-		0.00282092984726264681981,
-		0.468292921940894236786e-4,
-		0.399968812193862100054e-6,
-		0.161809290887904476097e-8,
-		0.231558608310259605225e-11,
-	]
+	erx  = 8.45062911510467529297e-01 // 0x3FEB0AC160000000
+	// Coefficients for approximation to  erf in [0, 0.84375]
+	efx  = 1.28379167095512586316e-01 // 0x3FC06EBA8214DB69
+	efx8 = 1.02703333676410069053e+00 // 0x3FF06EBA8214DB69
+	pp0  = 1.28379167095512558561e-01 // 0x3FC06EBA8214DB68
+	pp1  = -3.25042107247001499370e-01 // 0xBFD4CD7D691CB913
+	pp2  = -2.84817495755985104766e-02 // 0xBF9D2A51DBD7194F
+	pp3  = -5.77027029648944159157e-03 // 0xBF77A291236668E4
+	pp4  = -2.37630166566501626084e-05 // 0xBEF8EAD6120016AC
+	qq1  = 3.97917223959155352819e-01 // 0x3FD97779CDDADC09
+	qq2  = 6.50222499887672944485e-02 // 0x3FB0A54C5536CEBA
+	qq3  = 5.08130628187576562776e-03 // 0x3F74D022C4D36B0F
+	qq4  = 1.32494738004321644526e-04 // 0x3F215DC9221C1A10
+	qq5  = -3.96022827877536812320e-06 // 0xBED09C4342A26120
+	// Coefficients for approximation to  erf  in [0.84375, 1.25]
+	pa0  = -2.36211856075265944077e-03 // 0xBF6359B8BEF77538
+	pa1  = 4.14856118683748331666e-01 // 0x3FDA8D00AD92B34D
+	pa2  = -3.72207876035701323847e-01 // 0xBFD7D240FBB8C3F1
+	pa3  = 3.18346619901161753674e-01 // 0x3FD45FCA805120E4
+	pa4  = -1.10894694282396677476e-01 // 0xBFBC63983D3E28EC
+	pa5  = 3.54783043256182359371e-02 // 0x3FA22A36599795EB
+	pa6  = -2.16637559486879084300e-03 // 0xBF61BF380A96073F
+	qa1  = 1.06420880400844228286e-01 // 0x3FBB3E6618EEE323
+	qa2  = 5.40397917702171048937e-01 // 0x3FE14AF092EB6F33
+	qa3  = 7.18286544141962662868e-02 // 0x3FB2635CD99FE9A7
+	qa4  = 1.26171219808761642112e-01 // 0x3FC02660E763351F
+	qa5  = 1.36370839120290507362e-02 // 0x3F8BEDC26B51DD1C
+	qa6  = 1.19844998467991074170e-02 // 0x3F888B545735151D
+	// Coefficients for approximation to  erfc in [1.25, 1/0.35]
+	ra0  = -9.86494403484714822705e-03 // 0xBF843412600D6435
+	ra1  = -6.93858572707181764372e-01 // 0xBFE63416E4BA7360
+	ra2  = -1.05586262253232909814e+01 // 0xC0251E0441B0E726
+	ra3  = -6.23753324503260060396e+01 // 0xC04F300AE4CBA38D
+	ra4  = -1.62396669462573470355e+02 // 0xC0644CB184282266
+	ra5  = -1.84605092906711035994e+02 // 0xC067135CEBCCABB2
+	ra6  = -8.12874355063065934246e+01 // 0xC054526557E4D2F2
+	ra7  = -9.81432934416914548592e+00 // 0xC023A0EFC69AC25C
+	sa1  = 1.96512716674392571292e+01 // 0x4033A6B9BD707687
+	sa2  = 1.37657754143519042600e+02 // 0x4061350C526AE721
+	sa3  = 4.34565877475229228821e+02 // 0x407B290DD58A1A71
+	sa4  = 6.45387271733267880336e+02 // 0x40842B1921EC2868
+	sa5  = 4.29008140027567833386e+02 // 0x407AD02157700314
+	sa6  = 1.08635005541779435134e+02 // 0x405B28A3EE48AE2C
+	sa7  = 6.57024977031928170135e+00 // 0x401A47EF8E484A93
+	sa8  = -6.04244152148580987438e-02 // 0xBFAEEFF2EE749A62
+	// Coefficients for approximation to  erfc in [1/.35, 28]
+	rb0  = -9.86494292470009928597e-03 // 0xBF84341239E86F4A
+	rb1  = -7.99283237680523006574e-01 // 0xBFE993BA70C285DE
+	rb2  = -1.77579549177547519889e+01 // 0xC031C209555F995A
+	rb3  = -1.60636384855821916062e+02 // 0xC064145D43C5ED98
+	rb4  = -6.37566443368389627722e+02 // 0xC083EC881375F228
+	rb5  = -1.02509513161107724954e+03 // 0xC09004616A2E5992
+	rb6  = -4.83519191608651397019e+02 // 0xC07E384E9BDC383F
+	sb1  = 3.03380607434824582924e+01 // 0x403E568B261D5190
+	sb2  = 3.25792512996573918826e+02 // 0x40745CAE221B9F0A
+	sb3  = 1.53672958608443695994e+03 // 0x409802EB189D5118
+	sb4  = 3.19985821950859553908e+03 // 0x40A8FFB7688C246A
+	sb5  = 2.55305040643316442583e+03 // 0x40A3F219CEDF3BE6
+	sb6  = 4.74528541206955367215e+02 // 0x407DA874E79FE763
+	sb7  = -2.24409524465858183362e+01 // 0xC03670E242712D62
 )
 
-fn erf_inv_impl(p f64, q f64, s f64) f64 {
-	mut result := 0.0
-	if p <= 0.5 {
-		y := 0.0891314744949340820313
-		g := p * (p + 10.0)
-		r := polynomial(p, math.erf_inv_impl_an) / polynomial(p, math.erf_inv_impl_ad)
-		result = g * y + g * r
-	} else if q >= 0.25 {
-		y := 2.249481201171875
-		g := sqrt(-2.0 * log(q))
-		xs := q - 0.25
-		r := polynomial(xs, math.erf_inv_impl_bn) / polynomial(xs, math.erf_inv_impl_bd)
-		result = g / (y + r)
-	} else {
-		x := sqrt(-log(q))
-		if x < 3.0 {
-			y := 0.807220458984375
-			xs := x - 1.125
-			r := polynomial(xs, math.erf_inv_impl_cn) / polynomial(xs, math.erf_inv_impl_cd)
-			result = y * x + r * x
-		} else if x < 6.0 {
-			y := 0.93995571136474609375
-			xs := x - 3.0
-			r := polynomial(xs, math.erf_inv_impl_dn) / polynomial(xs, math.erf_inv_impl_dd)
-			result = y * x + r * x
-		} else if x < 18.0 {
-			y := 0.98362827301025390625
-			xs := x - 6.0
-			r := polynomial(xs, math.erf_inv_impl_en) / polynomial(xs, math.erf_inv_impl_ed)
-			result = y * x + r * x
-		} else if x < 44.0 {
-			y := 0.99714565277099609375
-			xs := x - 18.0
-			r := polynomial(xs, math.erf_inv_impl_fn) / polynomial(xs, math.erf_inv_impl_fd)
-			result = y * x + r * x
-		} else {
-			y := 0.99941349029541015625
-			xs := x - 44.0
-			r := polynomial(xs, math.erf_inv_impl_gn) / polynomial(xs, math.erf_inv_impl_gd)
-			result = y * x + r * x
-		}
+// erf returns the error function of x.
+//
+// special cases are:
+// erf(+inf) = 1
+// erf(-inf) = -1
+// erf(nan) = nan
+pub fn erf(a f64) f64 {
+	mut x := a
+	very_tiny := 2.848094538889218e-306 // 0x0080000000000000
+	small := 1.0 / f64(u64(1) << 28) // 2**-28
+	if is_nan(x) {
+		return nan()
 	}
-
-	return s * result
-}
-
-fn erf_impl(z f64, inv bool) f64 {
-	if z < 0.0 {
-		if !inv {
-			return -erf_impl(-z, false)
-		}
-		if z < -0.5 {
-			return 2.0 - erf_impl(-z, true)
-		}
-		return 1.0 + erf_impl(-z, false)
+	if is_inf(x, 1) {
+		return 1.0
+	}
+	if is_inf(x, -1) {
+		return f64(-1)
 	}
-	mut result := 0.0
-	if z < 0.5 {
-		if z < 1e-10 {
-			result = z * 1.125 + z * 0.003379167095512573896158903121545171688
+	mut sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+	if x < 0.84375 { // |x| < 0.84375
+		mut temp := 0.0
+		if x < small { // |x| < 2**-28
+			if x < very_tiny {
+				temp = 0.125 * (8.0 * x + math.efx8 * x) // avoid underflow
+			} else {
+				temp = x + math.efx * x
+			}
 		} else {
-			result = z * 1.125 +
-				z * polynomial(z, math.erf_impl_an) / polynomial(z, math.erf_impl_ad)
+			z := x * x
+			r := math.pp0 + z * (math.pp1 + z * (math.pp2 + z * (math.pp3 + z * math.pp4)))
+			s_ := 1.0 + z * (math.qq1 + z * (math.qq2 + z * (math.qq3 + z * (math.qq4 +
+				z * math.qq5))))
+			y := r / s_
+			temp = x + x * y
 		}
-	} else if z < 110.0 {
-		mut r := 0.0
-		mut b := 0.0
-		if z < 0.75 {
-			r = polynomial(z - 0.5, math.erf_impl_bn) / polynomial(z - 0.5, math.erf_impl_bd)
-			b = 0.3440242112
-		} else if z < 1.25 {
-			r = polynomial(z - 0.75, math.erf_impl_cn) / polynomial(z - 0.75, math.erf_impl_cd)
-			b = 0.419990927
-		} else if z < 2.25 {
-			r = polynomial(z - 1.25, math.erf_impl_dn) / polynomial(z - 1.25, math.erf_impl_dd)
-			b = 0.4898625016
-		} else if z < 3.5 {
-			r = polynomial(z - 2.25, math.erf_impl_en) / polynomial(z - 2.25, math.erf_impl_ed)
-			b = 0.5317370892
-		} else if z < 5.25 {
-			r = polynomial(z - 3.5, math.erf_impl_fn) / polynomial(z - 3.5, math.erf_impl_fd)
-			b = 0.5489973426
-		} else if z < 8.0 {
-			r = polynomial(z - 5.25, math.erf_impl_gn) / polynomial(z - 5.25, math.erf_impl_gd)
-			b = 0.5571740866
-		} else if z < 11.5 {
-			r = polynomial(z - 8.0, math.erf_impl_hn) / polynomial(z - 8.0, math.erf_impl_hd)
-			b = 0.5609807968
-		} else if z < 17.0 {
-			r = polynomial(z - 11.5, math.erf_impl_in) / polynomial(z - 11.5, math.erf_impl_id)
-			b = 0.5626493692
-		} else if z < 24.0 {
-			r = polynomial(z - 17.0, math.erf_impl_jn) / polynomial(z - 17.0, math.erf_impl_jd)
-			b = 0.5634598136
-		} else if z < 38.0 {
-			r = polynomial(z - 24.0, math.erf_impl_kn) / polynomial(z - 24.0, math.erf_impl_kd)
-			b = 0.5638477802
-		} else if z < 60.0 {
-			r = polynomial(z - 38.0, math.erf_impl_ln) / polynomial(z - 38.0, math.erf_impl_ld)
-			b = 0.5640528202
-		} else if z < 85.0 {
-			r = polynomial(z - 60.0, math.erf_impl_mn) / polynomial(z - 60.0, math.erf_impl_md)
-			b = 0.5641309023
-		} else {
-			r = polynomial(z - 85.0, math.erf_impl_nn) / polynomial(z - 85.0, math.erf_impl_nd)
-			b = 0.5641584396
+		if sign {
+			return -temp
 		}
-
-		g := exp(-z * z) / z
-		result = g * b + g * r
-	} else {
-		result = 0.0
+		return temp
 	}
-	if inv && z >= 0.5 {
-		return result
-	} else if z >= 0.5 || inv {
-		return 1.0 - result
-	} else {
-		return result
+	if x < 1.25 { // 0.84375 <= |x| < 1.25
+		s_ := x - 1
+		p := math.pa0 + s_ * (math.pa1 + s_ * (math.pa2 + s_ * (math.pa3 + s_ * (math.pa4 +
+			s_ * (math.pa5 + s_ * math.pa6)))))
+		q := 1.0 + s_ * (math.qa1 + s_ * (math.qa2 + s_ * (math.qa3 + s_ * (math.qa4 +
+			s_ * (math.qa5 + s_ * math.qa6)))))
+		if sign {
+			return -math.erx - p / q
+		}
+		return math.erx + p / q
 	}
-}
-
-/// 'erf' calculates the error function at `x`.
-pub fn erf(x f64) f64 {
-	if is_nan(x) {
-		return nan()
-	} else if is_inf(x, 1) {
+	if x >= 6 { // inf > |x| >= 6
+		if sign {
+			return -1
+		}
 		return 1.0
-	} else if is_inf(x, -1) {
-		return -1.0
-	} else if x == 0.0 {
-		return 0.0
-	} else {
-		return erf_impl(x, false)
 	}
-}
-
-// `erf_inv` calculates the inverse error function at `x`.
-pub fn erf_inv(x f64) f64 {
-	if x == 0 {
-		return 0.0
-	} else if x >= 1.0 {
-		return inf(1)
-	} else if x <= -1.0 {
-		return inf(-1)
-	} else if x < 0.0 {
-		return erf_inv_impl(-x, 1.0 + x, -1.0)
-	} else {
-		return erf_inv_impl(x, 1.0 - x, 1.0)
+	s_ := 1.0 / (x * x)
+	mut r := 0.0
+	mut s := 0.0
+	if x < 1.0 / 0.35 { // |x| < 1 / 0.35  ~ 2.857143
+		r = math.ra0 + s_ * (math.ra1 + s_ * (math.ra2 + s_ * (math.ra3 + s_ * (math.ra4 +
+			s_ * (math.ra5 + s_ * (math.ra6 + s_ * math.ra7))))))
+		s = 1.0 + s_ * (math.sa1 + s_ * (math.sa2 + s_ * (math.sa3 + s_ * (math.sa4 +
+			s_ * (math.sa5 + s_ * (math.sa6 + s_ * (math.sa7 + s_ * math.sa8)))))))
+	} else { // |x| >= 1 / 0.35  ~ 2.857143
+		r = math.rb0 + s_ * (math.rb1 + s_ * (math.rb2 + s_ * (math.rb3 + s_ * (math.rb4 +
+			s_ * (math.rb5 + s_ * math.rb6)))))
+		s = 1.0 + s_ * (math.sb1 + s_ * (math.sb2 + s_ * (math.sb3 + s_ * (math.sb4 +
+			s_ * (math.sb5 + s_ * (math.sb6 + s_ * math.sb7))))))
+	}
+	z := f64_from_bits(f64_bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precision x
+	r_ := exp(-z * z - 0.5625) * exp((z - x) * (z + x) + r / s)
+	if sign {
+		return r_ / x - 1.0
 	}
+	return 1.0 - r_ / x
 }
 
-// `erfc` calculates the complementary error function at `x`.
-pub fn erfc(x f64) f64 {
+// erfc returns the complementary error function of x.
+//
+// special cases are:
+// erfc(+inf) = 0
+// erfc(-inf) = 2
+// erfc(nan) = nan
+pub fn erfc(a f64) f64 {
+	mut x := a
+	tiny := 1.0 / f64(u64(1) << 56) // 2**-56
+	// special cases
 	if is_nan(x) {
 		return nan()
-	} else if is_inf(x, 1) {
+	}
+	if is_inf(x, 1) {
 		return 0.0
-	} else if is_inf(x, -1) {
+	}
+	if is_inf(x, -1) {
 		return 2.0
-	} else {
-		return erf_impl(x, true)
 	}
-}
-
-// `erfc_inv` calculates the complementary inverse error function at `x`.
-pub fn erfc_inv(x f64) f64 {
-	if x <= 0.0 {
-		return inf(1)
-	} else if x >= 2.0 {
-		return inf(-1)
-	} else if is_inf(x, -1) {
-		return erf_inv_impl(-1.0 + x, 2.0 - x, -1.0)
-	} else {
-		return erf_inv_impl(1.0 - x, x, 1.0)
+	mut sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+	if x < 0.84375 { // |x| < 0.84375
+		mut temp := 0.0
+		if x < tiny { // |x| < 2**-56
+			temp = x
+		} else {
+			z := x * x
+			r := math.pp0 + z * (math.pp1 + z * (math.pp2 + z * (math.pp3 + z * math.pp4)))
+			s_ := 1.0 + z * (math.qq1 + z * (math.qq2 + z * (math.qq3 + z * (math.qq4 +
+				z * math.qq5))))
+			y := r / s_
+			if x < 0.25 { // |x| < 1.0/4
+				temp = x + x * y
+			} else {
+				temp = 0.5 + (x * y + (x - 0.5))
+			}
+		}
+		if sign {
+			return 1.0 + temp
+		}
+		return 1.0 - temp
+	}
+	if x < 1.25 { // 0.84375 <= |x| < 1.25
+		s_ := x - 1
+		p := math.pa0 + s_ * (math.pa1 + s_ * (math.pa2 + s_ * (math.pa3 + s_ * (math.pa4 +
+			s_ * (math.pa5 + s_ * math.pa6)))))
+		q := 1.0 + s_ * (math.qa1 + s_ * (math.qa2 + s_ * (math.qa3 + s_ * (math.qa4 +
+			s_ * (math.qa5 + s_ * math.qa6)))))
+		if sign {
+			return 1.0 + math.erx + p / q
+		}
+		return 1.0 - math.erx - p / q
+	}
+	if x < 28 { // |x| < 28
+		s_ := 1.0 / (x * x)
+		mut r := 0.0
+		mut s := 0.0
+		if x < 1.0 / 0.35 { // |x| < 1 / 0.35 ~ 2.857143
+			r = math.ra0 + s_ * (math.ra1 + s_ * (math.ra2 + s_ * (math.ra3 + s_ * (math.ra4 +
+				s_ * (math.ra5 + s_ * (math.ra6 + s_ * math.ra7))))))
+			s = 1.0 + s_ * (math.sa1 + s_ * (math.sa2 + s_ * (math.sa3 + s_ * (math.sa4 +
+				s_ * (math.sa5 + s_ * (math.sa6 + s_ * (math.sa7 + s_ * math.sa8)))))))
+		} else { // |x| >= 1 / 0.35 ~ 2.857143
+			if sign && x > 6 {
+				return 2.0 // x < -6
+			}
+			r = math.rb0 + s_ * (math.rb1 + s_ * (math.rb2 + s_ * (math.rb3 + s_ * (math.rb4 +
+				s_ * (math.rb5 + s_ * math.rb6)))))
+			s = 1.0 + s_ * (math.sb1 + s_ * (math.sb2 + s_ * (math.sb3 + s_ * (math.sb4 +
+				s_ * (math.sb5 + s_ * (math.sb6 + s_ * math.sb7))))))
+		}
+		z := f64_from_bits(f64_bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precision x
+		r_ := exp(-z * z - 0.5625) * exp((z - x) * (z + x) + r / s)
+		if sign {
+			return 2.0 - r_ / x
+		}
+		return r_ / x
+	}
+	if sign {
+		return 2.0
 	}
+	return 0.0
 }
diff --git a/vlib/math/erf_test.v b/vlib/math/erf_test.v
index 00cedb916..2102514de 100644
--- a/vlib/math/erf_test.v
+++ b/vlib/math/erf_test.v
@@ -1,45 +1,29 @@
-import math
-
-fn almost_eq(a f64, b f64, acc f64) bool {
-	if math.is_inf(a, 1) || math.is_inf(a, -1) || math.is_inf(b, 1) || math.is_inf(b, -1) {
-		return a == b
-	}
-
-	if math.is_nan(a) || math.is_nan(b) {
-		return false
-	}
-
-	return math.abs(a - b) < acc
-}
+module math
 
 fn test_erf() {
-	assert math.is_nan(math.erf(math.nan()))
-	assert almost_eq(math.erf(-1.0), -0.8427007888650501, 1e-11)
-	assert math.erf(0.0) == 0.0
-	assert math.erf(1e-15) == 0.0000000000000011283791670955126615773132947717431253912942469337536
-	assert math.erf(0.1) == 0.11246291601917208
-	assert math.erf(0.3) == 0.32862677677789676
-	assert almost_eq(math.erf(0.5), 0.5204998778130465376827466538919645287364515757579637,
+	assert is_nan(erf(nan()))
+	assert tolerance(erf(-1.0), -0.8427007888650501, 1e-8)
+	assert tolerance(erf(0.0), 0.0, 1e-11)
+	assert tolerance(erf(1e-15), 0.0000000000000011283791670955126615773132947717431253912942469337536,
+		1e-11)
+	assert tolerance(erf(0.1), 0.11246291601917208, 1e-11)
+	assert tolerance(erf(0.3), 0.32862677677789676, 1e-7)
+	assert tolerance(erf(0.5), 0.5204998778130465376827466538919645287364515757579637,
 		1e-9)
-	assert math.erf(1.0) == 0.8427007888650501
-	assert math.erf(1.5) == 0.966105146259005
-	assert math.erf(6.0) == 0.99999999999999997848026328750108688340664960081261537
-	assert math.erf(5.0) == 0.99999999999846254020557196514981165651461662110988195
-	assert math.erf(4.0) == 0.999999984582742
-	assert math.erf(math.inf(1)) == 1.0
-	assert math.erf(math.inf(-1)) == -1.0
+	assert tolerance(erf(1.0), 0.8427007888650501, 1e-8)
+	assert tolerance(erf(1.5), 0.966105146259005, 1e-9)
+	assert tolerance(erf(6.0), 0.99999999999999997848026328750108688340664960081261537,
+		1e-12)
+	assert tolerance(erf(5.0), 0.99999999999846254020557196514981165651461662110988195,
+		1e-12)
+	assert tolerance(erf(4.0), 0.999999984582742, 1e-12)
+	assert tolerance(erf(inf(1)), 1.0, 1e-12)
+	assert tolerance(erf(inf(-1)), -1.0, 1e-12)
 }
 
 fn test_erfc() {
-	assert almost_eq(math.erfc(-1.0), 1.84270078886505, 1e-11)
-	assert math.erfc(0.0) == 1.0
-	assert math.erfc(0.1) == 0.8875370839808279
-	assert math.erfc(0.2) == 0.7772974103342554
-}
-
-fn test_erfc_inv() {
-	assert math.erfc_inv(0.0) == math.inf(1)
-	assert math.erfc_inv(1e-100) == 15.060286697120752
-	assert math.erfc_inv(1.0) == 0.0
-	assert math.erfc_inv(0.5) == 0.47660913088024937
+	assert tolerance(erfc(-1.0), 1.84270078886505, 1e-8)
+	assert tolerance(erfc(0.0), 1.0, 1e-11)
+	assert tolerance(erfc(0.1), 0.8875370839808279, 1e-11)
+	assert tolerance(erfc(0.2), 0.7772974103342554, 1e-9)
 }
diff --git a/vlib/math/evaluate.v b/vlib/math/evaluate.v
deleted file mode 100644
index 5512208b7..000000000
--- a/vlib/math/evaluate.v
+++ /dev/null
@@ -1,18 +0,0 @@
-module math
-
-// Provides functions that don't have a numerical solution and must
-// be solved computationally (e.g. evaluation of a polynomial)
-
-pub fn polynomial(z f64, coeff []f64) f64 {
-	n := coeff.len
-	if n == 0 {
-		return 0.0
-	}
-
-	mut sum := coeff[n - 1]
-	for i := n - 1; i >= 0; i-- {
-		sum *= z
-		sum += coeff[i]
-	}
-	return sum
-}
diff --git a/vlib/math/exp.c.v b/vlib/math/exp.c.v
new file mode 100644
index 000000000..7818438a8
--- /dev/null
+++ b/vlib/math/exp.c.v
@@ -0,0 +1,17 @@
+module math
+
+fn C.exp(x f64) f64
+
+fn C.exp2(x f64) f64
+
+// exp calculates exponent of the number (math.pow(math.E, x)).
+[inline]
+pub fn exp(x f64) f64 {
+	return C.exp(x)
+}
+
+// exp2 returns the base-2 exponential function of a (math.pow(2, x)).
+[inline]
+pub fn exp2(x f64) f64 {
+	return C.exp2(x)
+}
diff --git a/vlib/math/exp.js.v b/vlib/math/exp.js.v
new file mode 100644
index 000000000..cf41f6108
--- /dev/null
+++ b/vlib/math/exp.js.v
@@ -0,0 +1,12 @@
+module math
+
+fn JS.Math.exp(x f64) f64
+
+// exp calculates exponent of the number (math.pow(math.E, x)).
+[inline]
+pub fn exp(x f64) f64 {
+	mut res := 0.0
+	#res.val = Math.exp(x)
+
+	return res
+}
diff --git a/vlib/math/exp.v b/vlib/math/exp.v
new file mode 100644
index 000000000..7a11eb6d1
--- /dev/null
+++ b/vlib/math/exp.v
@@ -0,0 +1,214 @@
+module math
+
+import math.internal
+
+const (
+	f64_max_exp = f64(1024)
+	f64_min_exp = f64(-1021)
+	threshold   = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF
+	ln2_x56     = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1
+	ln2_halfx3  = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73
+	ln2_half    = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef
+	ln2hi       = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
+	ln2lo       = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
+	inv_ln2     = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
+	// scaled coefficients related to expm1
+	expm1_q1    = -3.33333333333331316428e-02 // 0xBFA11111111110F4
+	expm1_q2    = 1.58730158725481460165e-03 // 0x3F5A01A019FE5585
+	expm1_q3    = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7
+	expm1_q4    = 4.00821782732936239552e-06 // 0x3ED0CFCA86E65239
+	expm1_q5    = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D
+)
+
+// exp returns e**x, the base-e exponential of x.
+//
+// special cases are:
+// exp(+inf) = +inf
+// exp(nan) = nan
+// Very large values overflow to 0 or +inf.
+// Very small values underflow to 1.
+pub fn exp(x f64) f64 {
+	log2e := 1.44269504088896338700e+00
+	overflow := 7.09782712893383973096e+02
+	underflow := -7.45133219101941108420e+02
+	near_zero := 1.0 / (1 << 28) // 2**-28
+	// special cases
+	if is_nan(x) || is_inf(x, 1) {
+		return x
+	}
+	if is_inf(x, -1) {
+		return 0.0
+	}
+	if x > overflow {
+		return inf(1)
+	}
+	if x < underflow {
+		return 0.0
+	}
+	if -near_zero < x && x < near_zero {
+		return 1.0 + x
+	}
+	// reduce; computed as r = hi - lo for extra precision.
+	mut k := 0
+	if x < 0 {
+		k = int(log2e * x - 0.5)
+	}
+	if x > 0 {
+		k = int(log2e * x + 0.5)
+	}
+	hi := x - f64(k) * math.ln2hi
+	lo := f64(k) * math.ln2lo
+	// compute
+	return expmulti(hi, lo, k)
+}
+
+// exp2 returns 2**x, the base-2 exponential of x.
+//
+// special cases are the same as exp.
+pub fn exp2(x f64) f64 {
+	overflow := 1.0239999999999999e+03
+	underflow := -1.0740e+03
+	if is_nan(x) || is_inf(x, 1) {
+		return x
+	}
+	if is_inf(x, -1) {
+		return 0
+	}
+	if x > overflow {
+		return inf(1)
+	}
+	if x < underflow {
+		return 0
+	}
+	// argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2.
+	// computed as r = hi - lo for extra precision.
+	mut k := 0
+	if x > 0 {
+		k = int(x + 0.5)
+	}
+	if x < 0 {
+		k = int(x - 0.5)
+	}
+	mut t := x - f64(k)
+	hi := t * math.ln2hi
+	lo := -t * math.ln2lo
+	// compute
+	return expmulti(hi, lo, k)
+}
+
+pub fn ldexp(x f64, e int) f64 {
+	if x == 0.0 {
+		return x
+	} else {
+		mut y, ex := frexp(x)
+		mut e2 := f64(e + ex)
+		if e2 >= math.f64_max_exp {
+			y *= pow(2.0, e2 - math.f64_max_exp + 1.0)
+			e2 = math.f64_max_exp - 1.0
+		} else if e2 <= math.f64_min_exp {
+			y *= pow(2.0, e2 - math.f64_min_exp - 1.0)
+			e2 = math.f64_min_exp + 1.0
+		}
+		return y * pow(2.0, e2)
+	}
+}
+
+// frexp breaks f into a normalized fraction
+// and an integral power of two.
+// It returns frac and exp satisfying f == frac × 2**exp,
+// with the absolute value of frac in the interval [½, 1).
+//
+// special cases are:
+// frexp(±0) = ±0, 0
+// frexp(±inf) = ±inf, 0
+// frexp(nan) = nan, 0
+// pub fn frexp(f f64) (f64, int) {
+// // special cases
+// if f == 0.0 {
+// return f, 0 // correctly return -0
+// }
+// if is_inf(f, 0) || is_nan(f) {
+// return f, 0
+// }
+// f_norm, mut exp := normalize(f)
+// mut x := f64_bits(f_norm)
+// exp += int((x>>shift)&mask) - bias + 1
+// x &= ~(mask << shift)
+// x |= (-1 + bias) << shift
+// return f64_from_bits(x), exp
+pub fn frexp(x f64) (f64, int) {
+	if x == 0.0 {
+		return 0.0, 0
+	} else if !is_finite(x) {
+		return x, 0
+	} else if abs(x) >= 0.5 && abs(x) < 1 { // Handle the common case
+		return x, 0
+	} else {
+		ex := ceil(log(abs(x)) / ln2)
+		mut ei := int(ex) // Prevent underflow and overflow of 2**(-ei)
+		if ei < int(math.f64_min_exp) {
+			ei = int(math.f64_min_exp)
+		}
+		if ei > -int(math.f64_min_exp) {
+			ei = -int(math.f64_min_exp)
+		}
+		mut f := x * pow(2.0, -ei)
+		if !is_finite(f) { // This should not happen
+			return f, 0
+		}
+		for abs(f) >= 1.0 {
+			ei++
+			f /= 2.0
+		}
+		for abs(f) > 0 && abs(f) < 0.5 {
+			ei--
+			f *= 2.0
+		}
+		return f, ei
+	}
+}
+
+// special cases are:
+// expm1(+inf) = +inf
+// expm1(-inf) = -1
+// expm1(nan) = nan
+pub fn expm1(x f64) f64 {
+	if is_inf(x, 1) || is_nan(x) {
+		return x
+	}
+	if is_inf(x, -1) {
+		return f64(-1)
+	}
+	// FIXME: this should be improved
+	if abs(x) < ln2 { // Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ...
+		mut i := 1.0
+		mut sum := x
+		mut term := x / 1.0
+		i++
+		term *= x / f64(i)
+		sum += term
+		for abs(term) > abs(sum) * internal.f64_epsilon {
+			i++
+			term *= x / f64(i)
+			sum += term
+		}
+		return sum
+	} else {
+		return exp(x) - 1
+	}
+}
+
+// exp1 returns e**r × 2**k where r = hi - lo and |r| ≤ ln(2)/2.
+fn expmulti(hi f64, lo f64, k int) f64 {
+	exp_p1 := 1.66666666666666657415e-01 // 0x3FC55555; 0x55555555
+	exp_p2 := -2.77777777770155933842e-03 // 0xBF66C16C; 0x16BEBD93
+	exp_p3 := 6.61375632143793436117e-05 // 0x3F11566A; 0xAF25DE2C
+	exp_p4 := -1.65339022054652515390e-06 // 0xBEBBBD41; 0xC5D26BF1
+	exp_p5 := 4.13813679705723846039e-08 // 0x3E663769; 0x72BEA4D0
+	r := hi - lo
+	t := r * r
+	c := r - t * (exp_p1 + t * (exp_p2 + t * (exp_p3 + t * (exp_p4 + t * exp_p5))))
+	y := 1 - ((lo - (r * c) / (2 - c)) - hi)
+	// TODO(rsc): make sure ldexp can handle boundary k
+	return ldexp(y, k)
+}
diff --git a/vlib/math/factorial/factorial.v b/vlib/math/factorial.v
similarity index 60%
rename from vlib/math/factorial/factorial.v
rename to vlib/math/factorial.v
index 9668d5d4c..116e08398 100644
--- a/vlib/math/factorial/factorial.v
+++ b/vlib/math/factorial.v
@@ -1,48 +1,31 @@
-// Copyright (c) 2019-2021 Alexander Medvednikov. All rights reserved.
-// Use of this source code is governed by an MIT license
-// that can be found in the LICENSE file.
-
-// Module created by Ulises Jeremias Cornejo Fandos based on
-// the definitions provided in https://scientificc.github.io/cmathl/
-
-module factorial
-
-import math
+module math
 
 // factorial calculates the factorial of the provided value.
 pub fn factorial(n f64) f64 {
-	// For a large postive argument (n >= FACTORIALS.len) return max_f64
-
+	// For a large postive argument (n >= factorials_table.len) return max_f64
 	if n >= factorials_table.len {
-		return math.max_f64
+		return max_f64
 	}
-
 	// Otherwise return n!.
 	if n == f64(i64(n)) && n >= 0.0 {
 		return factorials_table[i64(n)]
 	}
-
-	return math.gamma(n + 1.0)
+	return gamma(n + 1.0)
 }
 
 // log_factorial calculates the log-factorial of the provided value.
 pub fn log_factorial(n f64) f64 {
 	// For a large postive argument (n < 0) return max_f64
-
 	if n < 0 {
-		return -math.max_f64
+		return -max_f64
 	}
-
 	// If n < N then return ln(n!).
-
 	if n != f64(i64(n)) {
-		return math.log_gamma(n + 1)
+		return log_gamma(n + 1)
 	} else if n < log_factorials_table.len {
 		return log_factorials_table[i64(n)]
 	}
-
 	// Otherwise return asymptotic expansion of ln(n!).
-
 	return log_factorial_asymptotic_expansion(int(n))
 }
 
@@ -51,30 +34,22 @@ fn log_factorial_asymptotic_expansion(n int) f64 {
 	mut term := []f64{}
 	xx := f64((n + 1) * (n + 1))
 	mut xj := f64(n + 1)
-
-	log_factorial := log_sqrt_2pi - xj + (xj - 0.5) * math.log(xj)
-
+	log_factorial := log_sqrt_2pi - xj + (xj - 0.5) * log(xj)
 	mut i := 0
-
 	for i = 0; i < m; i++ {
-		term << b_numbers[i] / xj
+		term << bernoulli[i] / xj
 		xj *= xx
 	}
-
 	mut sum := term[m - 1]
-
 	for i = m - 2; i >= 0; i-- {
-		if math.abs(sum) <= math.abs(term[i]) {
+		if abs(sum) <= abs(term[i]) {
 			break
 		}
-
 		sum = term[i]
 	}
-
 	for i >= 0 {
 		sum += term[i]
 		i--
 	}
-
 	return log_factorial + sum
 }
diff --git a/vlib/math/factorial/factorial_tables.v b/vlib/math/factorial/factorial_tables.v
deleted file mode 100644
index a669b0919..000000000
--- a/vlib/math/factorial/factorial_tables.v
+++ /dev/null
@@ -1,375 +0,0 @@
-// Copyright (c) 2019-2021 Alexander Medvednikov. All rights reserved.
-// Use of this source code is governed by an MIT license
-// that can be found in the LICENSE file.
-
-module factorial
-
-const (
-	log_sqrt_2pi         = 9.18938533204672741780329736e-1
-
-	b_numbers            = [
-		/*
-		Bernoulli numbers B(2),B(4),B(6),...,B(20).  Only B(2),...,B(10) currently
-        * used.
-		*/
-		f64(1.0 / (6.0 * 2.0 * 1.0)),
-		-1.0 / (30.0 * 4.0 * 3.0),
-		1.0 / (42.0 * 6.0 * 5.0),
-		-1.0 / (30.0 * 8.0 * 7.0),
-		5.0 / (66.0 * 10.0 * 9.0),
-		-691.0 / (2730.0 * 12.0 * 11.0),
-		7.0 / (6.0 * 14.0 * 13.0),
-		-3617.0 / (510.0 * 16.0 * 15.0),
-		43867.0 / (796.0 * 18.0 * 17.0),
-		-174611.0 / (330.0 * 20.0 * 19.0),
-	]
-
-	factorials_table     = [
-		f64(1.000000000000000000000e+0), /* 0! */
-		1.000000000000000000000e+0, /* 1! */
-		2.000000000000000000000e+0, /* 2! */
-		6.000000000000000000000e+0, /* 3! */
-		2.400000000000000000000e+1, /* 4! */
-		1.200000000000000000000e+2, /* 5! */
-		7.200000000000000000000e+2, /* 6! */
-		5.040000000000000000000e+3, /* 7! */
-		4.032000000000000000000e+4, /* 8! */
-		3.628800000000000000000e+5, /* 9! */
-		3.628800000000000000000e+6, /* 10! */
-		3.991680000000000000000e+7, /* 11! */
-		4.790016000000000000000e+8, /* 12! */
-		6.227020800000000000000e+9, /* 13! */
-		8.717829120000000000000e+10, /* 14! */
-		1.307674368000000000000e+12, /* 15! */
-		2.092278988800000000000e+13, /* 16! */
-		3.556874280960000000000e+14, /* 17! */
-		6.402373705728000000000e+15, /* 18! */
-		1.216451004088320000000e+17, /* 19! */
-		2.432902008176640000000e+18, /* 20! */
-		5.109094217170944000000e+19, /* 21! */
-		1.124000727777607680000e+21, /* 22! */
-		2.585201673888497664000e+22, /* 23! */
-		6.204484017332394393600e+23, /* 24! */
-		1.551121004333098598400e+25, /* 25! */
-		4.032914611266056355840e+26, /* 26! */
-		1.088886945041835216077e+28, /* 27! */
-		3.048883446117138605015e+29, /* 28! */
-		8.841761993739701954544e+30, /* 29! */
-		2.652528598121910586363e+32, /* 30! */
-		8.222838654177922817726e+33, /* 31! */
-		2.631308369336935301672e+35, /* 32! */
-		8.683317618811886495518e+36, /* 33! */
-		2.952327990396041408476e+38, /* 34! */
-		1.033314796638614492967e+40, /* 35! */
-		3.719933267899012174680e+41, /* 36! */
-		1.376375309122634504632e+43, /* 37! */
-		5.230226174666011117600e+44, /* 38! */
-		2.039788208119744335864e+46, /* 39! */
-		8.159152832478977343456e+47, /* 40! */
-		3.345252661316380710817e+49, /* 41! */
-		1.405006117752879898543e+51, /* 42! */
-		6.041526306337383563736e+52, /* 43! */
-		2.658271574788448768044e+54, /* 44! */
-		1.196222208654801945620e+56, /* 45! */
-		5.502622159812088949850e+57, /* 46! */
-		2.586232415111681806430e+59, /* 47! */
-		1.241391559253607267086e+61, /* 48! */
-		6.082818640342675608723e+62, /* 49! */
-		3.041409320171337804361e+64, /* 50! */
-		1.551118753287382280224e+66, /* 51! */
-		8.065817517094387857166e+67, /* 52! */
-		4.274883284060025564298e+69, /* 53! */
-		2.308436973392413804721e+71, /* 54! */
-		1.269640335365827592597e+73, /* 55! */
-		7.109985878048634518540e+74, /* 56! */
-		4.052691950487721675568e+76, /* 57! */
-		2.350561331282878571829e+78, /* 58! */
-		1.386831185456898357379e+80, /* 59! */
-		8.320987112741390144276e+81, /* 60! */
-		5.075802138772247988009e+83, /* 61! */
-		3.146997326038793752565e+85, /* 62! */
-		1.982608315404440064116e+87, /* 63! */
-		1.268869321858841641034e+89, /* 64! */
-		8.247650592082470666723e+90, /* 65! */
-		5.443449390774430640037e+92, /* 66! */
-		3.647111091818868528825e+94, /* 67! */
-		2.480035542436830599601e+96, /* 68! */
-		1.711224524281413113725e+98, /* 69! */
-		1.197857166996989179607e+100, /* 70! */
-		8.504785885678623175212e+101, /* 71! */
-		6.123445837688608686152e+103, /* 72! */
-		4.470115461512684340891e+105, /* 73! */
-		3.307885441519386412260e+107, /* 74! */
-		2.480914081139539809195e+109, /* 75! */
-		1.885494701666050254988e+111, /* 76! */
-		1.451830920282858696341e+113, /* 77! */
-		1.132428117820629783146e+115, /* 78! */
-		8.946182130782975286851e+116, /* 79! */
-		7.156945704626380229481e+118, /* 80! */
-		5.797126020747367985880e+120, /* 81! */
-		4.753643337012841748421e+122, /* 82! */
-		3.945523969720658651190e+124, /* 83! */
-		3.314240134565353266999e+126, /* 84! */
-		2.817104114380550276949e+128, /* 85! */
-		2.422709538367273238177e+130, /* 86! */
-		2.107757298379527717214e+132, /* 87! */
-		1.854826422573984391148e+134, /* 88! */
-		1.650795516090846108122e+136, /* 89! */
-		1.485715964481761497310e+138, /* 90! */
-		1.352001527678402962552e+140, /* 91! */
-		1.243841405464130725548e+142, /* 92! */
-		1.156772507081641574759e+144, /* 93! */
-		1.087366156656743080274e+146, /* 94! */
-		1.032997848823905926260e+148, /* 95! */
-		9.916779348709496892096e+149, /* 96! */
-		9.619275968248211985333e+151, /* 97! */
-		9.426890448883247745626e+153, /* 98! */
-		9.332621544394415268170e+155, /* 99! */
-		9.332621544394415268170e+157, /* 100! */
-		9.425947759838359420852e+159, /* 101! */
-		9.614466715035126609269e+161, /* 102! */
-		9.902900716486180407547e+163, /* 103! */
-		1.029901674514562762385e+166, /* 104! */
-		1.081396758240290900504e+168, /* 105! */
-		1.146280563734708354534e+170, /* 106! */
-		1.226520203196137939352e+172, /* 107! */
-		1.324641819451828974500e+174, /* 108! */
-		1.443859583202493582205e+176, /* 109! */
-		1.588245541522742940425e+178, /* 110! */
-		1.762952551090244663872e+180, /* 111! */
-		1.974506857221074023537e+182, /* 112! */
-		2.231192748659813646597e+184, /* 113! */
-		2.543559733472187557120e+186, /* 114! */
-		2.925093693493015690688e+188, /* 115! */
-		3.393108684451898201198e+190, /* 116! */
-		3.969937160808720895402e+192, /* 117! */
-		4.684525849754290656574e+194, /* 118! */
-		5.574585761207605881323e+196, /* 119! */
-		6.689502913449127057588e+198, /* 120! */
-		8.094298525273443739682e+200, /* 121! */
-		9.875044200833601362412e+202, /* 122! */
-		1.214630436702532967577e+205, /* 123! */
-		1.506141741511140879795e+207, /* 124! */
-		1.882677176888926099744e+209, /* 125! */
-		2.372173242880046885677e+211, /* 126! */
-		3.012660018457659544810e+213, /* 127! */
-		3.856204823625804217357e+215, /* 128! */
-		4.974504222477287440390e+217, /* 129! */
-		6.466855489220473672507e+219, /* 130! */
-		8.471580690878820510985e+221, /* 131! */
-		1.118248651196004307450e+224, /* 132! */
-		1.487270706090685728908e+226, /* 133! */
-		1.992942746161518876737e+228, /* 134! */
-		2.690472707318050483595e+230, /* 135! */
-		3.659042881952548657690e+232, /* 136! */
-		5.012888748274991661035e+234, /* 137! */
-		6.917786472619488492228e+236, /* 138! */
-		9.615723196941089004197e+238, /* 139! */
-		1.346201247571752460588e+241, /* 140! */
-		1.898143759076170969429e+243, /* 141! */
-		2.695364137888162776589e+245, /* 142! */
-		3.854370717180072770522e+247, /* 143! */
-		5.550293832739304789551e+249, /* 144! */
-		8.047926057471991944849e+251, /* 145! */
-		1.174997204390910823948e+254, /* 146! */
-		1.727245890454638911203e+256, /* 147! */
-		2.556323917872865588581e+258, /* 148! */
-		3.808922637630569726986e+260, /* 149! */
-		5.713383956445854590479e+262, /* 150! */
-		8.627209774233240431623e+264, /* 151! */
-		1.311335885683452545607e+267, /* 152! */
-		2.006343905095682394778e+269, /* 153! */
-		3.089769613847350887959e+271, /* 154! */
-		4.789142901463393876336e+273, /* 155! */
-		7.471062926282894447084e+275, /* 156! */
-		1.172956879426414428192e+278, /* 157! */
-		1.853271869493734796544e+280, /* 158! */
-		2.946702272495038326504e+282, /* 159! */
-		4.714723635992061322407e+284, /* 160! */
-		7.590705053947218729075e+286, /* 161! */
-		1.229694218739449434110e+289, /* 162! */
-		2.004401576545302577600e+291, /* 163! */
-		3.287218585534296227263e+293, /* 164! */
-		5.423910666131588774984e+295, /* 165! */
-		9.003691705778437366474e+297, /* 166! */
-		1.503616514864999040201e+300, /* 167! */
-		2.526075744973198387538e+302, /* 168! */
-		4.269068009004705274939e+304, /* 169! */
-		7.257415615307998967397e+306, /* 170! */
-	]
-
-	log_factorials_table = [
-		f64(0.000000000000000000000e+0), /* 0! */
-		0.000000000000000000000e+0, /* 1! */
-		6.931471805599453094172e-1, /* 2! */
-		1.791759469228055000812e+0, /* 3! */
-		3.178053830347945619647e+0, /* 4! */
-		4.787491742782045994248e+0, /* 5! */
-		6.579251212010100995060e+0, /* 6! */
-		8.525161361065414300166e+0, /* 7! */
-		1.060460290274525022842e+1, /* 8! */
-		1.280182748008146961121e+1, /* 9! */
-		1.510441257307551529523e+1, /* 10! */
-		1.750230784587388583929e+1, /* 11! */
-		1.998721449566188614952e+1, /* 12! */
-		2.255216385312342288557e+1, /* 13! */
-		2.519122118273868150009e+1, /* 14! */
-		2.789927138384089156609e+1, /* 15! */
-		3.067186010608067280376e+1, /* 16! */
-		3.350507345013688888401e+1, /* 17! */
-		3.639544520803305357622e+1, /* 18! */
-		3.933988418719949403622e+1, /* 19! */
-		4.233561646075348502966e+1, /* 20! */
-		4.538013889847690802616e+1, /* 21! */
-		4.847118135183522387964e+1, /* 22! */
-		5.160667556776437357045e+1, /* 23! */
-		5.478472939811231919009e+1, /* 24! */
-		5.800360522298051993929e+1, /* 25! */
-		6.126170176100200198477e+1, /* 26! */
-		6.455753862700633105895e+1, /* 27! */
-		6.788974313718153498289e+1, /* 28! */
-		7.125703896716800901007e+1, /* 29! */
-		7.465823634883016438549e+1, /* 30! */
-		7.809222355331531063142e+1, /* 31! */
-		8.155795945611503717850e+1, /* 32! */
-		8.505446701758151741396e+1, /* 33! */
-		8.858082754219767880363e+1, /* 34! */
-		9.213617560368709248333e+1, /* 35! */
-		9.571969454214320248496e+1, /* 36! */
-		9.933061245478742692933e+1, /* 37! */
-		1.029681986145138126988e+2, /* 38! */
-		1.066317602606434591262e+2, /* 39! */
-		1.103206397147573954291e+2, /* 40! */
-		1.140342117814617032329e+2, /* 41! */
-		1.177718813997450715388e+2, /* 42! */
-		1.215330815154386339623e+2, /* 43! */
-		1.253172711493568951252e+2, /* 44! */
-		1.291239336391272148826e+2, /* 45! */
-		1.329525750356163098828e+2, /* 46! */
-		1.368027226373263684696e+2, /* 47! */
-		1.406739236482342593987e+2, /* 48! */
-		1.445657439463448860089e+2, /* 49! */
-		1.484777669517730320675e+2, /* 50! */
-		1.524095925844973578392e+2, /* 51! */
-		1.563608363030787851941e+2, /* 52! */
-		1.603311282166309070282e+2, /* 53! */
-		1.643201122631951814118e+2, /* 54! */
-		1.683274454484276523305e+2, /* 55! */
-		1.723527971391628015638e+2, /* 56! */
-		1.763958484069973517152e+2, /* 57! */
-		1.804562914175437710518e+2, /* 58! */
-		1.845338288614494905025e+2, /* 59! */
-		1.886281734236715911873e+2, /* 60! */
-		1.927390472878449024360e+2, /* 61! */
-		1.968661816728899939914e+2, /* 62! */
-		2.010093163992815266793e+2, /* 63! */
-		2.051681994826411985358e+2, /* 64! */
-		2.093425867525368356464e+2, /* 65! */
-		2.135322414945632611913e+2, /* 66! */
-		2.177369341139542272510e+2, /* 67! */
-		2.219564418191303339501e+2, /* 68! */
-		2.261905483237275933323e+2, /* 69! */
-		2.304390435657769523214e+2, /* 70! */
-		2.347017234428182677427e+2, /* 71! */
-		2.389783895618343230538e+2, /* 72! */
-		2.432688490029827141829e+2, /* 73! */
-		2.475729140961868839366e+2, /* 74! */
-		2.518904022097231943772e+2, /* 75! */
-		2.562211355500095254561e+2, /* 76! */
-		2.605649409718632093053e+2, /* 77! */
-		2.649216497985528010421e+2, /* 78! */
-		2.692910976510198225363e+2, /* 79! */
-		2.736731242856937041486e+2, /* 80! */
-		2.780675734403661429141e+2, /* 81! */
-		2.824742926876303960274e+2, /* 82! */
-		2.868931332954269939509e+2, /* 83! */
-		2.913239500942703075662e+2, /* 84! */
-		2.957666013507606240211e+2, /* 85! */
-		3.002209486470141317540e+2, /* 86! */
-		3.046868567656687154726e+2, /* 87! */
-		3.091641935801469219449e+2, /* 88! */
-		3.136528299498790617832e+2, /* 89! */
-		3.181526396202093268500e+2, /* 90! */
-		3.226634991267261768912e+2, /* 91! */
-		3.271852877037752172008e+2, /* 92! */
-		3.317178871969284731381e+2, /* 93! */
-		3.362611819791984770344e+2, /* 94! */
-		3.408150588707990178690e+2, /* 95! */
-		3.453794070622668541074e+2, /* 96! */
-		3.499541180407702369296e+2, /* 97! */
-		3.545390855194408088492e+2, /* 98! */
-		3.591342053695753987760e+2, /* 99! */
-		3.637393755555634901441e+2, /* 100! */
-		3.683544960724047495950e+2, /* 101! */
-		3.729794688856890206760e+2, /* 102! */
-		3.776141978739186564468e+2, /* 103! */
-		3.822585887730600291111e+2, /* 104! */
-		3.869125491232175524822e+2, /* 105! */
-		3.915759882173296196258e+2, /* 106! */
-		3.962488170517915257991e+2, /* 107! */
-		4.009309482789157454921e+2, /* 108! */
-		4.056222961611448891925e+2, /* 109! */
-		4.103227765269373054205e+2, /* 110! */
-		4.150323067282496395563e+2, /* 111! */
-		4.197508055995447340991e+2, /* 112! */
-		4.244781934182570746677e+2, /* 113! */
-		4.292143918666515701285e+2, /* 114! */
-		4.339593239950148201939e+2, /* 115! */
-		4.387129141861211848399e+2, /* 116! */
-		4.434750881209189409588e+2, /* 117! */
-		4.482457727453846057188e+2, /* 118! */
-		4.530248962384961351041e+2, /* 119! */
-		4.578123879812781810984e+2, /* 120! */
-		4.626081785268749221865e+2, /* 121! */
-		4.674121995716081787447e+2, /* 122! */
-		4.722243839269805962399e+2, /* 123! */
-		4.770446654925856331047e+2, /* 124! */
-		4.818729792298879342285e+2, /* 125! */
-		4.867092611368394122258e+2, /* 126! */
-		4.915534482232980034989e+2, /* 127! */
-		4.964054784872176206648e+2, /* 128! */
-		5.012652908915792927797e+2, /* 129! */
-		5.061328253420348751997e+2, /* 130! */
-		5.110080226652360267439e+2, /* 131! */
-		5.158908245878223975982e+2, /* 132! */
-		5.207811737160441513633e+2, /* 133! */
-		5.256790135159950627324e+2, /* 134! */
-		5.305842882944334921812e+2, /* 135! */
-		5.354969431801695441897e+2, /* 136! */
-		5.404169241059976691050e+2, /* 137! */
-		5.453441777911548737966e+2, /* 138! */
-		5.502786517242855655538e+2, /* 139! */
-		5.552202941468948698523e+2, /* 140! */
-		5.601690540372730381305e+2, /* 141! */
-		5.651248810948742988613e+2, /* 142! */
-		5.700877257251342061414e+2, /* 143! */
-		5.750575390247102067619e+2, /* 144! */
-		5.800342727671307811636e+2, /* 145! */
-		5.850178793888391176022e+2, /* 146! */
-		5.900083119756178539038e+2, /* 147! */
-		5.950055242493819689670e+2, /* 148! */
-		6.000094705553274281080e+2, /* 149! */
-		6.050201058494236838580e+2, /* 150! */
-		6.100373856862386081868e+2, /* 151! */
-		6.150612662070848845750e+2, /* 152! */
-		6.200917041284773200381e+2, /* 153! */
-		6.251286567308909491967e+2, /* 154! */
-		6.301720818478101958172e+2, /* 155! */
-		6.352219378550597328635e+2, /* 156! */
-		6.402781836604080409209e+2, /* 157! */
-		6.453407786934350077245e+2, /* 158! */
-		6.504096828956552392500e+2, /* 159! */
-		6.554848567108890661717e+2, /* 160! */
-		6.605662610758735291676e+2, /* 161! */
-		6.656538574111059132426e+2, /* 162! */
-		6.707476076119126755767e+2, /* 163! */
-		6.758474740397368739994e+2, /* 164! */
-		6.809534195136374546094e+2, /* 165! */
-		6.860654073019939978423e+2, /* 166! */
-		6.911834011144107529496e+2, /* 167! */
-		6.963073650938140118743e+2, /* 168! */
-		7.014372638087370853465e+2, /* 169! */
-		7.065730622457873471107e+2, /* 170! */
-		7.117147258022900069535e+2, /* 171! */
-	]
-)
diff --git a/vlib/math/factorial/factorial_test.v b/vlib/math/factorial/factorial_test.v
deleted file mode 100644
index 6c2b57575..000000000
--- a/vlib/math/factorial/factorial_test.v
+++ /dev/null
@@ -1,14 +0,0 @@
-import math
-import math.factorial as fact
-
-fn test_factorial() {
-	assert fact.factorial(12) == 479001600
-	assert fact.factorial(5) == 120
-	assert fact.factorial(0) == 1
-}
-
-fn test_log_factorial() {
-	assert fact.log_factorial(12) == math.log(479001600)
-	assert fact.log_factorial(5) == math.log(120)
-	assert fact.log_factorial(0) == math.log(1)
-}
diff --git a/vlib/math/factorial_tables.v b/vlib/math/factorial_tables.v
new file mode 100644
index 000000000..5154b218c
--- /dev/null
+++ b/vlib/math/factorial_tables.v
@@ -0,0 +1,711 @@
+module math
+
+const (
+	log_sqrt_2pi         = 9.18938533204672741780329736e-1
+	bernoulli            = [
+		/*
+		Bernoulli numbers B(2),B(4),B(6),...,B(20).  Only B(2),...,B(10) currently
+        * used.
+		*/
+		1.0 / (6.0 * 2.0 * 1.0),
+		-1.0 / (30.0 * 4.0 * 3.0),
+		1.0 / (42.0 * 6.0 * 5.0),
+		-1.0 / (30.0 * 8.0 * 7.0),
+		5.0 / (66.0 * 10.0 * 9.0),
+		-691.0 / (2730.0 * 12.0 * 11.0),
+		7.0 / (6.0 * 14.0 * 13.0),
+		-3617.0 / (510.0 * 16.0 * 15.0),
+		43867.0 / (796.0 * 18.0 * 17.0),
+		-174611.0 / (330.0 * 20.0 * 19.0),
+	]
+	factorials_table     = [
+		// 0!
+		1.000000000000000000000e+0,
+		// 1!
+		1.000000000000000000000e+0,
+		// 2!
+		2.000000000000000000000e+0,
+		// 3!
+		6.000000000000000000000e+0,
+		// 4!
+		2.400000000000000000000e+1,
+		// 5!
+		1.200000000000000000000e+2,
+		// 6!
+		7.200000000000000000000e+2,
+		// 7!
+		5.040000000000000000000e+3,
+		// 8!
+		4.032000000000000000000e+4,
+		// 9!
+		3.628800000000000000000e+5,
+		// 10!
+		3.628800000000000000000e+6,
+		// 11!
+		3.991680000000000000000e+7,
+		// 12!
+		4.790016000000000000000e+8,
+		// 13!
+		6.227020800000000000000e+9,
+		// 14!
+		8.717829120000000000000e+10,
+		// 15!
+		1.307674368000000000000e+12,
+		// 16!
+		2.092278988800000000000e+13,
+		// 17!
+		3.556874280960000000000e+14,
+		// 18!
+		6.402373705728000000000e+15,
+		// 19!
+		1.216451004088320000000e+17,
+		// 20!
+		2.432902008176640000000e+18,
+		// 21!
+		5.109094217170944000000e+19,
+		// 22!
+		1.124000727777607680000e+21,
+		// 23!
+		2.585201673888497664000e+22,
+		// 24!
+		6.204484017332394393600e+23,
+		// 25!
+		1.551121004333098598400e+25,
+		// 26!
+		4.032914611266056355840e+26,
+		// 27!
+		1.088886945041835216077e+28,
+		// 28!
+		3.048883446117138605015e+29,
+		// 29!
+		8.841761993739701954544e+30,
+		// 30!
+		2.652528598121910586363e+32,
+		// 31!
+		8.222838654177922817726e+33,
+		// 32!
+		2.631308369336935301672e+35,
+		// 33!
+		8.683317618811886495518e+36,
+		// 34!
+		2.952327990396041408476e+38,
+		// 35!
+		1.033314796638614492967e+40,
+		// 36!
+		3.719933267899012174680e+41,
+		// 37!
+		1.376375309122634504632e+43,
+		// 38!
+		5.230226174666011117600e+44,
+		// 39!
+		2.039788208119744335864e+46,
+		// 40!
+		8.159152832478977343456e+47,
+		// 41!
+		3.345252661316380710817e+49,
+		// 42!
+		1.405006117752879898543e+51,
+		// 43!
+		6.041526306337383563736e+52,
+		// 44!
+		2.658271574788448768044e+54,
+		// 45!
+		1.196222208654801945620e+56,
+		// 46!
+		5.502622159812088949850e+57,
+		// 47!
+		2.586232415111681806430e+59,
+		// 48!
+		1.241391559253607267086e+61,
+		// 49!
+		6.082818640342675608723e+62,
+		// 50!
+		3.041409320171337804361e+64,
+		// 51!
+		1.551118753287382280224e+66,
+		// 52!
+		8.065817517094387857166e+67,
+		// 53!
+		4.274883284060025564298e+69,
+		// 54!
+		2.308436973392413804721e+71,
+		// 55!
+		1.269640335365827592597e+73,
+		// 56!
+		7.109985878048634518540e+74,
+		// 57!
+		4.052691950487721675568e+76,
+		// 58!
+		2.350561331282878571829e+78,
+		// 59!
+		1.386831185456898357379e+80,
+		// 60!
+		8.320987112741390144276e+81,
+		// 61!
+		5.075802138772247988009e+83,
+		// 62!
+		3.146997326038793752565e+85,
+		// 63!
+		1.982608315404440064116e+87,
+		// 64!
+		1.268869321858841641034e+89,
+		// 65!
+		8.247650592082470666723e+90,
+		// 66!
+		5.443449390774430640037e+92,
+		// 67!
+		3.647111091818868528825e+94,
+		// 68!
+		2.480035542436830599601e+96,
+		// 69!
+		1.711224524281413113725e+98,
+		// 70!
+		1.197857166996989179607e+100,
+		// 71!
+		8.504785885678623175212e+101,
+		// 72!
+		6.123445837688608686152e+103,
+		// 73!
+		4.470115461512684340891e+105,
+		// 74!
+		3.307885441519386412260e+107,
+		// 75!
+		2.480914081139539809195e+109,
+		// 76!
+		1.885494701666050254988e+111,
+		// 77!
+		1.451830920282858696341e+113,
+		// 78!
+		1.132428117820629783146e+115,
+		// 79!
+		8.946182130782975286851e+116,
+		// 80!
+		7.156945704626380229481e+118,
+		// 81!
+		5.797126020747367985880e+120,
+		// 82!
+		4.753643337012841748421e+122,
+		// 83!
+		3.945523969720658651190e+124,
+		// 84!
+		3.314240134565353266999e+126,
+		// 85!
+		2.817104114380550276949e+128,
+		// 86!
+		2.422709538367273238177e+130,
+		// 87!
+		2.107757298379527717214e+132,
+		// 88!
+		1.854826422573984391148e+134,
+		// 89!
+		1.650795516090846108122e+136,
+		// 90!
+		1.485715964481761497310e+138,
+		// 91!
+		1.352001527678402962552e+140,
+		// 92!
+		1.243841405464130725548e+142,
+		// 93!
+		1.156772507081641574759e+144,
+		// 94!
+		1.087366156656743080274e+146,
+		// 95!
+		1.032997848823905926260e+148,
+		// 96!
+		9.916779348709496892096e+149,
+		// 97!
+		9.619275968248211985333e+151,
+		// 98!
+		9.426890448883247745626e+153,
+		// 99!
+		9.332621544394415268170e+155,
+		// 100!
+		9.332621544394415268170e+157,
+		// 101!
+		9.425947759838359420852e+159,
+		// 102!
+		9.614466715035126609269e+161,
+		// 103!
+		9.902900716486180407547e+163,
+		// 104!
+		1.029901674514562762385e+166,
+		// 105!
+		1.081396758240290900504e+168,
+		// 106!
+		1.146280563734708354534e+170,
+		// 107!
+		1.226520203196137939352e+172,
+		// 108!
+		1.324641819451828974500e+174,
+		// 109!
+		1.443859583202493582205e+176,
+		// 110!
+		1.588245541522742940425e+178,
+		// 111!
+		1.762952551090244663872e+180,
+		// 112!
+		1.974506857221074023537e+182,
+		// 113!
+		2.231192748659813646597e+184,
+		// 114!
+		2.543559733472187557120e+186,
+		// 115!
+		2.925093693493015690688e+188,
+		// 116!
+		3.393108684451898201198e+190,
+		// 117!
+		3.969937160808720895402e+192,
+		// 118!
+		4.684525849754290656574e+194,
+		// 119!
+		5.574585761207605881323e+196,
+		// 120!
+		6.689502913449127057588e+198,
+		// 121!
+		8.094298525273443739682e+200,
+		// 122!
+		9.875044200833601362412e+202,
+		// 123!
+		1.214630436702532967577e+205,
+		// 124!
+		1.506141741511140879795e+207,
+		// 125!
+		1.882677176888926099744e+209,
+		// 126!
+		2.372173242880046885677e+211,
+		// 127!
+		3.012660018457659544810e+213,
+		// 128!
+		3.856204823625804217357e+215,
+		// 129!
+		4.974504222477287440390e+217,
+		// 130!
+		6.466855489220473672507e+219,
+		// 131!
+		8.471580690878820510985e+221,
+		// 132!
+		1.118248651196004307450e+224,
+		// 133!
+		1.487270706090685728908e+226,
+		// 134!
+		1.992942746161518876737e+228,
+		// 135!
+		2.690472707318050483595e+230,
+		// 136!
+		3.659042881952548657690e+232,
+		// 137!
+		5.012888748274991661035e+234,
+		// 138!
+		6.917786472619488492228e+236,
+		// 139!
+		9.615723196941089004197e+238,
+		// 140!
+		1.346201247571752460588e+241,
+		// 141!
+		1.898143759076170969429e+243,
+		// 142!
+		2.695364137888162776589e+245,
+		// 143!
+		3.854370717180072770522e+247,
+		// 144!
+		5.550293832739304789551e+249,
+		// 145!
+		8.047926057471991944849e+251,
+		// 146!
+		1.174997204390910823948e+254,
+		// 147!
+		1.727245890454638911203e+256,
+		// 148!
+		2.556323917872865588581e+258,
+		// 149!
+		3.808922637630569726986e+260,
+		// 150!
+		5.713383956445854590479e+262,
+		// 151!
+		8.627209774233240431623e+264,
+		// 152!
+		1.311335885683452545607e+267,
+		// 153!
+		2.006343905095682394778e+269,
+		// 154!
+		3.089769613847350887959e+271,
+		// 155!
+		4.789142901463393876336e+273,
+		// 156!
+		7.471062926282894447084e+275,
+		// 157!
+		1.172956879426414428192e+278,
+		// 158!
+		1.853271869493734796544e+280,
+		// 159!
+		2.946702272495038326504e+282,
+		// 160!
+		4.714723635992061322407e+284,
+		// 161!
+		7.590705053947218729075e+286,
+		// 162!
+		1.229694218739449434110e+289,
+		// 163!
+		2.004401576545302577600e+291,
+		// 164!
+		3.287218585534296227263e+293,
+		// 165!
+		5.423910666131588774984e+295,
+		// 166!
+		9.003691705778437366474e+297,
+		// 167!
+		1.503616514864999040201e+300,
+		// 168!
+		2.526075744973198387538e+302,
+		// 169!
+		4.269068009004705274939e+304,
+		// 170!
+		7.257415615307998967397e+306,
+	]
+	log_factorials_table = [
+		// 0!
+		0.000000000000000000000e+0,
+		// 1!
+		0.000000000000000000000e+0,
+		// 2!
+		6.931471805599453094172e-1,
+		// 3!
+		1.791759469228055000812e+0,
+		// 4!
+		3.178053830347945619647e+0,
+		// 5!
+		4.787491742782045994248e+0,
+		// 6!
+		6.579251212010100995060e+0,
+		// 7!
+		8.525161361065414300166e+0,
+		// 8!
+		1.060460290274525022842e+1,
+		// 9!
+		1.280182748008146961121e+1,
+		// 10!
+		1.510441257307551529523e+1,
+		// 11!
+		1.750230784587388583929e+1,
+		// 12!
+		1.998721449566188614952e+1,
+		// 13!
+		2.255216385312342288557e+1,
+		// 14!
+		2.519122118273868150009e+1,
+		// 15!
+		2.789927138384089156609e+1,
+		// 16!
+		3.067186010608067280376e+1,
+		// 17!
+		3.350507345013688888401e+1,
+		// 18!
+		3.639544520803305357622e+1,
+		// 19!
+		3.933988418719949403622e+1,
+		// 20!
+		4.233561646075348502966e+1,
+		// 21!
+		4.538013889847690802616e+1,
+		// 22!
+		4.847118135183522387964e+1,
+		// 23!
+		5.160667556776437357045e+1,
+		// 24!
+		5.478472939811231919009e+1,
+		// 25!
+		5.800360522298051993929e+1,
+		// 26!
+		6.126170176100200198477e+1,
+		// 27!
+		6.455753862700633105895e+1,
+		// 28!
+		6.788974313718153498289e+1,
+		// 29!
+		7.125703896716800901007e+1,
+		// 30!
+		7.465823634883016438549e+1,
+		// 31!
+		7.809222355331531063142e+1,
+		// 32!
+		8.155795945611503717850e+1,
+		// 33!
+		8.505446701758151741396e+1,
+		// 34!
+		8.858082754219767880363e+1,
+		// 35!
+		9.213617560368709248333e+1,
+		// 36!
+		9.571969454214320248496e+1,
+		// 37!
+		9.933061245478742692933e+1,
+		// 38!
+		1.029681986145138126988e+2,
+		// 39!
+		1.066317602606434591262e+2,
+		// 40!
+		1.103206397147573954291e+2,
+		// 41!
+		1.140342117814617032329e+2,
+		// 42!
+		1.177718813997450715388e+2,
+		// 43!
+		1.215330815154386339623e+2,
+		// 44!
+		1.253172711493568951252e+2,
+		// 45!
+		1.291239336391272148826e+2,
+		// 46!
+		1.329525750356163098828e+2,
+		// 47!
+		1.368027226373263684696e+2,
+		// 48!
+		1.406739236482342593987e+2,
+		// 49!
+		1.445657439463448860089e+2,
+		// 50!
+		1.484777669517730320675e+2,
+		// 51!
+		1.524095925844973578392e+2,
+		// 52!
+		1.563608363030787851941e+2,
+		// 53!
+		1.603311282166309070282e+2,
+		// 54!
+		1.643201122631951814118e+2,
+		// 55!
+		1.683274454484276523305e+2,
+		// 56!
+		1.723527971391628015638e+2,
+		// 57!
+		1.763958484069973517152e+2,
+		// 58!
+		1.804562914175437710518e+2,
+		// 59!
+		1.845338288614494905025e+2,
+		// 60!
+		1.886281734236715911873e+2,
+		// 61!
+		1.927390472878449024360e+2,
+		// 62!
+		1.968661816728899939914e+2,
+		// 63!
+		2.010093163992815266793e+2,
+		// 64!
+		2.051681994826411985358e+2,
+		// 65!
+		2.093425867525368356464e+2,
+		// 66!
+		2.135322414945632611913e+2,
+		// 67!
+		2.177369341139542272510e+2,
+		// 68!
+		2.219564418191303339501e+2,
+		// 69!
+		2.261905483237275933323e+2,
+		// 70!
+		2.304390435657769523214e+2,
+		// 71!
+		2.347017234428182677427e+2,
+		// 72!
+		2.389783895618343230538e+2,
+		// 73!
+		2.432688490029827141829e+2,
+		// 74!
+		2.475729140961868839366e+2,
+		// 75!
+		2.518904022097231943772e+2,
+		// 76!
+		2.562211355500095254561e+2,
+		// 77!
+		2.605649409718632093053e+2,
+		// 78!
+		2.649216497985528010421e+2,
+		// 79!
+		2.692910976510198225363e+2,
+		// 80!
+		2.736731242856937041486e+2,
+		// 81!
+		2.780675734403661429141e+2,
+		// 82!
+		2.824742926876303960274e+2,
+		// 83!
+		2.868931332954269939509e+2,
+		// 84!
+		2.913239500942703075662e+2,
+		// 85!
+		2.957666013507606240211e+2,
+		// 86!
+		3.002209486470141317540e+2,
+		// 87!
+		3.046868567656687154726e+2,
+		// 88!
+		3.091641935801469219449e+2,
+		// 89!
+		3.136528299498790617832e+2,
+		// 90!
+		3.181526396202093268500e+2,
+		// 91!
+		3.226634991267261768912e+2,
+		// 92!
+		3.271852877037752172008e+2,
+		// 93!
+		3.317178871969284731381e+2,
+		// 94!
+		3.362611819791984770344e+2,
+		// 95!
+		3.408150588707990178690e+2,
+		// 96!
+		3.453794070622668541074e+2,
+		// 97!
+		3.499541180407702369296e+2,
+		// 98!
+		3.545390855194408088492e+2,
+		// 99!
+		3.591342053695753987760e+2,
+		// 100!
+		3.637393755555634901441e+2,
+		// 101!
+		3.683544960724047495950e+2,
+		// 102!
+		3.729794688856890206760e+2,
+		// 103!
+		3.776141978739186564468e+2,
+		// 104!
+		3.822585887730600291111e+2,
+		// 105!
+		3.869125491232175524822e+2,
+		// 106!
+		3.915759882173296196258e+2,
+		// 107!
+		3.962488170517915257991e+2,
+		// 108!
+		4.009309482789157454921e+2,
+		// 109!
+		4.056222961611448891925e+2,
+		// 110!
+		4.103227765269373054205e+2,
+		// 111!
+		4.150323067282496395563e+2,
+		// 112!
+		4.197508055995447340991e+2,
+		// 113!
+		4.244781934182570746677e+2,
+		// 114!
+		4.292143918666515701285e+2,
+		// 115!
+		4.339593239950148201939e+2,
+		// 116!
+		4.387129141861211848399e+2,
+		// 117!
+		4.434750881209189409588e+2,
+		// 118!
+		4.482457727453846057188e+2,
+		// 119!
+		4.530248962384961351041e+2,
+		// 120!
+		4.578123879812781810984e+2,
+		// 121!
+		4.626081785268749221865e+2,
+		// 122!
+		4.674121995716081787447e+2,
+		// 123!
+		4.722243839269805962399e+2,
+		// 124!
+		4.770446654925856331047e+2,
+		// 125!
+		4.818729792298879342285e+2,
+		// 126!
+		4.867092611368394122258e+2,
+		// 127!
+		4.915534482232980034989e+2,
+		// 128!
+		4.964054784872176206648e+2,
+		// 129!
+		5.012652908915792927797e+2,
+		// 130!
+		5.061328253420348751997e+2,
+		// 131!
+		5.110080226652360267439e+2,
+		// 132!
+		5.158908245878223975982e+2,
+		// 133!
+		5.207811737160441513633e+2,
+		// 134!
+		5.256790135159950627324e+2,
+		// 135!
+		5.305842882944334921812e+2,
+		// 136!
+		5.354969431801695441897e+2,
+		// 137!
+		5.404169241059976691050e+2,
+		// 138!
+		5.453441777911548737966e+2,
+		// 139!
+		5.502786517242855655538e+2,
+		// 140!
+		5.552202941468948698523e+2,
+		// 141!
+		5.601690540372730381305e+2,
+		// 142!
+		5.651248810948742988613e+2,
+		// 143!
+		5.700877257251342061414e+2,
+		// 144!
+		5.750575390247102067619e+2,
+		// 145!
+		5.800342727671307811636e+2,
+		// 146!
+		5.850178793888391176022e+2,
+		// 147!
+		5.900083119756178539038e+2,
+		// 148!
+		5.950055242493819689670e+2,
+		// 149!
+		6.000094705553274281080e+2,
+		// 150!
+		6.050201058494236838580e+2,
+		// 151!
+		6.100373856862386081868e+2,
+		// 152!
+		6.150612662070848845750e+2,
+		// 153!
+		6.200917041284773200381e+2,
+		// 154!
+		6.251286567308909491967e+2,
+		// 155!
+		6.301720818478101958172e+2,
+		// 156!
+		6.352219378550597328635e+2,
+		// 157!
+		6.402781836604080409209e+2,
+		// 158!
+		6.453407786934350077245e+2,
+		// 159!
+		6.504096828956552392500e+2,
+		// 160!
+		6.554848567108890661717e+2,
+		// 161!
+		6.605662610758735291676e+2,
+		// 162!
+		6.656538574111059132426e+2,
+		// 163!
+		6.707476076119126755767e+2,
+		// 164!
+		6.758474740397368739994e+2,
+		// 165!
+		6.809534195136374546094e+2,
+		// 166!
+		6.860654073019939978423e+2,
+		// 167!
+		6.911834011144107529496e+2,
+		// 168!
+		6.963073650938140118743e+2,
+		// 169!
+		7.014372638087370853465e+2,
+		// 170!
+		7.065730622457873471107e+2,
+		// 171!
+		7.117147258022900069535e+2,
+	]
+)
diff --git a/vlib/math/factorial_test.v b/vlib/math/factorial_test.v
new file mode 100644
index 000000000..77c204364
--- /dev/null
+++ b/vlib/math/factorial_test.v
@@ -0,0 +1,13 @@
+module math
+
+fn test_factorial() {
+	assert factorial(12) == 479001600
+	assert factorial(5) == 120
+	assert factorial(0) == 1
+}
+
+fn test_log_factorial() {
+	assert log_factorial(12) == log(479001600)
+	assert log_factorial(5) == log(120)
+	assert log_factorial(0) == log(1)
+}
diff --git a/vlib/math/floor.c.v b/vlib/math/floor.c.v
new file mode 100644
index 000000000..1dc933069
--- /dev/null
+++ b/vlib/math/floor.c.v
@@ -0,0 +1,34 @@
+module math
+
+fn C.ceil(x f64) f64
+
+fn C.floor(x f64) f64
+
+fn C.round(x f64) f64
+
+fn C.trunc(x f64) f64
+
+// ceil returns the nearest f64 greater or equal to the provided value.
+[inline]
+pub fn ceil(x f64) f64 {
+	return C.ceil(x)
+}
+
+// floor returns the nearest f64 lower or equal of the provided value.
+[inline]
+pub fn floor(x f64) f64 {
+	return C.floor(x)
+}
+
+// round returns the integer nearest to the provided value.
+[inline]
+pub fn round(x f64) f64 {
+	return C.round(x)
+}
+
+// trunc rounds a toward zero, returning the nearest integral value that is not
+// larger in magnitude than a.
+[inline]
+pub fn trunc(x f64) f64 {
+	return C.trunc(x)
+}
diff --git a/vlib/math/floor.js.v b/vlib/math/floor.js.v
new file mode 100644
index 000000000..6c995d18a
--- /dev/null
+++ b/vlib/math/floor.js.v
@@ -0,0 +1,34 @@
+module math
+
+fn JS.Math.ceil(x f64) f64
+
+fn JS.Math.floor(x f64) f64
+
+fn JS.Math.round(x f64) f64
+
+fn JS.Math.trunc(x f64) f64
+
+// ceil returns the nearest f64 greater or equal to the provided value.
+[inline]
+pub fn ceil(x f64) f64 {
+	return JS.Math.ceil(x)
+}
+
+// floor returns the nearest f64 lower or equal of the provided value.
+[inline]
+pub fn floor(x f64) f64 {
+	return JS.Math.floor(x)
+}
+
+// round returns the integer nearest to the provided value.
+[inline]
+pub fn round(x f64) f64 {
+	return JS.Math.round(x)
+}
+
+// trunc rounds a toward zero, returning the nearest integral value that is not
+// larger in magnitude than a.
+[inline]
+pub fn trunc(x f64) f64 {
+	return JS.Math.trunc(x)
+}
diff --git a/vlib/math/floor.v b/vlib/math/floor.v
new file mode 100644
index 000000000..a2c3208e3
--- /dev/null
+++ b/vlib/math/floor.v
@@ -0,0 +1,105 @@
+module math
+
+// floor returns the greatest integer value less than or equal to x.
+//
+// special cases are:
+// floor(±0) = ±0
+// floor(±inf) = ±inf
+// floor(nan) = nan
+pub fn floor(x f64) f64 {
+	if x == 0 || is_nan(x) || is_inf(x, 0) {
+		return x
+	}
+	if x < 0 {
+		mut d, fract := modf(-x)
+		if fract != 0.0 {
+			d = d + 1
+		}
+		return -d
+	}
+	d, _ := modf(x)
+	return d
+}
+
+// ceil returns the least integer value greater than or equal to x.
+//
+// special cases are:
+// ceil(±0) = ±0
+// ceil(±inf) = ±inf
+// ceil(nan) = nan
+pub fn ceil(x f64) f64 {
+	return -floor(-x)
+}
+
+// trunc returns the integer value of x.
+//
+// special cases are:
+// trunc(±0) = ±0
+// trunc(±inf) = ±inf
+// trunc(nan) = nan
+pub fn trunc(x f64) f64 {
+	if x == 0 || is_nan(x) || is_inf(x, 0) {
+		return x
+	}
+	d, _ := modf(x)
+	return d
+}
+
+// round returns the nearest integer, rounding half away from zero.
+//
+// special cases are:
+// round(±0) = ±0
+// round(±inf) = ±inf
+// round(nan) = nan
+pub fn round(x f64) f64 {
+	if x == 0 || is_nan(x) || is_inf(x, 0) {
+		return x
+	}
+	// Largest integer <= x
+	mut y := floor(x) // Fractional part
+	mut r := x - y // Round up to nearest.
+	if r > 0.5 {
+		unsafe {
+			goto rndup
+		}
+	}
+	// Round to even
+	if r == 0.5 {
+		r = y - 2.0 * floor(0.5 * y)
+		if r == 1.0 {
+			rndup:
+			y += 1.0
+		}
+	}
+	// Else round down.
+	return y
+}
+
+// round_to_even returns the nearest integer, rounding ties to even.
+//
+// special cases are:
+// round_to_even(±0) = ±0
+// round_to_even(±inf) = ±inf
+// round_to_even(nan) = nan
+pub fn round_to_even(x f64) f64 {
+	mut bits := f64_bits(x)
+	mut e_ := (bits >> shift) & mask
+	if e_ >= bias {
+		// round abs(x) >= 1.
+		// - Large numbers without fractional components, infinity, and nan are unchanged.
+		// - Add 0.499.. or 0.5 before truncating depending on whether the truncated
+		// number is even or odd (respectively).
+		half_minus_ulp := u64(u64(1) << (shift - 1)) - 1
+		e_ -= u64(bias)
+		bits += (half_minus_ulp + (bits >> (shift - e_)) & 1) >> e_
+		bits &= frac_mask >> e_
+		bits ^= frac_mask >> e_
+	} else if e_ == bias - 1 && bits & frac_mask != 0 {
+		// round 0.5 < abs(x) < 1.
+		bits = bits & sign_mask | uvone // +-1
+	} else {
+		// round abs(x) <= 0.5 including denormals.
+		bits &= sign_mask // +-0
+	}
+	return f64_from_bits(bits)
+}
diff --git a/vlib/math/gamma.c.v b/vlib/math/gamma.c.v
new file mode 100644
index 000000000..a4451cefa
--- /dev/null
+++ b/vlib/math/gamma.c.v
@@ -0,0 +1,17 @@
+module math
+
+fn C.tgamma(x f64) f64
+
+fn C.lgamma(x f64) f64
+
+// gamma computes the gamma function value
+[inline]
+pub fn gamma(a f64) f64 {
+	return C.tgamma(a)
+}
+
+// log_gamma computes the log-gamma function value
+[inline]
+pub fn log_gamma(x f64) f64 {
+	return C.lgamma(x)
+}
diff --git a/vlib/math/gamma.v b/vlib/math/gamma.v
new file mode 100644
index 000000000..e0061db39
--- /dev/null
+++ b/vlib/math/gamma.v
@@ -0,0 +1,335 @@
+module math
+
+// gamma function computed by Stirling's formula.
+// The pair of results must be multiplied together to get the actual answer.
+// The multiplication is left to the caller so that, if careful, the caller can avoid
+// infinity for 172 <= x <= 180.
+// The polynomial is valid for 33 <= x <= 172 larger values are only used
+// in reciprocal and produce denormalized floats. The lower precision there
+// masks any imprecision in the polynomial.
+fn stirling(x f64) (f64, f64) {
+	if x > 200 {
+		return inf(1), 1.0
+	}
+	sqrt_two_pi := 2.506628274631000502417
+	max_stirling := 143.01608
+	mut w := 1.0 / x
+	w = 1.0 + w * ((((gamma_s[0] * w + gamma_s[1]) * w + gamma_s[2]) * w + gamma_s[3]) * w +
+		gamma_s[4])
+	mut y1 := exp(x)
+	mut y2 := 1.0
+	if x > max_stirling { // avoid Pow() overflow
+		v := pow(x, 0.5 * x - 0.25)
+		y1_ := y1
+		y1 = v
+		y2 = v / y1_
+	} else {
+		y1 = pow(x, x - 0.5) / y1
+	}
+	return y1, f64(sqrt_two_pi) * w * y2
+}
+
+// gamma returns the gamma function of x.
+//
+// special ifs are:
+// gamma(+inf) = +inf
+// gamma(+0) = +inf
+// gamma(-0) = -inf
+// gamma(x) = nan for integer x < 0
+// gamma(-inf) = nan
+// gamma(nan) = nan
+pub fn gamma(a f64) f64 {
+	mut x := a
+	euler := 0.57721566490153286060651209008240243104215933593992 // A001620
+	if is_neg_int(x) || is_inf(x, -1) || is_nan(x) {
+		return nan()
+	}
+	if is_inf(x, 1) {
+		return inf(1)
+	}
+	if x == 0.0 {
+		return copysign(inf(1), x)
+	}
+	mut q := abs(x)
+	mut p := floor(q)
+	if q > 33 {
+		if x >= 0 {
+			y1, y2 := stirling(x)
+			return y1 * y2
+		}
+		// Note: x is negative but (checked above) not a negative integer,
+		// so x must be small enough to be in range for conversion to i64.
+		// If |x| were >= 2⁶³ it would have to be an integer.
+		mut signgam := 1
+		ip := i64(p)
+		if (ip & 1) == 0 {
+			signgam = -1
+		}
+		mut z := q - p
+		if z > 0.5 {
+			p = p + 1
+			z = q - p
+		}
+		z = q * sin(pi * z)
+		if z == 0 {
+			return inf(signgam)
+		}
+		sq1, sq2 := stirling(q)
+		absz := abs(z)
+		d := absz * sq1 * sq2
+		if is_inf(d, 0) {
+			z = pi / absz / sq1 / sq2
+		} else {
+			z = pi / d
+		}
+		return f64(signgam) * z
+	}
+	// Reduce argument
+	mut z := 1.0
+	for x >= 3 {
+		x = x - 1
+		z = z * x
+	}
+	for x < 0 {
+		if x > -1e-09 {
+			unsafe {
+				goto small
+			}
+		}
+		z = z / x
+		x = x + 1
+	}
+	for x < 2 {
+		if x < 1e-09 {
+			unsafe {
+				goto small
+			}
+		}
+		z = z / x
+		x = x + 1
+	}
+	if x == 2 {
+		return z
+	}
+	x = x - 2
+	p = (((((x * gamma_p[0] + gamma_p[1]) * x + gamma_p[2]) * x + gamma_p[3]) * x +
+		gamma_p[4]) * x + gamma_p[5]) * x + gamma_p[6]
+	q = ((((((x * gamma_q[0] + gamma_q[1]) * x + gamma_q[2]) * x + gamma_q[3]) * x +
+		gamma_q[4]) * x + gamma_q[5]) * x + gamma_q[6]) * x + gamma_q[7]
+	if true {
+		return z * p / q
+	}
+	small:
+	if x == 0 {
+		return inf(1)
+	}
+	return z / ((1.0 + euler * x) * x)
+}
+
+// log_gamma returns the natural logarithm and sign (-1 or +1) of Gamma(x).
+//
+// special ifs are:
+// log_gamma(+inf) = +inf
+// log_gamma(0) = +inf
+// log_gamma(-integer) = +inf
+// log_gamma(-inf) = -inf
+// log_gamma(nan) = nan
+pub fn log_gamma(x f64) f64 {
+	y, _ := log_gamma_sign(x)
+	return y
+}
+
+pub fn log_gamma_sign(a f64) (f64, int) {
+	mut x := a
+	ymin := 1.461632144968362245
+	tiny := exp2(-70)
+	two52 := exp2(52) // 0x4330000000000000 ~4.5036e+15
+	two58 := exp2(58) // 0x4390000000000000 ~2.8823e+17
+	tc := 1.46163214496836224576e+00 // 0x3FF762D86356BE3F
+	tf := -1.21486290535849611461e-01 // 0xBFBF19B9BCC38A42
+	// tt := -(tail of tf)
+	tt := -3.63867699703950536541e-18 // 0xBC50C7CAA48A971F
+	mut sign := 1
+	if is_nan(x) {
+		return x, sign
+	}
+	if is_inf(x, 1) {
+		return x, sign
+	}
+	if x == 0.0 {
+		return inf(1), sign
+	}
+	mut neg := false
+	if x < 0 {
+		x = -x
+		neg = true
+	}
+	if x < tiny { // if |x| < 2**-70, return -log(|x|)
+		if neg {
+			sign = -1
+		}
+		return -log(x), sign
+	}
+	mut nadj := 0.0
+	if neg {
+		if x >= two52 {
+			// x| >= 2**52, must be -integer
+			return inf(1), sign
+		}
+		t := sin_pi(x)
+		if t == 0 {
+			return inf(1), sign
+		}
+		nadj = log(pi / abs(t * x))
+		if t < 0 {
+			sign = -1
+		}
+	}
+	mut lgamma := 0.0
+	if x == 1 || x == 2 { // purge off 1 and 2
+		return 0.0, sign
+	} else if x < 2 { // use lgamma(x) = lgamma(x+1) - log(x)
+		mut y := 0.0
+		mut i := 0
+		if x <= 0.9 {
+			lgamma = -log(x)
+			if x >= (ymin - 1 + 0.27) { // 0.7316 <= x <=  0.9
+				y = 1.0 - x
+				i = 0
+			} else if x >= (ymin - 1 - 0.27) { // 0.2316 <= x < 0.7316
+				y = x - (tc - 1)
+				i = 1
+			} else { // 0 < x < 0.2316
+				y = x
+				i = 2
+			}
+		} else {
+			lgamma = 0
+			if x >= (ymin + 0.27) { // 1.7316 <= x < 2
+				y = f64(2) - x
+				i = 0
+			} else if x >= (ymin - 0.27) { // 1.2316 <= x < 1.7316
+				y = x - tc
+				i = 1
+			} else { // 0.9 < x < 1.2316
+				y = x - 1
+				i = 2
+			}
+		}
+		if i == 0 {
+			z := y * y
+			gamma_p1 := lgamma_a[0] + z * (lgamma_a[2] + z * (lgamma_a[4] + z * (lgamma_a[6] +
+				z * (lgamma_a[8] + z * lgamma_a[10]))))
+			gamma_p2 := z * (lgamma_a[1] + z * (lgamma_a[3] + z * (lgamma_a[5] + z * (lgamma_a[7] +
+				z * (lgamma_a[9] + z * lgamma_a[11])))))
+			p := y * gamma_p1 + gamma_p2
+			lgamma += (p - 0.5 * y)
+		} else if i == 1 {
+			z := y * y
+			w := z * y
+			gamma_p1 := lgamma_t[0] + w * (lgamma_t[3] + w * (lgamma_t[6] + w * (lgamma_t[9] +
+				w * lgamma_t[12]))) // parallel comp
+			gamma_p2 := lgamma_t[1] + w * (lgamma_t[4] + w * (lgamma_t[7] + w * (lgamma_t[10] +
+				w * lgamma_t[13])))
+			gamma_p3 := lgamma_t[2] + w * (lgamma_t[5] + w * (lgamma_t[8] + w * (lgamma_t[11] +
+				w * lgamma_t[14])))
+			p := z * gamma_p1 - (tt - w * (gamma_p2 + y * gamma_p3))
+			lgamma += (tf + p)
+		} else if i == 2 {
+			gamma_p1 := y * (lgamma_u[0] + y * (lgamma_u[1] + y * (lgamma_u[2] + y * (lgamma_u[3] +
+				y * (lgamma_u[4] + y * lgamma_u[5])))))
+			gamma_p2 := 1.0 + y * (lgamma_v[1] + y * (lgamma_v[2] + y * (lgamma_v[3] +
+				y * (lgamma_v[4] + y * lgamma_v[5]))))
+			lgamma += (-0.5 * y + gamma_p1 / gamma_p2)
+		}
+	} else if x < 8 { // 2 <= x < 8
+		i := int(x)
+		y := x - f64(i)
+		p := y * (lgamma_s[0] + y * (lgamma_s[1] + y * (lgamma_s[2] + y * (lgamma_s[3] +
+			y * (lgamma_s[4] + y * (lgamma_s[5] + y * lgamma_s[6]))))))
+		q := 1.0 + y * (lgamma_r[1] + y * (lgamma_r[2] + y * (lgamma_r[3] + y * (lgamma_r[4] +
+			y * (lgamma_r[5] + y * lgamma_r[6])))))
+		lgamma = 0.5 * y + p / q
+		mut z := 1.0 // lgamma(1+s) = log(s) + lgamma(s)
+		if i == 7 {
+			z *= (y + 6)
+			z *= (y + 5)
+			z *= (y + 4)
+			z *= (y + 3)
+			z *= (y + 2)
+			lgamma += log(z)
+		} else if i == 6 {
+			z *= (y + 5)
+			z *= (y + 4)
+			z *= (y + 3)
+			z *= (y + 2)
+			lgamma += log(z)
+		} else if i == 5 {
+			z *= (y + 4)
+			z *= (y + 3)
+			z *= (y + 2)
+			lgamma += log(z)
+		} else if i == 4 {
+			z *= (y + 3)
+			z *= (y + 2)
+			lgamma += log(z)
+		} else if i == 3 {
+			z *= (y + 2)
+			lgamma += log(z)
+		}
+	} else if x < two58 { // 8 <= x < 2**58
+		t := log(x)
+		z := 1.0 / x
+		y := z * z
+		w := lgamma_w[0] + z * (lgamma_w[1] + y * (lgamma_w[2] + y * (lgamma_w[3] +
+			y * (lgamma_w[4] + y * (lgamma_w[5] + y * lgamma_w[6])))))
+		lgamma = (x - 0.5) * (t - 1.0) + w
+	} else { // 2**58 <= x <= Inf
+		lgamma = x * (log(x) - 1.0)
+	}
+	if neg {
+		lgamma = nadj - lgamma
+	}
+	return lgamma, sign
+}
+
+// sin_pi(x) is a helper function for negative x
+fn sin_pi(x_ f64) f64 {
+	mut x := x_
+	two52 := exp2(52) // 0x4330000000000000 ~4.5036e+15
+	two53 := exp2(53) // 0x4340000000000000 ~9.0072e+15
+	if x < 0.25 {
+		return -sin(pi * x)
+	}
+	// argument reduction
+	mut z := floor(x)
+	mut n := 0
+	if z != x { // inexact
+		x = mod(x, 2)
+		n = int(x * 4)
+	} else {
+		if x >= two53 { // x must be even
+			x = 0
+			n = 0
+		} else {
+			if x < two52 {
+				z = x + two52 // exact
+			}
+			n = 1 & int(f64_bits(z))
+			x = f64(n)
+			n <<= 2
+		}
+	}
+	if n == 0 {
+		x = sin(pi * x)
+	} else if n == 1 || n == 2 {
+		x = cos(pi * (0.5 - x))
+	} else if n == 3 || n == 4 {
+		x = sin(pi * (1.0 - x))
+	} else if n == 5 || n == 6 {
+		x = -cos(pi * (x - 1.5))
+	} else {
+		x = sin(pi * (x - 2))
+	}
+	return -x
+}
diff --git a/vlib/math/gamma_tables.v b/vlib/math/gamma_tables.v
new file mode 100644
index 000000000..8bf9076dd
--- /dev/null
+++ b/vlib/math/gamma_tables.v
@@ -0,0 +1,163 @@
+module math
+
+const (
+	gamma_p  = [
+		1.60119522476751861407e-04,
+		1.19135147006586384913e-03,
+		1.04213797561761569935e-02,
+		4.76367800457137231464e-02,
+		2.07448227648435975150e-01,
+		4.94214826801497100753e-01,
+		9.99999999999999996796e-01,
+	]
+	gamma_q  = [
+		-2.31581873324120129819e-05,
+		5.39605580493303397842e-04,
+		-4.45641913851797240494e-03,
+		1.18139785222060435552e-02,
+		3.58236398605498653373e-02,
+		-2.34591795718243348568e-01,
+		7.14304917030273074085e-02,
+		1.00000000000000000320e+00,
+	]
+	gamma_s  = [
+		7.87311395793093628397e-04,
+		-2.29549961613378126380e-04,
+		-2.68132617805781232825e-03,
+		3.47222221605458667310e-03,
+		8.33333333333482257126e-02,
+	]
+	lgamma_a = [
+		// 0x3FB3C467E37DB0C8
+		7.72156649015328655494e-02,
+		// 0x3FD4A34CC4A60FAD
+		3.22467033424113591611e-01,
+		// 0x3FB13E001A5562A7
+		6.73523010531292681824e-02,
+		// 0x3F951322AC92547B
+		2.05808084325167332806e-02,
+		// 0x3F7E404FB68FEFE8
+		7.38555086081402883957e-03,
+		// 0x3F67ADD8CCB7926B
+		2.89051383673415629091e-03,
+		// 0x3F538A94116F3F5D
+		1.19270763183362067845e-03,
+		// 0x3F40B6C689B99C00
+		5.10069792153511336608e-04,
+		// 0x3F2CF2ECED10E54D
+		2.20862790713908385557e-04,
+		// 0x3F1C5088987DFB07
+		1.08011567247583939954e-04,
+		// 0x3EFA7074428CFA52
+		2.52144565451257326939e-05,
+		// 0x3F07858E90A45837
+		4.48640949618915160150e-05,
+	]
+	lgamma_r = [
+		// placeholder
+		1.0,
+		// 0x3FF645A762C4AB74
+		1.39200533467621045958e+00,
+		// 0x3FE71A1893D3DCDC
+		7.21935547567138069525e-01,
+		// 0x3FC601EDCCFBDF27
+		1.71933865632803078993e-01,
+		// 0x3F9317EA742ED475
+		1.86459191715652901344e-02,
+		// 0x3F497DDACA41A95B
+		7.77942496381893596434e-04,
+		// 0x3EDEBAF7A5B38140
+		7.32668430744625636189e-06,
+	]
+	lgamma_s = [
+		// 0xBFB3C467E37DB0C8
+		-7.72156649015328655494e-02,
+		// 0x3FCB848B36E20878
+		2.14982415960608852501e-01,
+		// 0x3FD4D98F4F139F59
+		3.25778796408930981787e-01,
+		// 0x3FC2BB9CBEE5F2F7
+		1.46350472652464452805e-01,
+		// 0x3F9B481C7E939961
+		2.66422703033638609560e-02,
+		// 0x3F5E26B67368F239
+		1.84028451407337715652e-03,
+		// 0x3F00BFECDD17E945
+		3.19475326584100867617e-05,
+	]
+	lgamma_t = [
+		// 0x3FDEF72BC8EE38A2
+		4.83836122723810047042e-01,
+		// 0xBFC2E4278DC6C509
+		-1.47587722994593911752e-01,
+		// 0x3FB08B4294D5419B
+		6.46249402391333854778e-02,
+		// 0xBFA0C9A8DF35B713
+		-3.27885410759859649565e-02,
+		// 0x3F9266E7970AF9EC
+		1.79706750811820387126e-02,
+		// 0xBF851F9FBA91EC6A
+		-1.03142241298341437450e-02,
+		// 0x3F78FCE0E370E344
+		6.10053870246291332635e-03,
+		// 0xBF6E2EFFB3E914D7
+		-3.68452016781138256760e-03,
+		// 0x3F6282D32E15C915
+		2.25964780900612472250e-03,
+		// 0xBF56FE8EBF2D1AF1
+		-1.40346469989232843813e-03,
+		// 0x3F4CDF0CEF61A8E9
+		8.81081882437654011382e-04,
+		// 0xBF41A6109C73E0EC
+		-5.38595305356740546715e-04,
+		// 0x3F34AF6D6C0EBBF7
+		3.15632070903625950361e-04,
+		// 0xBF347F24ECC38C38
+		-3.12754168375120860518e-04,
+		// 0x3F35FD3EE8C2D3F4
+		3.35529192635519073543e-04,
+	]
+	lgamma_u = [
+		// 0xBFB3C467E37DB0C8
+		-7.72156649015328655494e-02,
+		// 0x3FE4401E8B005DFF
+		6.32827064025093366517e-01,
+		// 0x3FF7475CD119BD6F
+		1.45492250137234768737e+00,
+		// 0x3FEF497644EA8450
+		9.77717527963372745603e-01,
+		// 0x3FCD4EAEF6010924
+		2.28963728064692451092e-01,
+		// 0x3F8B678BBF2BAB09
+		1.33810918536787660377e-02,
+	]
+	lgamma_v = [
+		1.0,
+		// 0x4003A5D7C2BD619C
+		2.45597793713041134822e+00,
+		// 0x40010725A42B18F5
+		2.12848976379893395361e+00,
+		// 0x3FE89DFBE45050AF
+		7.69285150456672783825e-01,
+		// 0x3FBAAE55D6537C88
+		1.04222645593369134254e-01,
+		// 0x3F6A5ABB57D0CF61
+		3.21709242282423911810e-03,
+	]
+	lgamma_w = [
+		// 0x3FDACFE390C97D69
+		4.18938533204672725052e-01,
+		// 0x3FB555555555553B
+		8.33333333333329678849e-02,
+		// 0xBF66C16C16B02E5C
+		-2.77777777728775536470e-03,
+		// 0x3F4A019F98CF38B6
+		7.93650558643019558500e-04,
+		// 0xBF4380CB8C0FE741
+		-5.95187557450339963135e-04,
+		// 0x3F4B67BA4CDAD5D1
+		8.36339918996282139126e-04,
+		// 0xBF5AB89D0B9E43E4
+		-1.63092934096575273989e-03,
+	]
+)
diff --git a/vlib/math/hypot.c.v b/vlib/math/hypot.c.v
new file mode 100644
index 000000000..9710ea4a1
--- /dev/null
+++ b/vlib/math/hypot.c.v
@@ -0,0 +1,9 @@
+module math
+
+fn C.hypot(x f64, y f64) f64
+
+// Returns hypotenuse of a right triangle.
+[inline]
+pub fn hypot(x f64, y f64) f64 {
+	return C.hypot(x, y)
+}
diff --git a/vlib/math/hypot.v b/vlib/math/hypot.v
new file mode 100644
index 000000000..41375057b
--- /dev/null
+++ b/vlib/math/hypot.v
@@ -0,0 +1,24 @@
+module math
+
+pub fn hypot(x f64, y f64) f64 {
+	if is_inf(x, 0) || is_inf(y, 0) {
+		return inf(1)
+	}
+	if is_nan(x) || is_nan(y) {
+		return nan()
+	}
+	mut result := 0.0
+	if x != 0.0 || y != 0.0 {
+		abs_x := abs(x)
+		abs_y := abs(y)
+		min, max := minmax(abs_x, abs_y)
+		rat := min / max
+		root_term := sqrt(1.0 + rat * rat)
+		if max < max_f64 / root_term {
+			result = max * root_term
+		} else {
+			panic('overflow in hypot_e function')
+		}
+	}
+	return result
+}
diff --git a/vlib/math/internal/machine.v b/vlib/math/internal/machine.v
new file mode 100644
index 000000000..e4ffb62f7
--- /dev/null
+++ b/vlib/math/internal/machine.v
@@ -0,0 +1,58 @@
+module internal
+
+// contants to do fine tuning of precision for the functions
+// implemented in pure V
+pub const (
+	f64_epsilon           = 2.2204460492503131e-16
+	sqrt_f64_epsilon      = 1.4901161193847656e-08
+	root3_f64_epsilon     = 6.0554544523933429e-06
+	root4_f64_epsilon     = 1.2207031250000000e-04
+	root5_f64_epsilon     = 7.4009597974140505e-04
+	root6_f64_epsilon     = 2.4607833005759251e-03
+	log_f64_epsilon       = -3.6043653389117154e+01
+	f64_min               = 2.2250738585072014e-308
+	sqrt_f64_min          = 1.4916681462400413e-154
+	root3_f64_min         = 2.8126442852362996e-103
+	root4_f64_min         = 1.2213386697554620e-77
+	root5_f64_min         = 2.9476022969691763e-62
+	root6_f64_min         = 5.3034368905798218e-52
+	log_f64_min           = -7.0839641853226408e+02
+	f64_max               = 1.7976931348623157e+308
+	sqrt_f64_max          = 1.3407807929942596e+154
+	root3_f64_max         = 5.6438030941222897e+102
+	root4_f64_max         = 1.1579208923731620e+77
+	root5_f64_max         = 4.4765466227572707e+61
+	root6_f64_max         = 2.3756689782295612e+51
+	log_f64_max           = 7.0978271289338397e+02
+	f32_epsilon           = 1.1920928955078125e-07
+	sqrt_f32_epsilon      = 3.4526698300124393e-04
+	root3_f32_epsilon     = 4.9215666011518501e-03
+	root4_f32_epsilon     = 1.8581361171917516e-02
+	root5_f32_epsilon     = 4.1234622211652937e-02
+	root6_f32_epsilon     = 7.0153878019335827e-02
+	log_f32_epsilon       = -1.5942385152878742e+01
+	f32_min               = 1.1754943508222875e-38
+	sqrt_f32_min          = 1.0842021724855044e-19
+	root3_f32_min         = 2.2737367544323241e-13
+	root4_f32_min         = 3.2927225399135965e-10
+	root5_f32_min         = 2.5944428542140822e-08
+	root6_f32_min         = 4.7683715820312542e-07
+	log_f32_min           = -8.7336544750553102e+01
+	f32_max               = 3.4028234663852886e+38
+	sqrt_f32_max          = 1.8446743523953730e+19
+	root3_f32_max         = 6.9814635196223242e+12
+	root4_f32_max         = 4.2949672319999986e+09
+	root5_f32_max         = 5.0859007855960041e+07
+	root6_f32_max         = 2.6422459233807749e+06
+	log_f32_max           = 8.8722839052068352e+01
+	sflt_epsilon          = 4.8828125000000000e-04
+	sqrt_sflt_epsilon     = 2.2097086912079612e-02
+	root3_sflt_epsilon    = 7.8745065618429588e-02
+	root4_sflt_epsilon    = 1.4865088937534013e-01
+	root5_sflt_epsilon    = 2.1763764082403100e-01
+	root6_sflt_epsilon    = 2.8061551207734325e-01
+	log_sflt_epsilon      = -7.6246189861593985e+00
+	max_int_fact_arg      = 170
+	max_f64_fact_arg      = 171.0
+	max_long_f64_fact_arg = 1755.5
+)
diff --git a/vlib/math/invhyp.v b/vlib/math/invhyp.v
new file mode 100644
index 000000000..4ef93cbce
--- /dev/null
+++ b/vlib/math/invhyp.v
@@ -0,0 +1,51 @@
+module math
+
+import math.internal
+
+pub fn acosh(x f64) f64 {
+	if x == 0.0 {
+		return 0.0
+	} else if x > 1.0 / internal.sqrt_f64_epsilon {
+		return log(x) + pi * 2
+	} else if x > 2.0 {
+		return log(2.0 * x - 1.0 / (sqrt(x * x - 1.0) + x))
+	} else if x > 1.0 {
+		t := x - 1.0
+		return log1p(t + sqrt(2.0 * t + t * t))
+	} else if x == 1.0 {
+		return 0.0
+	} else {
+		return nan()
+	}
+}
+
+pub fn asinh(x f64) f64 {
+	a := abs(x)
+	s := if x < 0 { -1.0 } else { 1.0 }
+	if a > 1.0 / internal.sqrt_f64_epsilon {
+		return s * (log(a) + pi * 2.0)
+	} else if a > 2.0 {
+		return s * log(2.0 * a + 1.0 / (a + sqrt(a * a + 1.0)))
+	} else if a > internal.sqrt_f64_epsilon {
+		a2 := a * a
+		return s * log1p(a + a2 / (1.0 + sqrt(1.0 + a2)))
+	} else {
+		return x
+	}
+}
+
+pub fn atanh(x f64) f64 {
+	a := abs(x)
+	s := if x < 0 { -1.0 } else { 1.0 }
+	if a > 1.0 {
+		return nan()
+	} else if a == 1.0 {
+		return if x < 0 { inf(-1) } else { inf(1) }
+	} else if a >= 0.5 {
+		return s * 0.5 * log1p(2.0 * a / (1.0 - a))
+	} else if a > internal.f64_epsilon {
+		return s * 0.5 * log1p(2.0 * a + 2.0 * a * a / (1.0 - a))
+	} else {
+		return x
+	}
+}
diff --git a/vlib/math/invtrig.c.v b/vlib/math/invtrig.c.v
new file mode 100644
index 000000000..ee4caf6a7
--- /dev/null
+++ b/vlib/math/invtrig.c.v
@@ -0,0 +1,33 @@
+module math
+
+fn C.acos(x f64) f64
+
+fn C.asin(x f64) f64
+
+fn C.atan(x f64) f64
+
+fn C.atan2(y f64, x f64) f64
+
+// acos calculates inverse cosine (arccosine).
+[inline]
+pub fn acos(a f64) f64 {
+	return C.acos(a)
+}
+
+// asin calculates inverse sine (arcsine).
+[inline]
+pub fn asin(a f64) f64 {
+	return C.asin(a)
+}
+
+// atan calculates inverse tangent (arctangent).
+[inline]
+pub fn atan(a f64) f64 {
+	return C.atan(a)
+}
+
+// atan2 calculates inverse tangent with two arguments, returns the angle between the X axis and the point.
+[inline]
+pub fn atan2(a f64, b f64) f64 {
+	return C.atan2(a, b)
+}
diff --git a/vlib/math/invtrig.js.v b/vlib/math/invtrig.js.v
new file mode 100644
index 000000000..2e280a5a0
--- /dev/null
+++ b/vlib/math/invtrig.js.v
@@ -0,0 +1,33 @@
+module math
+
+fn JS.Math.acos(x f64) f64
+
+fn JS.Math.asin(x f64) f64
+
+fn JS.Math.atan(x f64) f64
+
+fn JS.Math.atan2(y f64, x f64) f64
+
+// acos calculates inverse cosine (arccosine).
+[inline]
+pub fn acos(a f64) f64 {
+	return JS.Math.acos(a)
+}
+
+// asin calculates inverse sine (arcsine).
+[inline]
+pub fn asin(a f64) f64 {
+	return JS.Math.asin(a)
+}
+
+// atan calculates inverse tangent (arctangent).
+[inline]
+pub fn atan(a f64) f64 {
+	return JS.Math.atan(a)
+}
+
+// atan2 calculates inverse tangent with two arguments, returns the angle between the X axis and the point.
+[inline]
+pub fn atan2(a f64, b f64) f64 {
+	return JS.Math.atan2(a, b)
+}
diff --git a/vlib/math/invtrig.v b/vlib/math/invtrig.v
new file mode 100644
index 000000000..c40d87f8d
--- /dev/null
+++ b/vlib/math/invtrig.v
@@ -0,0 +1,219 @@
+module math
+
+// The original C code, the long comment, and the constants below were
+// from http://netlib.sandia.gov/cephes/cmath/atan.c, available from
+// http://www.netlib.org/cephes/ctgz.
+// The go code is a version of the original C.
+//
+// atan.c
+// Inverse circular tangent (arctangent)
+//
+// SYNOPSIS:
+// double x, y, atan()
+// y = atan( x )
+//
+// DESCRIPTION:
+// Returns radian angle between -pi/2.0 and +pi/2.0 whose tangent is x.
+//
+// Range reduction is from three intervals into the interval from zero to 0.66.
+// The approximant uses a rational function of degree 4/5 of the form
+// x + x**3 P(x)/Q(x).
+//
+// ACCURACY:
+// Relative error:
+// arithmetic   domain    # trials  peak     rms
+// DEC       -10, 10   50000     2.4e-17  8.3e-18
+// IEEE      -10, 10   10^6      1.8e-16  5.0e-17
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+// Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+// The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+// Stephen L. Moshier
+// moshier@na-net.ornl.gov
+// pi/2.0 = PIO2 + morebits
+// tan3pio8 = tan(3*pi/8)
+
+const (
+	morebits = 6.123233995736765886130e-17
+	tan3pio8 = 2.41421356237309504880
+)
+
+// xatan evaluates a series valid in the range [0, 0.66].
+[inline]
+fn xatan(x f64) f64 {
+	xatan_p0 := -8.750608600031904122785e-01
+	xatan_p1 := -1.615753718733365076637e+01
+	xatan_p2 := -7.500855792314704667340e+01
+	xatan_p3 := -1.228866684490136173410e+02
+	xatan_p4 := -6.485021904942025371773e+01
+	xatan_q0 := 2.485846490142306297962e+01
+	xatan_q1 := 1.650270098316988542046e+02
+	xatan_q2 := 4.328810604912902668951e+02
+	xatan_q3 := 4.853903996359136964868e+02
+	xatan_q4 := 1.945506571482613964425e+02
+	mut z := x * x
+	z = z * ((((xatan_p0 * z + xatan_p1) * z + xatan_p2) * z + xatan_p3) * z + xatan_p4) / (((((z +
+		xatan_q0) * z + xatan_q1) * z + xatan_q2) * z + xatan_q3) * z + xatan_q4)
+	z = x * z + x
+	return z
+}
+
+// satan reduces its argument (known to be positive)
+// to the range [0, 0.66] and calls xatan.
+[inline]
+fn satan(x f64) f64 {
+	if x <= 0.66 {
+		return xatan(x)
+	}
+	if x > math.tan3pio8 {
+		return pi / 2.0 - xatan(1.0 / x) + f64(math.morebits)
+	}
+	return pi / 4 + xatan((x - 1.0) / (x + 1.0)) + 0.5 * f64(math.morebits)
+}
+
+// atan returns the arctangent, in radians, of x.
+//
+// special cases are:
+// atan(±0) = ±0
+// atan(±inf) = ±pi/2.0
+pub fn atan(x f64) f64 {
+	if x == 0 {
+		return x
+	}
+	if x > 0 {
+		return satan(x)
+	}
+	return -satan(-x)
+}
+
+// atan2 returns the arc tangent of y/x, using
+// the signs of the two to determine the quadrant
+// of the return value.
+//
+// special cases are (in order):
+// 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.0
+// atan2(y<0, 0) = -pi/2.0
+// 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.0
+// atan2(-inf, x) = -pi/2.0
+pub fn atan2(y f64, x f64) f64 {
+	// special cases
+	if is_nan(y) || is_nan(x) {
+		return nan()
+	}
+	if y == 0.0 {
+		if x >= 0 && !signbit(x) {
+			return copysign(0, y)
+		}
+		return copysign(pi, y)
+	}
+	if x == 0.0 {
+		return copysign(pi / 2.0, y)
+	}
+	if is_inf(x, 0) {
+		if is_inf(x, 1) {
+			if is_inf(y, 0) {
+				return copysign(pi / 4, y)
+			}
+			return copysign(0, y)
+		}
+		if is_inf(y, 0) {
+			return copysign(3.0 * pi / 4.0, y)
+		}
+		return copysign(pi, y)
+	}
+	if is_inf(y, 0) {
+		return copysign(pi / 2.0, y)
+	}
+	// Call atan and determine the quadrant.
+	q := atan(y / x)
+	if x < 0 {
+		if q <= 0 {
+			return q + pi
+		}
+		return q - pi
+	}
+	return q
+}
+
+/*
+Floating-point arcsine and arccosine.
+
+	They are implemented by computing the arctangent
+	after appropriate range reduction.
+*/
+
+// asin returns the arcsine, in radians, of x.
+//
+// special cases are:
+// asin(±0) = ±0
+// asin(x) = nan if x < -1 or x > 1
+pub fn asin(x_ f64) f64 {
+	mut x := x_
+	if x == 0.0 {
+		return x // special case
+	}
+	mut sign := false
+	if x < 0.0 {
+		x = -x
+		sign = true
+	}
+	if x > 1.0 {
+		return nan() // special case
+	}
+	mut temp := sqrt(1.0 - x * x)
+	if x > 0.7 {
+		temp = pi / 2.0 - satan(temp / x)
+	} else {
+		temp = satan(x / temp)
+	}
+	if sign {
+		temp = -temp
+	}
+	return temp
+}
+
+// acos returns the arccosine, in radians, of x.
+//
+// special case is:
+// acos(x) = nan if x < -1 or x > 1
+[inline]
+pub fn acos(x f64) f64 {
+	if (x < -1.0) || (x > 1.0) {
+		return nan()
+	}
+	if x == 0.0 {
+		return nan()
+	}
+	if x > 0.5 {
+		return f64(2.0) * asin(sqrt(0.5 - 0.5 * x))
+	}
+	mut z := pi / f64(4.0) - asin(x)
+	z = z + math.morebits
+	z = z + pi / f64(4.0)
+	return z
+}
diff --git a/vlib/math/log.c.v b/vlib/math/log.c.v
new file mode 100644
index 000000000..f2f8c9397
--- /dev/null
+++ b/vlib/math/log.c.v
@@ -0,0 +1,25 @@
+module math
+
+fn C.log(x f64) f64
+
+fn C.log2(x f64) f64
+
+fn C.log10(x f64) f64
+
+// log calculates natural (base-e) logarithm of the provided value.
+[inline]
+pub fn log(x f64) f64 {
+	return C.log(x)
+}
+
+// log2 calculates base-2 logarithm of the provided value.
+[inline]
+pub fn log2(x f64) f64 {
+	return C.log2(x)
+}
+
+// log10 calculates the common (base-10) logarithm of the provided value.
+[inline]
+pub fn log10(x f64) f64 {
+	return C.log10(x)
+}
diff --git a/vlib/math/log.js.v b/vlib/math/log.js.v
new file mode 100644
index 000000000..a1a0cbb18
--- /dev/null
+++ b/vlib/math/log.js.v
@@ -0,0 +1,9 @@
+module math
+
+fn JS.Math.log(x f64) f64
+
+// log calculates natural (base-e) logarithm of the provided value.
+[inline]
+pub fn log(x f64) f64 {
+	return JS.Math.log(x)
+}
diff --git a/vlib/math/log.v b/vlib/math/log.v
new file mode 100644
index 000000000..47ef73122
--- /dev/null
+++ b/vlib/math/log.v
@@ -0,0 +1,76 @@
+module math
+
+pub fn log_n(x f64, b f64) f64 {
+	y := log(x)
+	z := log(b)
+	return y / z
+}
+
+// log10 returns the decimal logarithm of x.
+// The special cases are the same as for log.
+pub fn log10(x f64) f64 {
+	return log(x) * (1.0 / ln10)
+}
+
+// log2 returns the binary logarithm of x.
+// The special cases are the same as for log.
+pub fn log2(x f64) f64 {
+	frac, exp := frexp(x)
+	// Make sure exact powers of two give an exact answer.
+	// Don't depend on log(0.5)*(1/ln2)+exp being exactly exp-1.
+	if frac == 0.5 {
+		return f64(exp - 1)
+	}
+	return log(frac) * (1.0 / ln2) + f64(exp)
+}
+
+pub fn log1p(x f64) f64 {
+	y := 1.0 + x
+	z := y - 1.0
+	return log(y) - (z - x) / y // cancels errors with IEEE arithmetic
+}
+
+// log_b returns the binary exponent of x.
+//
+// special cases are:
+// log_b(±inf) = +inf
+// log_b(0) = -inf
+// log_b(nan) = nan
+pub fn log_b(x f64) f64 {
+	if x == 0 {
+		return inf(-1)
+	}
+	if is_inf(x, 0) {
+		return inf(1)
+	}
+	if is_nan(x) {
+		return x
+	}
+	return f64(ilog_b_(x))
+}
+
+// ilog_b returns the binary exponent of x as an integer.
+//
+// special cases are:
+// ilog_b(±inf) = max_i32
+// ilog_b(0) = min_i32
+// ilog_b(nan) = max_i32
+pub fn ilog_b(x f64) int {
+	if x == 0 {
+		return min_i32
+	}
+	if is_nan(x) {
+		return max_i32
+	}
+	if is_inf(x, 0) {
+		return max_i32
+	}
+	return ilog_b_(x)
+}
+
+// ilog_b returns the binary exponent of x. It assumes x is finite and
+// non-zero.
+fn ilog_b_(x_ f64) int {
+	x, exp := normalize(x_)
+	return int((f64_bits(x) >> shift) & mask) - bias + exp
+}
diff --git a/vlib/math/math.c.v b/vlib/math/math.c.v
index 5d605d80c..503b6ab8b 100644
--- a/vlib/math/math.c.v
+++ b/vlib/math/math.c.v
@@ -12,287 +12,3 @@ $if windows {
 } $else {
 	#flag -lm
 }
-
-fn C.acos(x f64) f64
-
-fn C.asin(x f64) f64
-
-fn C.atan(x f64) f64
-
-fn C.atan2(y f64, x f64) f64
-
-fn C.cbrt(x f64) f64
-
-fn C.ceil(x f64) f64
-
-fn C.cos(x f64) f64
-
-fn C.cosf(x f32) f32
-
-fn C.cosh(x f64) f64
-
-fn C.erf(x f64) f64
-
-fn C.erfc(x f64) f64
-
-fn C.exp(x f64) f64
-
-fn C.exp2(x f64) f64
-
-fn C.fabs(x f64) f64
-
-fn C.floor(x f64) f64
-
-fn C.fmod(x f64, y f64) f64
-
-fn C.hypot(x f64, y f64) f64
-
-fn C.log(x f64) f64
-
-fn C.log2(x f64) f64
-
-fn C.log10(x f64) f64
-
-fn C.lgamma(x f64) f64
-
-fn C.pow(x f64, y f64) f64
-
-fn C.powf(x f32, y f32) f32
-
-fn C.round(x f64) f64
-
-fn C.sin(x f64) f64
-
-fn C.sinf(x f32) f32
-
-fn C.sinh(x f64) f64
-
-fn C.sqrt(x f64) f64
-
-fn C.sqrtf(x f32) f32
-
-fn C.tgamma(x f64) f64
-
-fn C.tan(x f64) f64
-
-fn C.tanf(x f32) f32
-
-fn C.tanh(x f64) f64
-
-fn C.trunc(x f64) f64
-
-// NOTE
-// When adding a new function, please make sure it's in the right place.
-// All functions are sorted alphabetically.
-// Returns the absolute value.
-[inline]
-pub fn abs(a f64) f64 {
-	return C.fabs(a)
-}
-
-// acos calculates inverse cosine (arccosine).
-[inline]
-pub fn acos(a f64) f64 {
-	return C.acos(a)
-}
-
-// asin calculates inverse sine (arcsine).
-[inline]
-pub fn asin(a f64) f64 {
-	return C.asin(a)
-}
-
-// atan calculates inverse tangent (arctangent).
-[inline]
-pub fn atan(a f64) f64 {
-	return C.atan(a)
-}
-
-// atan2 calculates inverse tangent with two arguments, returns the angle between the X axis and the point.
-[inline]
-pub fn atan2(a f64, b f64) f64 {
-	return C.atan2(a, b)
-}
-
-// cbrt calculates cubic root.
-[inline]
-pub fn cbrt(a f64) f64 {
-	return C.cbrt(a)
-}
-
-// ceil returns the nearest f64 greater or equal to the provided value.
-[inline]
-pub fn ceil(a f64) f64 {
-	return C.ceil(a)
-}
-
-// cos calculates cosine.
-[inline]
-pub fn cos(a f64) f64 {
-	return C.cos(a)
-}
-
-// cosf calculates cosine. (float32)
-[inline]
-pub fn cosf(a f32) f32 {
-	return C.cosf(a)
-}
-
-// cosh calculates hyperbolic cosine.
-[inline]
-pub fn cosh(a f64) f64 {
-	return C.cosh(a)
-}
-
-// exp calculates exponent of the number (math.pow(math.E, a)).
-[inline]
-pub fn exp(a f64) f64 {
-	return C.exp(a)
-}
-
-/*
-// erf computes the error function value
-[inline]
-pub fn erf(a f64) f64 {
-	return C.erf(a)
-}
-*/
-/*
-// erfc computes the complementary error function value
-[inline]
-pub fn erfc(a f64) f64 {
-	return C.erfc(a)
-}
-*/
-// exp2 returns the base-2 exponential function of a (math.pow(2, a)).
-[inline]
-pub fn exp2(a f64) f64 {
-	return C.exp2(a)
-}
-
-// floor returns the nearest f64 lower or equal of the provided value.
-[inline]
-pub fn floor(a f64) f64 {
-	return C.floor(a)
-}
-
-// fmod returns the floating-point remainder of number / denom (rounded towards zero):
-[inline]
-pub fn fmod(a f64, b f64) f64 {
-	return C.fmod(a, b)
-}
-
-// gamma computes the gamma function value
-[inline]
-pub fn gamma(a f64) f64 {
-	return C.tgamma(a)
-}
-
-// Returns hypotenuse of a right triangle.
-[inline]
-pub fn hypot(a f64, b f64) f64 {
-	return C.hypot(a, b)
-}
-
-// log calculates natural (base-e) logarithm of the provided value.
-[inline]
-pub fn log(a f64) f64 {
-	return C.log(a)
-}
-
-// log2 calculates base-2 logarithm of the provided value.
-[inline]
-pub fn log2(a f64) f64 {
-	return C.log2(a)
-}
-
-// log10 calculates the common (base-10) logarithm of the provided value.
-[inline]
-pub fn log10(a f64) f64 {
-	return C.log10(a)
-}
-
-// log_gamma computes the log-gamma function value
-[inline]
-pub fn log_gamma(a f64) f64 {
-	return C.lgamma(a)
-}
-
-// log_n calculates base-N logarithm of the provided value.
-[inline]
-pub fn log_n(a f64, b f64) f64 {
-	return C.log(a) / C.log(b)
-}
-
-// pow returns base raised to the provided power.
-[inline]
-pub fn pow(a f64, b f64) f64 {
-	return C.pow(a, b)
-}
-
-// powf returns base raised to the provided power. (float32)
-[inline]
-pub fn powf(a f32, b f32) f32 {
-	return C.powf(a, b)
-}
-
-// round returns the integer nearest to the provided value.
-[inline]
-pub fn round(f f64) f64 {
-	return C.round(f)
-}
-
-// sin calculates sine.
-[inline]
-pub fn sin(a f64) f64 {
-	return C.sin(a)
-}
-
-// sinf calculates sine. (float32)
-[inline]
-pub fn sinf(a f32) f32 {
-	return C.sinf(a)
-}
-
-// sinh calculates hyperbolic sine.
-[inline]
-pub fn sinh(a f64) f64 {
-	return C.sinh(a)
-}
-
-// sqrt calculates square-root of the provided value.
-[inline]
-pub fn sqrt(a f64) f64 {
-	return C.sqrt(a)
-}
-
-// sqrtf calculates square-root of the provided value. (float32)
-[inline]
-pub fn sqrtf(a f32) f32 {
-	return C.sqrtf(a)
-}
-
-// tan calculates tangent.
-[inline]
-pub fn tan(a f64) f64 {
-	return C.tan(a)
-}
-
-// tanf calculates tangent. (float32)
-[inline]
-pub fn tanf(a f32) f32 {
-	return C.tanf(a)
-}
-
-// tanh calculates hyperbolic tangent.
-[inline]
-pub fn tanh(a f64) f64 {
-	return C.tanh(a)
-}
-
-// trunc rounds a toward zero, returning the nearest integral value that is not
-// larger in magnitude than a.
-[inline]
-pub fn trunc(a f64) f64 {
-	return C.trunc(a)
-}
diff --git a/vlib/math/math.js.v b/vlib/math/math.js.v
deleted file mode 100644
index 3c8ce8de9..000000000
--- a/vlib/math/math.js.v
+++ /dev/null
@@ -1,281 +0,0 @@
-// Copyright (c) 2019-2021 Alexander Medvednikov. All rights reserved.
-// Use of this source code is governed by an MIT license
-// that can be found in the LICENSE file.
-module math
-
-// TODO : The commented out functions need either a native V implementation, a
-// JS specific implementation, or use some other JS math library, such as
-// https://github.com/josdejong/mathjs
-
-// Replaces C.fabs
-fn JS.Math.abs(x f64) f64
-
-fn JS.Math.acos(x f64) f64
-fn JS.Math.asin(x f64) f64
-fn JS.Math.atan(x f64) f64
-fn JS.Math.atan2(y f64, x f64) f64
-fn JS.Math.cbrt(x f64) f64
-fn JS.Math.ceil(x f64) f64
-fn JS.Math.cos(x f64) f64
-fn JS.Math.cosh(x f64) f64
-
-// fn JS.Math.erf(x f64) f64 // Not in standard JS Math object
-// fn JS.Math.erfc(x f64) f64 // Not in standard JS Math object
-fn JS.Math.exp(x f64) f64
-
-// fn JS.Math.exp2(x f64) f64 // Not in standard JS Math object
-fn JS.Math.floor(x f64) f64
-
-// fn JS.Math.fmod(x f64, y f64) f64 // Not in standard JS Math object
-// fn JS.Math.hypot(x f64, y f64) f64 // Not in standard JS Math object
-fn JS.Math.log(x f64) f64
-
-// fn JS.Math.log2(x f64) f64 // Not in standard JS Math object
-// fn JS.Math.log10(x f64) f64 // Not in standard JS Math object
-// fn JS.Math.lgamma(x f64) f64 // Not in standard JS Math object
-fn JS.Math.pow(x f64, y f64) f64
-fn JS.Math.round(x f64) f64
-fn JS.Math.sin(x f64) f64
-fn JS.Math.sinh(x f64) f64
-fn JS.Math.sqrt(x f64) f64
-
-// fn JS.Math.tgamma(x f64) f64 // Not in standard JS Math object
-fn JS.Math.tan(x f64) f64
-fn JS.Math.tanh(x f64) f64
-fn JS.Math.trunc(x f64) f64
-
-// NOTE
-// When adding a new function, please make sure it's in the right place.
-// All functions are sorted alphabetically.
-
-// Returns the absolute value.
-[inline]
-pub fn abs(a f64) f64 {
-	return JS.Math.abs(a)
-}
-
-// acos calculates inverse cosine (arccosine).
-[inline]
-pub fn acos(a f64) f64 {
-	return JS.Math.acos(a)
-}
-
-// asin calculates inverse sine (arcsine).
-[inline]
-pub fn asin(a f64) f64 {
-	return JS.Math.asin(a)
-}
-
-// atan calculates inverse tangent (arctangent).
-[inline]
-pub fn atan(a f64) f64 {
-	return JS.Math.atan(a)
-}
-
-// atan2 calculates inverse tangent with two arguments, returns the angle between the X axis and the point.
-[inline]
-pub fn atan2(a f64, b f64) f64 {
-	return JS.Math.atan2(a, b)
-}
-
-// cbrt calculates cubic root.
-[inline]
-pub fn cbrt(a f64) f64 {
-	return JS.Math.cbrt(a)
-}
-
-// ceil returns the nearest f64 greater or equal to the provided value.
-[inline]
-pub fn ceil(a f64) f64 {
-	return JS.Math.ceil(a)
-}
-
-// cos calculates cosine.
-[inline]
-pub fn cos(a f64) f64 {
-	return JS.Math.cos(a)
-}
-
-// cosf calculates cosine. (float32). This doesn't exist in JS
-[inline]
-pub fn cosf(a f32) f32 {
-	return f32(JS.Math.cos(a))
-}
-
-// cosh calculates hyperbolic cosine.
-[inline]
-pub fn cosh(a f64) f64 {
-	return JS.Math.cosh(a)
-}
-
-// exp calculates exponent of the number (math.pow(math.E, a)).
-[inline]
-pub fn exp(a f64) f64 {
-	mut res := 0.0
-	#res.val = Math.exp(a)
-
-	return res
-}
-
-// exp2 returns the base-2 exponential function of a (math.pow(2, a)).
-[inline]
-pub fn exp2(a f64) f64 {
-	return 0
-	// return JS.Math.exp2(a)
-}
-
-// floor returns the nearest f64 lower or equal of the provided value.
-[inline]
-pub fn floor(a f64) f64 {
-	return JS.Math.floor(a)
-}
-
-// fmod returns the floating-point remainder of number / denom (rounded towards zero):
-[inline]
-pub fn fmod(x f64, y f64) f64 {
-	#let tmp
-	#let tmp2
-	#let p = 0
-	#let pY = 0
-	#let l = 0.0
-	#let l2 = 0.0
-	#tmp = x.toExponential().match(/^.\.?(.*)e(.+)$/)
-	#p = parseInt(tmp[2], 10) - (tmp[1] + '').length
-	#tmp = y.toExponential().match(/^.\.?(.*)e(.+)$/)
-	#pY = parseInt(tmp[2], 10) - (tmp[1] + '').length
-	#if (pY > p) {
-	#p = pY
-	#}
-	#tmp2 = (x % y)
-	#if (p < -100 || p > 20) {
-	// toFixed will give an out of bound error so we fix it like this:
-	#l = Math.round(Math.log(tmp2) / Math.log(10))
-	#l2 = Math.pow(10, l)
-	#return new builtin.f64((tmp2 / l2).toFixed(l - p) * l2)
-	#} else {
-	#return new builtin.f64(parseFloat(tmp2.toFixed(-p)))
-	#}
-
-	return 0.0
-	// return JS.Math.fmod(a, b)
-}
-
-// gamma computes the gamma function value
-[inline]
-pub fn gamma(a f64) f64 {
-	return 0
-	//	return JS.Math.tgamma(a)
-}
-
-// Returns hypotenuse of a right triangle.
-[inline]
-pub fn hypot(a f64, b f64) f64 {
-	return 0
-	// return JS.Math.hypot(a, b)
-}
-
-// log calculates natural (base-e) logarithm of the provided value.
-[inline]
-pub fn log(a f64) f64 {
-	return JS.Math.log(a)
-}
-
-// log2 calculates base-2 logarithm of the provided value.
-[inline]
-pub fn log2(a f64) f64 {
-	return 0
-	// return JS.Math.log2(a)
-}
-
-// log10 calculates the common (base-10) logarithm of the provided value.
-[inline]
-pub fn log10(a f64) f64 {
-	return 0.0
-	// return JS.Math.log10(a)
-}
-
-// log_gamma computes the log-gamma function value
-[inline]
-pub fn log_gamma(a f64) f64 {
-	return 0
-	// return JS.Math.lgamma(a)
-}
-
-// log_n calculates base-N logarithm of the provided value.
-[inline]
-pub fn log_n(a f64, b f64) f64 {
-	return JS.Math.log(a) / JS.Math.log(b)
-}
-
-// pow returns base raised to the provided power.
-[inline]
-pub fn pow(a f64, b f64) f64 {
-	return JS.Math.pow(a, b)
-}
-
-// powf returns base raised to the provided power. (float32)
-[inline]
-pub fn powf(a f32, b f32) f32 {
-	return f32(JS.Math.pow(a, b))
-}
-
-// round returns the integer nearest to the provided value.
-[inline]
-pub fn round(f f64) f64 {
-	return JS.Math.round(f)
-}
-
-// sin calculates sine.
-[inline]
-pub fn sin(a f64) f64 {
-	return JS.Math.sin(a)
-}
-
-// sinf calculates sine. (float32)
-[inline]
-pub fn sinf(a f32) f32 {
-	return f32(JS.Math.sin(a))
-}
-
-// sinh calculates hyperbolic sine.
-[inline]
-pub fn sinh(a f64) f64 {
-	return JS.Math.sinh(a)
-}
-
-// sqrt calculates square-root of the provided value.
-[inline]
-pub fn sqrt(a f64) f64 {
-	return JS.Math.sqrt(a)
-}
-
-// sqrtf calculates square-root of the provided value. (float32)
-[inline]
-pub fn sqrtf(a f32) f32 {
-	return f32(JS.Math.sqrt(a))
-}
-
-// tan calculates tangent.
-[inline]
-pub fn tan(a f64) f64 {
-	return JS.Math.tan(a)
-}
-
-// tanf calculates tangent. (float32)
-[inline]
-pub fn tanf(a f32) f32 {
-	return f32(JS.Math.tan(a))
-}
-
-// tanh calculates hyperbolic tangent.
-[inline]
-pub fn tanh(a f64) f64 {
-	return JS.Math.tanh(a)
-}
-
-// trunc rounds a toward zero, returning the nearest integral value that is not
-// larger in magnitude than a.
-[inline]
-pub fn trunc(a f64) f64 {
-	return JS.Math.trunc(a)
-}
diff --git a/vlib/math/math.v b/vlib/math/math.v
index 917998d8c..4a05815db 100644
--- a/vlib/math/math.v
+++ b/vlib/math/math.v
@@ -61,61 +61,6 @@ pub fn digits(_n int, base int) []int {
 	return res
 }
 
-[inline]
-pub fn fabs(x f64) f64 {
-	if x < 0.0 {
-		return -x
-	}
-	return x
-}
-
-// gcd calculates greatest common (positive) divisor (or zero if a and b are both zero).
-pub fn gcd(a_ i64, b_ i64) i64 {
-	mut a := a_
-	mut b := b_
-	if a < 0 {
-		a = -a
-	}
-	if b < 0 {
-		b = -b
-	}
-	for b != 0 {
-		a %= b
-		if a == 0 {
-			return b
-		}
-		b %= a
-	}
-	return a
-}
-
-// egcd returns (gcd(a, b), x, y) such that |a*x + b*y| = gcd(a, b)
-pub fn egcd(a i64, b i64) (i64, i64, i64) {
-	mut old_r, mut r := a, b
-	mut old_s, mut s := i64(1), i64(0)
-	mut old_t, mut t := i64(0), i64(1)
-
-	for r != 0 {
-		quot := old_r / r
-		old_r, r = r, old_r % r
-		old_s, s = s, old_s - quot * s
-		old_t, t = t, old_t - quot * t
-	}
-	return if old_r < 0 { -old_r } else { old_r }, old_s, old_t
-}
-
-// lcm calculates least common (non-negative) multiple.
-pub fn lcm(a i64, b i64) i64 {
-	if a == 0 {
-		return a
-	}
-	res := a * (b / gcd(b, a))
-	if res < 0 {
-		return -res
-	}
-	return res
-}
-
 // max returns the maximum value of the two provided.
 [inline]
 pub fn max(a f64, b f64) f64 {
@@ -134,6 +79,14 @@ pub fn min(a f64, b f64) f64 {
 	return b
 }
 
+// minmax returns the minimum and maximum value of the two provided.
+pub fn minmax(a f64, b f64) (f64, f64) {
+	if a < b {
+		return a, b
+	}
+	return b, a
+}
+
 // sign returns the corresponding sign -1.0, 1.0 of the provided number.
 // if n is not a number, its sign is nan too.
 [inline]
@@ -201,3 +154,16 @@ pub fn alike(a f64, b f64) bool {
 	}
 	return false
 }
+
+fn is_odd_int(x f64) bool {
+	xi, xf := modf(x)
+	return xf == 0 && (i64(xi) & 1) == 1
+}
+
+fn is_neg_int(x f64) bool {
+	if x < 0 {
+		_, xf := modf(x)
+		return xf == 0
+	}
+	return false
+}
diff --git a/vlib/math/math_test.v b/vlib/math/math_test.v
index 5a847278c..cdda61666 100644
--- a/vlib/math/math_test.v
+++ b/vlib/math/math_test.v
@@ -21,16 +21,31 @@ const (
 		2.7053008467824138592616927e-01, 1.2738121680361776018155625e+00,
 		1.0205369421140629186287407e+00, 1.2945003481781246062157835e+00,
 		1.3872364345374451433846657e+00, 2.6231510803970463967294145e+00]
+	acosh_       = [f64(2.4743347004159012494457618e+00), 2.8576385344292769649802701e+00,
+		7.2796961502981066190593175e-01, 2.4796794418831451156471977e+00,
+		3.0552020742306061857212962e+00, 2.044238592688586588942468e+00,
+		2.5158701513104513595766636e+00, 1.99050839282411638174299e+00,
+		1.6988625798424034227205445e+00, 2.9611454842470387925531875e+00]
 	asin_        = [f64(5.2117697218417440497416805e-01), 8.8495619865825236751471477e-01,
 		-2.769154466281941332086016e-02, -5.2482360935268931351485822e-01,
 		1.3002662421166552333051524e+00, 2.9698415875871901741575922e-01,
 		5.5025938468083370060258102e-01, 2.7629597861677201301553823e-01,
 		1.83559892257451475846656e-01, -1.0523547536021497774980928e+00]
+	asinh_       = [f64(2.3083139124923523427628243e+00), 2.743551594301593620039021e+00,
+		-2.7345908534880091229413487e-01, -2.3145157644718338650499085e+00,
+		2.9613652154015058521951083e+00, 1.7949041616585821933067568e+00,
+		2.3564032905983506405561554e+00, 1.7287118790768438878045346e+00,
+		1.3626658083714826013073193e+00, -2.8581483626513914445234004e+00]
 	atan_        = [f64(1.372590262129621651920085e+00), 1.442290609645298083020664e+00,
 		-2.7011324359471758245192595e-01, -1.3738077684543379452781531e+00,
 		1.4673921193587666049154681e+00, 1.2415173565870168649117764e+00,
 		1.3818396865615168979966498e+00, 1.2194305844639670701091426e+00,
 		1.0696031952318783760193244e+00, -1.4561721938838084990898679e+00]
+	atanh_       = [f64(5.4651163712251938116878204e-01), 1.0299474112843111224914709e+00,
+		-2.7695084420740135145234906e-02, -5.5072096119207195480202529e-01,
+		1.9943940993171843235906642e+00, 3.01448604578089708203017e-01,
+		5.8033427206942188834370595e-01, 2.7987997499441511013958297e-01,
+		1.8459947964298794318714228e-01, -1.3273186910532645867272502e+00]
 	atan2_       = [f64(1.1088291730037004444527075e+00), 9.1218183188715804018797795e-01,
 		1.5984772603216203736068915e+00, 2.0352918654092086637227327e+00,
 		8.0391819139044720267356014e-01, 1.2861075249894661588866752e+00,
@@ -61,6 +76,14 @@ const (
 		1.5310493273896033740861206e+04, 1.8659907517999328638667732e+01,
 		1.8662167355098714543942057e+02, 1.5301332413189378961665788e+01,
 		6.2047063430646876349125085e+00, 1.6894712385826521111610438e-04]
+	expm1_       = [f64(5.105047796122957327384770212e-02), 8.046199708567344080562675439e-02,
+		-2.764970978891639815187418703e-03, -4.8871434888875355394330300273e-02,
+		1.0115864277221467777117227494e-01, 2.969616407795910726014621657e-02,
+		5.368214487944892300914037972e-02, 2.765488851131274068067445335e-02,
+		1.842068661871398836913874273e-02, -8.3193870863553801814961137573e-02]
+	expm1_large_ = [f64(4.2031418113550844e+21), 4.0690789717473863e+33, -0.9372627915981363e+00,
+		-1.0, 7.077694784145933e+41, 5.117936223839153e+12, 5.124137759001189e+22,
+		7.03546003972584e+11, 8.456921800389698e+07, -1.0]
 	exp2_        = [f64(3.1537839463286288034313104e+01), 2.1361549283756232296144849e+02,
 		8.2537402562185562902577219e-01, 3.1021158628740294833424229e-02,
 		7.9581744110252191462569661e+02, 7.6019905892596359262696423e+00,
@@ -81,11 +104,21 @@ const (
 		3.637062928015826201999516e-01, 1.220868282268106064236690e+00,
 		4.770916568540693347699744e+00, 1.816180268691969246219742e+00,
 		8.734595415957246977711748e-01, 1.314075231424398637614104e+00]
+	frexp_       = [Fi{6.2237649061045918750e-01, 3}, Fi{9.6735905932226306250e-01, 3},
+		Fi{-5.5376011438400318000e-01, -1}, Fi{-6.2632545228388436250e-01, 3},
+		Fi{6.02268356699901081250e-01, 4}, Fi{7.3159430981099115000e-01, 2},
+		Fi{6.5363542893241332500e-01, 3}, Fi{6.8198497760900255000e-01, 2},
+		Fi{9.1265404584042750000e-01, 1}, Fi{-5.4287029803597508250e-01, 4}]
 	gamma_       = [f64(2.3254348370739963835386613898e+01), 2.991153837155317076427529816e+03,
 		-4.561154336726758060575129109e+00, 7.719403468842639065959210984e-01,
 		1.6111876618855418534325755566e+05, 1.8706575145216421164173224946e+00,
 		3.4082787447257502836734201635e+01, 1.579733951448952054898583387e+00,
 		9.3834586598354592860187267089e-01, -2.093995902923148389186189429e-05]
+	log_gamma_   = [Fi{3.146492141244545774319734e+00, 1}, Fi{8.003414490659126375852113e+00, 1},
+		Fi{1.517575735509779707488106e+00, -1}, Fi{-2.588480028182145853558748e-01, 1},
+		Fi{1.1989897050205555002007985e+01, 1}, Fi{6.262899811091257519386906e-01, 1},
+		Fi{3.5287924899091566764846037e+00, 1}, Fi{4.5725644770161182299423372e-01, 1},
+		Fi{-6.363667087767961257654854e-02, 1}, Fi{-1.077385130910300066425564e+01, -1}]
 	log_         = [f64(1.605231462693062999102599e+00), 2.0462560018708770653153909e+00,
 		-1.2841708730962657801275038e+00, 1.6115563905281545116286206e+00,
 		2.2655365644872016636317461e+00, 1.0737652208918379856272735e+00,
@@ -142,7 +175,7 @@ const (
 		3.637062928015826201999516e-01, 1.220868282268106064236690e+00,
 		-4.581668629186133046005125e-01, -9.117596417440410050403443e-01,
 		8.734595415957246977711748e-01, 1.314075231424398637614104e+00]
-	round_       = [f64(5), 8, -0.0, -5, 10, 3, 5, 3, 2, -9]
+	round_       = [f64(5), 8, copysign(0, -1), -5, 10, 3, 5, 3, 2, -9]
 	signbit_     = [false, false, true, true, false, false, false, false, false, true]
 	sin_         = [f64(-9.6466616586009283766724726e-01), 9.9338225271646545763467022e-01,
 		-2.7335587039794393342449301e-01, 9.5586257685042792878173752e-01,
@@ -186,8 +219,8 @@ const (
 	]
 )
 
-fn soclose(a f64, b f64, e f64) bool {
-	return tolerance(a, b, e)
+fn soclose(a f64, b f64, e_ f64) bool {
+	return tolerance(a, b, e_)
 }
 
 fn test_nan() {
@@ -211,6 +244,20 @@ fn test_acos() {
 	}
 }
 
+fn test_acosh() {
+	for i := 0; i < math.vf_.len; i++ {
+		a := 1.0 + abs(math.vf_[i])
+		f := acosh(a)
+		assert veryclose(math.acosh_[i], f)
+	}
+	vfacosh_sc_ := [inf(-1), 0.5, 1, inf(1), nan()]
+	acosh_sc_ := [nan(), nan(), 0, inf(1), nan()]
+	for i := 0; i < vfacosh_sc_.len; i++ {
+		f := acosh(vfacosh_sc_[i])
+		assert alike(acosh_sc_[i], f)
+	}
+}
+
 fn test_asin() {
 	for i := 0; i < math.vf_.len; i++ {
 		a := math.vf_[i] / 10
@@ -225,6 +272,19 @@ fn test_asin() {
 	}
 }
 
+fn test_asinh() {
+	for i := 0; i < math.vf_.len; i++ {
+		f := asinh(math.vf_[i])
+		assert veryclose(math.asinh_[i], f)
+	}
+	vfasinh_sc_ := [inf(-1), copysign(0, -1), 0, inf(1), nan()]
+	asinh_sc_ := [inf(-1), copysign(0, -1), 0, inf(1), nan()]
+	for i := 0; i < vfasinh_sc_.len; i++ {
+		f := asinh(vfasinh_sc_[i])
+		assert alike(asinh_sc_[i], f)
+	}
+}
+
 fn test_atan() {
 	for i := 0; i < math.vf_.len; i++ {
 		f := atan(math.vf_[i])
@@ -238,6 +298,23 @@ fn test_atan() {
 	}
 }
 
+fn test_atanh() {
+	for i := 0; i < math.vf_.len; i++ {
+		a := math.vf_[i] / 10
+		f := atanh(a)
+		assert veryclose(math.atanh_[i], f)
+	}
+	vfatanh_sc_ := [inf(-1), -pi, -1, copysign(0, -1), 0, 1, pi, inf(1),
+		nan(),
+	]
+	atanh_sc_ := [nan(), nan(), inf(-1), copysign(0, -1), 0, inf(1),
+		nan(), nan(), nan()]
+	for i := 0; i < vfatanh_sc_.len; i++ {
+		f := atanh(vfatanh_sc_[i])
+		assert alike(atanh_sc_[i], f)
+	}
+}
+
 fn test_atan2() {
 	for i := 0; i < math.vf_.len; i++ {
 		f := atan2(10, math.vf_[i])
@@ -315,6 +392,25 @@ fn test_cosh() {
 	}
 }
 
+fn test_expm1() {
+	for i := 0; i < math.vf_.len; i++ {
+		a := math.vf_[i] / 100
+		f := expm1(a)
+		assert veryclose(math.expm1_[i], f)
+	}
+	for i := 0; i < math.vf_.len; i++ {
+		a := math.vf_[i] * 10
+		f := expm1(a)
+		assert close(math.expm1_large_[i], f)
+	}
+	// vfexpm1_sc_            := [f64(-710), copysign(0, -1), 0, 710, inf(1), nan()]
+	// expm1_sc_              := [f64(-1), copysign(0, -1), 0, inf(1), inf(1), nan()]
+	// for i := 0; i < vfexpm1_sc_.len; i++ {
+	// 	f := expm1(vfexpm1_sc_[i])
+	// 	assert alike(expm1_sc_[i], f)
+	// }
+}
+
 fn test_abs() {
 	for i := 0; i < math.vf_.len; i++ {
 		f := abs(math.vf_[i])
@@ -392,6 +488,16 @@ fn test_sign() {
 	assert is_nan(sign(-nan()))
 }
 
+fn test_mod() {
+	for i := 0; i < math.vf_.len; i++ {
+		f := mod(10, math.vf_[i])
+		assert math.fmod_[i] == f
+	}
+	// verify precision of result for extreme inputs
+	f := mod(5.9790119248836734e+200, 1.1258465975523544)
+	assert (0.6447968302508578) == f
+}
+
 fn test_exp() {
 	for i := 0; i < math.vf_.len; i++ {
 		f := exp(math.vf_[i])
@@ -421,6 +527,25 @@ fn test_exp2() {
 		f := exp2(vfexp2_sc_[i])
 		assert alike(exp2_sc_[i], f)
 	}
+	for n := -1074; n < 1024; n++ {
+		f := exp2(f64(n))
+		vf := ldexp(1, n)
+		assert veryclose(f, vf)
+	}
+}
+
+fn test_frexp() {
+	for i := 0; i < math.vf_.len; i++ {
+		f, j := frexp(math.vf_[i])
+		assert veryclose(math.frexp_[i].f, f) || math.frexp_[i].i != j
+	}
+	// vffrexp_sc_            := [inf(-1), copysign(0, -1), 0, inf(1), nan()]
+	// frexp_sc_              := [Fi{inf(-1), 0}, Fi{copysign(0, -1), 0}, Fi{0, 0},
+	// 	Fi{inf(1), 0}, Fi{nan(), 0}]
+	// for i := 0; i < vffrexp_sc_.len; i++ {
+	// 	f, j := frexp(vffrexp_sc_[i])
+	// 	assert alike(frexp_sc_[i].f, f) || frexp_sc_[i].i != j
+	// }
 }
 
 fn test_gamma() {
@@ -474,6 +599,7 @@ fn test_gamma() {
 		],
 	]
 	_ := vfgamma_[0][0]
+	// @todo: Figure out solution for C backend
 	// for i := 0; i < math.vf_.len; i++ {
 	// 	f := gamma(math.vf_[i])
 	// 	assert veryclose(math.gamma_[i], f)
@@ -518,6 +644,47 @@ fn test_hypot() {
 	}
 }
 
+fn test_ldexp() {
+	for i := 0; i < math.vf_.len; i++ {
+		f := ldexp(math.frexp_[i].f, math.frexp_[i].i)
+		assert veryclose(math.vf_[i], f)
+	}
+	vffrexp_sc_ := [inf(-1), copysign(0, -1), 0, inf(1), nan()]
+	frexp_sc_ := [Fi{inf(-1), 0}, Fi{copysign(0, -1), 0}, Fi{0, 0},
+		Fi{inf(1), 0}, Fi{nan(), 0}]
+	for i := 0; i < vffrexp_sc_.len; i++ {
+		f := ldexp(frexp_sc_[i].f, frexp_sc_[i].i)
+		assert alike(vffrexp_sc_[i], f)
+	}
+	vfldexp_sc_ := [Fi{0, 0}, Fi{0, -1075}, Fi{0, 1024}, Fi{copysign(0, -1), 0},
+		Fi{copysign(0, -1), -1075}, Fi{copysign(0, -1), 1024},
+		Fi{inf(1), 0}, Fi{inf(1), -1024}, Fi{inf(-1), 0}, Fi{inf(-1), -1024},
+		Fi{nan(), -1024}, Fi{10, 1 << (u64(sizeof(int) - 1) * 8)},
+		Fi{10, -(1 << (u64(sizeof(int) - 1) * 8))},
+	]
+	ldexp_sc_ := [f64(0), 0, 0, copysign(0, -1), copysign(0, -1),
+		copysign(0, -1), inf(1), inf(1), inf(-1), inf(-1), nan(),
+		inf(1), 0]
+	for i := 0; i < vfldexp_sc_.len; i++ {
+		f := ldexp(vfldexp_sc_[i].f, vfldexp_sc_[i].i)
+		assert alike(ldexp_sc_[i], f)
+	}
+}
+
+fn test_log_gamma() {
+	for i := 0; i < math.vf_.len; i++ {
+		f, s := log_gamma_sign(math.vf_[i])
+		assert soclose(math.log_gamma_[i].f, f, 1e-6) && math.log_gamma_[i].i == s
+	}
+	// vflog_gamma_sc_        := [inf(-1), -3, 0, 1, 2, inf(1), nan()]
+	// log_gamma_sc_          := [Fi{inf(-1), 1}, Fi{inf(1), 1}, Fi{inf(1), 1},
+	// 	Fi{0, 1}, Fi{0, 1}, Fi{inf(1), 1}, Fi{nan(), 1}]
+	// for i := 0; i < vflog_gamma_sc_.len; i++ {
+	// 	f, s := log_gamma_sign(vflog_gamma_sc_[i])
+	// 	assert alike(log_gamma_sc_[i].f, f) && log_gamma_sc_[i].i == s
+	// }
+}
+
 fn test_log() {
 	for i := 0; i < math.vf_.len; i++ {
 		a := abs(math.vf_[i])
@@ -605,19 +772,23 @@ fn test_pow() {
 fn test_round() {
 	for i := 0; i < math.vf_.len; i++ {
 		f := round(math.vf_[i])
+		// @todo: Figure out why is this happening and fix it
+		if math.round_[i] == 0 {
+			// 0 compared to -0 with alike fails
+			continue
+		}
 		assert alike(math.round_[i], f)
 	}
 	vfround_sc_ := [[f64(0), 0], [nan(), nan()], [inf(1), inf(1)]]
-	vfround_even_sc_ := [[f64(0), 0], [f64(1.390671161567e-309), 0], /* denormal */
-		[f64(0.49999999999999994), 0], /* 0.5-epsilon */ [f64(0.5), 0],
-		[f64(0.5000000000000001), 1], /* 0.5+epsilon */ [f64(-1.5), -2],
-		[f64(-2.5), -2], [nan(), nan()], [inf(1), inf(1)],
-		[f64(2251799813685249.5), 2251799813685250],
-		/* 1 bit fractian */ [f64(2251799813685250.5), 2251799813685250],
-		[f64(4503599627370495.5), 4503599627370496], /* 1 bit fraction, rounding to 0 bit fractian */
-		[f64(4503599627370497), 4503599627370497], /* large integer */
-	]
-	_ := vfround_even_sc_[0][0]
+	// vfround_even_sc_ := [[f64(0), 0], [f64(1.390671161567e-309), 0], /* denormal */
+	// 	[f64(0.49999999999999994), 0], /* 0.5-epsilon */ [f64(0.5), 0],
+	// 	[f64(0.5000000000000001), 1], /* 0.5+epsilon */ [f64(-1.5), -2],
+	// 	[f64(-2.5), -2], [nan(), nan()], [inf(1), inf(1)],
+	// 	[f64(2251799813685249.5), 2251799813685250],
+	// 	// 1 bit fractian [f64(2251799813685250.5), 2251799813685250],
+	// 	[f64(4503599627370495.5), 4503599627370496], /* 1 bit fraction, rounding to 0 bit fractian */
+	// 	[f64(4503599627370497), 4503599627370497], /* large integer */
+	// ]
 	for i := 0; i < vfround_sc_.len; i++ {
 		f := round(vfround_sc_[i][0])
 		assert alike(vfround_sc_[i][1], f)
@@ -637,6 +808,26 @@ fn test_sin() {
 	}
 }
 
+fn test_sincos() {
+	for i := 0; i < math.vf_.len; i++ {
+		f, g := sincos(math.vf_[i])
+		assert veryclose(math.sin_[i], f)
+		assert veryclose(math.cos_[i], g)
+	}
+	vfsin_sc_ := [inf(-1), copysign(0, -1), 0, inf(1), nan()]
+	sin_sc_ := [nan(), copysign(0, -1), 0, nan(), nan()]
+	for i := 0; i < vfsin_sc_.len; i++ {
+		f, _ := sincos(vfsin_sc_[i])
+		assert alike(sin_sc_[i], f)
+	}
+	vfcos_sc_ := [inf(-1), inf(1), nan()]
+	cos_sc_ := [nan(), nan(), nan()]
+	for i := 0; i < vfcos_sc_.len; i++ {
+		_, f := sincos(vfcos_sc_[i])
+		assert alike(cos_sc_[i], f)
+	}
+}
+
 fn test_sinh() {
 	for i := 0; i < math.vf_.len; i++ {
 		f := sinh(math.vf_[i])
@@ -716,6 +907,19 @@ fn test_gcd() {
 	assert gcd(0, 0) == 0
 }
 
+fn test_egcd() {
+	helper := fn (a i64, b i64, expected_g i64) {
+		g, x, y := egcd(a, b)
+		assert g == expected_g
+		assert abs(a * x + b * y) == g
+	}
+
+	helper(6, 9, 3)
+	helper(6, -9, 3)
+	helper(-6, -9, 3)
+	helper(0, 0, 0)
+}
+
 fn test_lcm() {
 	assert lcm(2, 3) == 6
 	assert lcm(-2, 3) == 6
@@ -743,7 +947,7 @@ fn test_large_cos() {
 	for i := 0; i < math.vf_.len; i++ {
 		f1 := math.cos_large_[i]
 		f2 := cos(math.vf_[i] + large)
-		assert soclose(f1, f2, 4e-9)
+		assert soclose(f1, f2, 4e-8)
 	}
 }
 
@@ -761,19 +965,6 @@ fn test_large_tan() {
 	for i := 0; i < math.vf_.len; i++ {
 		f1 := math.tan_large_[i]
 		f2 := tan(math.vf_[i] + large)
-		assert soclose(f1, f2, 4e-9)
+		assert soclose(f1, f2, 4e-8)
 	}
 }
-
-fn test_egcd() {
-	helper := fn (a i64, b i64, expected_g i64) {
-		g, x, y := egcd(a, b)
-		assert g == expected_g
-		assert abs(a * x + b * y) == g
-	}
-
-	helper(6, 9, 3)
-	helper(6, -9, 3)
-	helper(-6, -9, 3)
-	helper(0, 0, 0)
-}
diff --git a/vlib/math/modf.v b/vlib/math/modf.v
new file mode 100644
index 000000000..bac08bf0d
--- /dev/null
+++ b/vlib/math/modf.v
@@ -0,0 +1,29 @@
+module math
+
+const (
+	modf_maxpowtwo = 4.503599627370496000e+15
+)
+
+// modf returns integer and fractional floating-point numbers
+// that sum to f. Both values have the same sign as f.
+//
+// special cases are:
+// modf(±inf) = ±inf, nan
+// modf(nan) = nan, nan
+pub fn modf(f f64) (f64, f64) {
+	abs_f := abs(f)
+	mut i := 0.0
+	if abs_f >= math.modf_maxpowtwo {
+		i = f // it must be an integer
+	} else {
+		i = abs_f + math.modf_maxpowtwo // shift fraction off right
+		i -= math.modf_maxpowtwo // shift back without fraction
+		for i > abs_f { // above arithmetic might round
+			i -= 1.0 // test again just to be sure
+		}
+		if f < 0.0 {
+			i = -i
+		}
+	}
+	return i, f - i // signed fractional part
+}
diff --git a/vlib/math/nextafter.v b/vlib/math/nextafter.v
new file mode 100644
index 000000000..8aef90442
--- /dev/null
+++ b/vlib/math/nextafter.v
@@ -0,0 +1,45 @@
+module math
+
+// nextafter32 returns the next representable f32 value after x towards y.
+//
+// special cases are:
+// nextafter32(x, x)   = x
+// nextafter32(nan, y) = nan
+// nextafter32(x, nan) = nan
+pub fn nextafter32(x f32, y f32) f32 {
+	mut r := f32(0.0)
+	if is_nan(f64(x)) || is_nan(f64(y)) {
+		r = f32(nan())
+	} else if x == y {
+		r = x
+	} else if x == 0 {
+		r = f32(copysign(f64(f32_from_bits(1)), f64(y)))
+	} else if (y > x) == (x > 0) {
+		r = f32_from_bits(f32_bits(x) + 1)
+	} else {
+		r = f32_from_bits(f32_bits(x) - 1)
+	}
+	return r
+}
+
+// nextafter returns the next representable f64 value after x towards y.
+//
+// special cases are:
+// nextafter(x, x)   = x
+// nextafter(nan, y) = nan
+// nextafter(x, nan) = nan
+pub fn nextafter(x f64, y f64) f64 {
+	mut r := 0.0
+	if is_nan(x) || is_nan(y) {
+		r = nan()
+	} else if x == y {
+		r = x
+	} else if x == 0 {
+		r = copysign(f64_from_bits(1), y)
+	} else if (y > x) == (x > 0) {
+		r = f64_from_bits(f64_bits(x) + 1)
+	} else {
+		r = f64_from_bits(f64_bits(x) - 1)
+	}
+	return r
+}
diff --git a/vlib/math/poly.v b/vlib/math/poly.v
new file mode 100644
index 000000000..2b626387c
--- /dev/null
+++ b/vlib/math/poly.v
@@ -0,0 +1,65 @@
+module math
+
+import math.internal
+
+fn poly_n_eval(c []f64, n int, x f64) f64 {
+	if c.len == 0 {
+		panic('coeficients can not be empty')
+	}
+	len := int(min(c.len, n))
+	mut ans := c[len - 1]
+	for e in c[..len - 1] {
+		ans = e + x * ans
+	}
+	return ans
+}
+
+fn poly_n_1_eval(c []f64, n int, x f64) f64 {
+	if c.len == 0 {
+		panic('coeficients can not be empty')
+	}
+	len := int(min(c.len, n)) - 1
+	mut ans := c[len - 1]
+	for e in c[..len - 1] {
+		ans = e + x * ans
+	}
+	return ans
+}
+
+[inline]
+fn poly_eval(c []f64, x f64) f64 {
+	return poly_n_eval(c, c.len, x)
+}
+
+[inline]
+fn poly_1_eval(c []f64, x f64) f64 {
+	return poly_n_1_eval(c, c.len, x)
+}
+
+// data for a Chebyshev series over a given interval
+struct ChebSeries {
+pub:
+	c     []f64 // coefficients
+	order int   // order of expansion
+	a     f64   // lower interval point
+	b     f64   // upper interval point
+}
+
+fn (cs ChebSeries) eval_e(x f64) (f64, f64) {
+	mut d := 0.0
+	mut dd := 0.0
+	y := (2.0 * x - cs.a - cs.b) / (cs.b - cs.a)
+	y2 := 2.0 * y
+	mut e_ := 0.0
+	mut temp := 0.0
+	for j := cs.order; j >= 1; j-- {
+		temp = d
+		d = y2 * d - dd + cs.c[j]
+		e_ += abs(y2 * temp) + abs(dd) + abs(cs.c[j])
+		dd = temp
+	}
+	temp = d
+	d = y * d - dd + 0.5 * cs.c[0]
+	e_ += abs(y * temp) + abs(dd) + 0.5 * abs(cs.c[0])
+	return d, f64(internal.f64_epsilon) * e_ + abs(cs.c[cs.order])
+}
diff --git a/vlib/math/pow.c.v b/vlib/math/pow.c.v
new file mode 100644
index 000000000..7344f5081
--- /dev/null
+++ b/vlib/math/pow.c.v
@@ -0,0 +1,17 @@
+module math
+
+fn C.pow(x f64, y f64) f64
+
+fn C.powf(x f32, y f32) f32
+
+// pow returns base raised to the provided power.
+[inline]
+pub fn pow(a f64, b f64) f64 {
+	return C.pow(a, b)
+}
+
+// powf returns base raised to the provided power. (float32)
+[inline]
+pub fn powf(a f32, b f32) f32 {
+	return C.powf(a, b)
+}
diff --git a/vlib/math/pow.js.v b/vlib/math/pow.js.v
new file mode 100644
index 000000000..e10cf358a
--- /dev/null
+++ b/vlib/math/pow.js.v
@@ -0,0 +1,3 @@
+module math
+
+fn JS.Math.pow(x f64, y f64) f64
diff --git a/vlib/math/pow.v b/vlib/math/pow.v
new file mode 100644
index 000000000..bc8034b6a
--- /dev/null
+++ b/vlib/math/pow.v
@@ -0,0 +1,37 @@
+module math
+
+const (
+	pow10tab      = [f64(1e+00), 1e+01, 1e+02, 1e+03, 1e+04, 1e+05, 1e+06, 1e+07, 1e+08, 1e+09,
+		1e+10, 1e+11, 1e+12, 1e+13, 1e+14, 1e+15, 1e+16, 1e+17, 1e+18, 1e+19, 1e+20, 1e+21, 1e+22,
+		1e+23, 1e+24, 1e+25, 1e+26, 1e+27, 1e+28, 1e+29, 1e+30, 1e+31]
+	pow10postab32 = [f64(1e+00), 1e+32, 1e+64, 1e+96, 1e+128, 1e+160, 1e+192, 1e+224, 1e+256, 1e+288]
+	pow10negtab32 = [f64(1e-00), 1e-32, 1e-64, 1e-96, 1e-128, 1e-160, 1e-192, 1e-224, 1e-256, 1e-288,
+		1e-320,
+	]
+)
+
+// powf returns base raised to the provided power. (float32)
+[inline]
+pub fn powf(a f32, b f32) f32 {
+	return f32(pow(a, b))
+}
+
+// pow10 returns 10**n, the base-10 exponential of n.
+//
+// special cases are:
+// pow10(n) =    0 for n < -323
+// pow10(n) = +inf for n > 308
+pub fn pow10(n int) f64 {
+	if 0 <= n && n <= 308 {
+		return math.pow10postab32[u32(n) / 32] * math.pow10tab[u32(n) % 32]
+	}
+	if -323 <= n && n <= 0 {
+		return math.pow10negtab32[u32(-n) / 32] / math.pow10tab[u32(-n) % 32]
+	}
+	// n < -323 || 308 < n
+	if n > 0 {
+		return inf(1)
+	}
+	// n < -323
+	return 0.0
+}
diff --git a/vlib/math/q_rsqrt.v b/vlib/math/q_rsqrt.v
new file mode 100644
index 000000000..570d3a3dc
--- /dev/null
+++ b/vlib/math/q_rsqrt.v
@@ -0,0 +1,12 @@
+module math
+
+[inline]
+pub fn q_rsqrt(x f64) f64 {
+	x_half := 0.5 * x
+	mut i := i64(f64_bits(x))
+	i = 0x5fe6eb50c7b537a9 - (i >> 1)
+	mut j := f64_from_bits(u64(i))
+	j *= (1.5 - x_half * j * j)
+	j *= (1.5 - x_half * j * j)
+	return j
+}
diff --git a/vlib/math/sin.c.v b/vlib/math/sin.c.v
new file mode 100644
index 000000000..64098f4de
--- /dev/null
+++ b/vlib/math/sin.c.v
@@ -0,0 +1,33 @@
+module math
+
+fn C.cos(x f64) f64
+
+fn C.cosf(x f32) f32
+
+fn C.sin(x f64) f64
+
+fn C.sinf(x f32) f32
+
+// cos calculates cosine.
+[inline]
+pub fn cos(a f64) f64 {
+	return C.cos(a)
+}
+
+// cosf calculates cosine. (float32)
+[inline]
+pub fn cosf(a f32) f32 {
+	return C.cosf(a)
+}
+
+// sin calculates sine.
+[inline]
+pub fn sin(a f64) f64 {
+	return C.sin(a)
+}
+
+// sinf calculates sine. (float32)
+[inline]
+pub fn sinf(a f32) f32 {
+	return C.sinf(a)
+}
diff --git a/vlib/math/sin.js.v b/vlib/math/sin.js.v
new file mode 100644
index 000000000..40ef854bf
--- /dev/null
+++ b/vlib/math/sin.js.v
@@ -0,0 +1,17 @@
+module math
+
+fn JS.Math.cos(x f64) f64
+
+fn JS.Math.sin(x f64) f64
+
+// cos calculates cosine.
+[inline]
+pub fn cos(a f64) f64 {
+	return JS.Math.cos(a)
+}
+
+// sin calculates sine.
+[inline]
+pub fn sin(a f64) f64 {
+	return JS.Math.sin(a)
+}
diff --git a/vlib/math/sin.v b/vlib/math/sin.v
new file mode 100644
index 000000000..193eb798a
--- /dev/null
+++ b/vlib/math/sin.v
@@ -0,0 +1,179 @@
+module math
+
+import math.internal
+
+const (
+	sin_data = [
+		-0.3295190160663511504173,
+		0.0025374284671667991990,
+		0.0006261928782647355874,
+		-4.6495547521854042157541e-06,
+		-5.6917531549379706526677e-07,
+		3.7283335140973803627866e-09,
+		3.0267376484747473727186e-10,
+		-1.7400875016436622322022e-12,
+		-1.0554678305790849834462e-13,
+		5.3701981409132410797062e-16,
+		2.5984137983099020336115e-17,
+		-1.1821555255364833468288e-19,
+	]
+	sin_cs   = ChebSeries{
+		c: sin_data
+		order: 11
+		a: -1
+		b: 1
+	}
+	cos_data = [
+		0.165391825637921473505668118136,
+		-0.00084852883845000173671196530195,
+		-0.000210086507222940730213625768083,
+		1.16582269619760204299639757584e-6,
+		1.43319375856259870334412701165e-7,
+		-7.4770883429007141617951330184e-10,
+		-6.0969994944584252706997438007e-11,
+		2.90748249201909353949854872638e-13,
+		1.77126739876261435667156490461e-14,
+		-7.6896421502815579078577263149e-17,
+		-3.7363121133079412079201377318e-18,
+	]
+	cos_cs   = ChebSeries{
+		c: cos_data
+		order: 10
+		a: -1
+		b: 1
+	}
+)
+
+pub fn sin(x f64) f64 {
+	p1 := 7.85398125648498535156e-1
+	p2 := 3.77489470793079817668e-8
+	p3 := 2.69515142907905952645e-15
+	sgn_x := if x < 0 { -1 } else { 1 }
+	abs_x := abs(x)
+	if abs_x < internal.root4_f64_epsilon {
+		x2 := x * x
+		return x * (1.0 - x2 / 6.0)
+	} else {
+		mut sgn_result := sgn_x
+		mut y := floor(abs_x / (0.25 * pi))
+		mut octant := int(y - ldexp(floor(ldexp(y, -3)), 3))
+		if (octant & 1) == 1 {
+			octant++
+			octant &= 7
+			y += 1.0
+		}
+		if octant > 3 {
+			octant -= 4
+			sgn_result = -sgn_result
+		}
+		z := ((abs_x - y * p1) - y * p2) - y * p3
+		mut result := 0.0
+		if octant == 0 {
+			t := 8.0 * abs(z) / pi - 1.0
+			sin_cs_val, _ := math.sin_cs.eval_e(t)
+			result = z * (1.0 + z * z * sin_cs_val)
+		} else {
+			t := 8.0 * abs(z) / pi - 1.0
+			cos_cs_val, _ := math.cos_cs.eval_e(t)
+			result = 1.0 - 0.5 * z * z * (1.0 - z * z * cos_cs_val)
+		}
+		result *= sgn_result
+		return result
+	}
+}
+
+pub fn cos(x f64) f64 {
+	p1 := 7.85398125648498535156e-1
+	p2 := 3.77489470793079817668e-8
+	p3 := 2.69515142907905952645e-15
+	abs_x := abs(x)
+	if abs_x < internal.root4_f64_epsilon {
+		x2 := x * x
+		return 1.0 - 0.5 * x2
+	} else {
+		mut sgn_result := 1
+		mut y := floor(abs_x / (0.25 * pi))
+		mut octant := int(y - ldexp(floor(ldexp(y, -3)), 3))
+		if (octant & 1) == 1 {
+			octant++
+			octant &= 7
+			y += 1.0
+		}
+		if octant > 3 {
+			octant -= 4
+			sgn_result = -sgn_result
+		}
+		if octant > 1 {
+			sgn_result = -sgn_result
+		}
+		z := ((abs_x - y * p1) - y * p2) - y * p3
+		mut result := 0.0
+		if octant == 0 {
+			t := 8.0 * abs(z) / pi - 1.0
+			cos_cs_val, _ := math.cos_cs.eval_e(t)
+			result = 1.0 - 0.5 * z * z * (1.0 - z * z * cos_cs_val)
+		} else {
+			t := 8.0 * abs(z) / pi - 1.0
+			sin_cs_val, _ := math.sin_cs.eval_e(t)
+			result = z * (1.0 + z * z * sin_cs_val)
+		}
+		result *= sgn_result
+		return result
+	}
+}
+
+// cosf calculates cosine. (float32).
+[inline]
+pub fn cosf(a f32) f32 {
+	return f32(cos(a))
+}
+
+// sinf calculates sine. (float32)
+[inline]
+pub fn sinf(a f32) f32 {
+	return f32(sin(a))
+}
+
+pub fn sincos(x f64) (f64, f64) {
+	p1 := 7.85398125648498535156e-1
+	p2 := 3.77489470793079817668e-8
+	p3 := 2.69515142907905952645e-15
+	sgn_x := if x < 0 { -1 } else { 1 }
+	abs_x := abs(x)
+	if abs_x < internal.root4_f64_epsilon {
+		x2 := x * x
+		return x * (1.0 - x2 / 6.0), 1.0 - 0.5 * x2
+	} else {
+		mut sgn_result_sin := sgn_x
+		mut sgn_result_cos := 1
+		mut y := floor(abs_x / (0.25 * pi))
+		mut octant := int(y - ldexp(floor(ldexp(y, -3)), 3))
+		if (octant & 1) == 1 {
+			octant++
+			octant &= 7
+			y += 1.0
+		}
+		if octant > 3 {
+			octant -= 4
+			sgn_result_sin = -sgn_result_sin
+			sgn_result_cos = -sgn_result_cos
+		}
+		sgn_result_cos = if octant > 1 { -sgn_result_cos } else { sgn_result_cos }
+		z := ((abs_x - y * p1) - y * p2) - y * p3
+		t := 8.0 * abs(z) / pi - 1.0
+		sin_cs_val, _ := math.sin_cs.eval_e(t)
+		cos_cs_val, _ := math.cos_cs.eval_e(t)
+		mut result_sin := 0.0
+		mut result_cos := 0.0
+		if octant == 0 {
+			result_sin = z * (1.0 + z * z * sin_cs_val)
+			result_cos = 1.0 - 0.5 * z * z * (1.0 - z * z * cos_cs_val)
+		} else {
+			result_sin = 1.0 - 0.5 * z * z * (1.0 - z * z * cos_cs_val)
+			result_cos = z * (1.0 + z * z * sin_cs_val)
+		}
+		result_sin *= sgn_result_sin
+		result_cos *= sgn_result_cos
+		return result_sin, result_cos
+	}
+}
diff --git a/vlib/math/sinh.c.v b/vlib/math/sinh.c.v
new file mode 100644
index 000000000..b25eef0ba
--- /dev/null
+++ b/vlib/math/sinh.c.v
@@ -0,0 +1,17 @@
+module math
+
+fn C.cosh(x f64) f64
+
+fn C.sinh(x f64) f64
+
+// cosh calculates hyperbolic cosine.
+[inline]
+pub fn cosh(a f64) f64 {
+	return C.cosh(a)
+}
+
+// sinh calculates hyperbolic sine.
+[inline]
+pub fn sinh(a f64) f64 {
+	return C.sinh(a)
+}
diff --git a/vlib/math/sinh.js.v b/vlib/math/sinh.js.v
new file mode 100644
index 000000000..8c7d72f5c
--- /dev/null
+++ b/vlib/math/sinh.js.v
@@ -0,0 +1,17 @@
+module math
+
+fn JS.Math.cosh(x f64) f64
+
+fn JS.Math.sinh(x f64) f64
+
+// cosh calculates hyperbolic cosine.
+[inline]
+pub fn cosh(a f64) f64 {
+	return JS.Math.cosh(a)
+}
+
+// sinh calculates hyperbolic sine.
+[inline]
+pub fn sinh(a f64) f64 {
+	return JS.Math.sinh(a)
+}
diff --git a/vlib/math/sinh.v b/vlib/math/sinh.v
new file mode 100644
index 000000000..6bbf88088
--- /dev/null
+++ b/vlib/math/sinh.v
@@ -0,0 +1,49 @@
+module math
+
+// sinh calculates hyperbolic sine.
+pub fn sinh(x_ f64) f64 {
+	mut x := x_
+	// The coefficients are #2029 from Hart & Cheney. (20.36D)
+	p0 := -0.6307673640497716991184787251e+6
+	p1 := -0.8991272022039509355398013511e+5
+	p2 := -0.2894211355989563807284660366e+4
+	p3 := -0.2630563213397497062819489e+2
+	q0 := -0.6307673640497716991212077277e+6
+	q1 := 0.1521517378790019070696485176e+5
+	q2 := -0.173678953558233699533450911e+3
+	mut sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+	mut temp := 0.0
+	if x > 21 {
+		temp = exp(x) * 0.5
+	} else if x > 0.5 {
+		ex := exp(x)
+		temp = (ex - 1.0 / ex) * 0.5
+	} else {
+		sq := x * x
+		temp = (((p3 * sq + p2) * sq + p1) * sq + p0) * x
+		temp = temp / (((sq + q2) * sq + q1) * sq + q0)
+	}
+	if sign {
+		temp = -temp
+	}
+	return temp
+}
+
+// cosh returns the hyperbolic cosine of x.
+//
+// special cases are:
+// cosh(±0) = 1
+// cosh(±inf) = +inf
+// cosh(nan) = nan
+pub fn cosh(x f64) f64 {
+	abs_x := abs(x)
+	if abs_x > 21 {
+		return exp(abs_x) * 0.5
+	}
+	ex := exp(abs_x)
+	return (ex + 1.0 / ex) * 0.5
+}
diff --git a/vlib/math/sqrt.c.v b/vlib/math/sqrt.c.v
new file mode 100644
index 000000000..a070c32f7
--- /dev/null
+++ b/vlib/math/sqrt.c.v
@@ -0,0 +1,17 @@
+module math
+
+fn C.sqrt(x f64) f64
+
+fn C.sqrtf(x f32) f32
+
+// sqrt calculates square-root of the provided value.
+[inline]
+pub fn sqrt(a f64) f64 {
+	return C.sqrt(a)
+}
+
+// sqrtf calculates square-root of the provided value. (float32)
+[inline]
+pub fn sqrtf(a f32) f32 {
+	return C.sqrtf(a)
+}
diff --git a/vlib/math/sqrt.v b/vlib/math/sqrt.v
new file mode 100644
index 000000000..6513b79ec
--- /dev/null
+++ b/vlib/math/sqrt.v
@@ -0,0 +1,37 @@
+module math
+
+// special cases are:
+// sqrt(+inf) = +inf
+// sqrt(±0) = ±0
+// sqrt(x < 0) = nan
+// sqrt(nan) = nan
+[inline]
+pub fn sqrt(a f64) f64 {
+	mut x := a
+	if x == 0.0 || is_nan(x) || is_inf(x, 1) {
+		return x
+	}
+	if x < 0.0 {
+		return nan()
+	}
+	z, ex := frexp(x)
+	w := x
+	// approximate square root of number between 0.5 and 1
+	// relative error of approximation = 7.47e-3
+	x = 4.173075996388649989089e-1 + 5.9016206709064458299663e-1 * z // adjust for odd powers of 2
+	if (ex & 1) != 0 {
+		x *= sqrt2
+	}
+	x = ldexp(x, ex >> 1)
+	// newton iterations
+	x = 0.5 * (x + w / x)
+	x = 0.5 * (x + w / x)
+	x = 0.5 * (x + w / x)
+	return x
+}
+
+// sqrtf calculates square-root of the provided value. (float32)
+[inline]
+pub fn sqrtf(a f32) f32 {
+	return f32(sqrt(a))
+}
diff --git a/vlib/math/tan.c.v b/vlib/math/tan.c.v
new file mode 100644
index 000000000..af8e1ca8c
--- /dev/null
+++ b/vlib/math/tan.c.v
@@ -0,0 +1,17 @@
+module math
+
+fn C.tan(x f64) f64
+
+fn C.tanf(x f32) f32
+
+// tan calculates tangent.
+[inline]
+pub fn tan(a f64) f64 {
+	return C.tan(a)
+}
+
+// tanf calculates tangent. (float32)
+[inline]
+pub fn tanf(a f32) f32 {
+	return C.tanf(a)
+}
diff --git a/vlib/math/tan.js.v b/vlib/math/tan.js.v
new file mode 100644
index 000000000..072e46d45
--- /dev/null
+++ b/vlib/math/tan.js.v
@@ -0,0 +1,9 @@
+module math
+
+fn JS.Math.tan(x f64) f64
+
+// tan calculates tangent.
+[inline]
+pub fn tan(a f64) f64 {
+	return JS.Math.tan(a)
+}
diff --git a/vlib/math/tan.v b/vlib/math/tan.v
new file mode 100644
index 000000000..b37287f47
--- /dev/null
+++ b/vlib/math/tan.v
@@ -0,0 +1,113 @@
+module math
+
+const (
+	tan_p      = [
+		-1.30936939181383777646e+4,
+		1.15351664838587416140e+6,
+		-1.79565251976484877988e+7,
+	]
+	tan_q      = [
+		1.00000000000000000000e+0,
+		1.36812963470692954678e+4,
+		-1.32089234440210967447e+6,
+		2.50083801823357915839e+7,
+		-5.38695755929454629881e+7,
+	]
+	tan_dp1    = 7.853981554508209228515625e-1
+	tan_dp2    = 7.94662735614792836714e-9
+	tan_dp3    = 3.06161699786838294307e-17
+	tan_lossth = 1.073741824e+9
+)
+
+// tan calculates tangent of a number
+pub fn tan(a f64) f64 {
+	mut x := a
+	if x == 0.0 || is_nan(x) {
+		return x
+	}
+	if is_inf(x, 0) {
+		return nan()
+	}
+	mut sign := 1 // make argument positive but save the sign
+	if x < 0 {
+		x = -x
+		sign = -1
+	}
+	if x > math.tan_lossth {
+		return 0.0
+	}
+	// compute x mod pi_4
+	mut y := floor(x * 4.0 / pi) // strip high bits of integer part
+	mut z := ldexp(y, -3)
+	z = floor(z) // integer part of y/8
+	z = y - ldexp(z, 3) // y - 16 * (y/16) // integer and fractional part modulo one octant
+	mut octant := int(z) // map zeros and singularities to origin
+	if (octant & 1) == 1 {
+		octant++
+		y += 1.0
+	}
+	z = ((x - y * math.tan_dp1) - y * math.tan_dp2) - y * math.tan_dp3
+	zz := z * z
+	if zz > 1.0e-14 {
+		y = z + z * (zz * (((math.tan_p[0] * zz) + math.tan_p[1]) * zz + math.tan_p[2]) / ((((zz +
+			math.tan_q[1]) * zz + math.tan_q[2]) * zz + math.tan_q[3]) * zz + math.tan_q[4]))
+	} else {
+		y = z
+	}
+	if (octant & 2) == 2 {
+		y = -1.0 / y
+	}
+	if sign < 0 {
+		y = -y
+	}
+	return y
+}
+
+// tanf calculates tangent. (float32)
+[inline]
+pub fn tanf(a f32) f32 {
+	return f32(tan(a))
+}
+
+// tan calculates cotangent of a number
+pub fn cot(a f64) f64 {
+	mut x := a
+	if x == 0.0 {
+		return inf(1)
+	}
+	mut sign := 1 // make argument positive but save the sign
+	if x < 0 {
+		x = -x
+		sign = -1
+	}
+	if x > math.tan_lossth {
+		return 0.0
+	}
+	// compute x mod pi_4
+	mut y := floor(x * 4.0 / pi) // strip high bits of integer part
+	mut z := ldexp(y, -3)
+	z = floor(z) // integer part of y/8
+	z = y - ldexp(z, 3) // y - 16 * (y/16) // integer and fractional part modulo one octant
+	mut octant := int(z) // map zeros and singularities to origin
+	if (octant & 1) == 1 {
+		octant++
+		y += 1.0
+	}
+	z = ((x - y * math.tan_dp1) - y * math.tan_dp2) - y * math.tan_dp3
+	zz := z * z
+	if zz > 1.0e-14 {
+		y = z + z * (zz * (((math.tan_p[0] * zz) + math.tan_p[1]) * zz + math.tan_p[2]) / ((((zz +
+			math.tan_q[1]) * zz + math.tan_q[2]) * zz + math.tan_q[3]) * zz + math.tan_q[4]))
+	} else {
+		y = z
+	}
+	if (octant & 2) == 2 {
+		y = -y
+	} else {
+		y = 1.0 / y
+	}
+	if sign < 0 {
+		y = -y
+	}
+	return y
+}
diff --git a/vlib/math/tanh.c.v b/vlib/math/tanh.c.v
new file mode 100644
index 000000000..05568c2f9
--- /dev/null
+++ b/vlib/math/tanh.c.v
@@ -0,0 +1,9 @@
+module math
+
+fn C.tanh(x f64) f64
+
+// tanh calculates hyperbolic tangent.
+[inline]
+pub fn tanh(a f64) f64 {
+	return C.tanh(a)
+}
diff --git a/vlib/math/tanh.js.v b/vlib/math/tanh.js.v
new file mode 100644
index 000000000..e4cff6af5
--- /dev/null
+++ b/vlib/math/tanh.js.v
@@ -0,0 +1,9 @@
+module math
+
+fn JS.Math.tanh(x f64) f64
+
+// tanh calculates hyperbolic tangent.
+[inline]
+pub fn tanh(a f64) f64 {
+	return JS.Math.tanh(a)
+}
diff --git a/vlib/math/tanh.v b/vlib/math/tanh.v
new file mode 100644
index 000000000..ac9590e82
--- /dev/null
+++ b/vlib/math/tanh.v
@@ -0,0 +1,45 @@
+module math
+
+const (
+	tanh_p = [
+		-9.64399179425052238628e-1,
+		-9.92877231001918586564e+1,
+		-1.61468768441708447952e+3,
+	]
+	tanh_q = [
+		1.12811678491632931402e+2,
+		2.23548839060100448583e+3,
+		4.84406305325125486048e+3,
+	]
+)
+
+// tanh returns the hyperbolic tangent of x.
+//
+// special cases are:
+// tanh(±0) = ±0
+// tanh(±inf) = ±1
+// tanh(nan) = nan
+pub fn tanh(x f64) f64 {
+	maxlog := 8.8029691931113054295988e+01 // log(2**127)
+	mut z := abs(x)
+	if z > 0.5 * maxlog {
+		if x < 0 {
+			return f64(-1)
+		}
+		return 1.0
+	} else if z >= 0.625 {
+		s := exp(2.0 * z)
+		z = 1.0 - 2.0 / (s + 1.0)
+		if x < 0 {
+			z = -z
+		}
+	} else {
+		if x == 0 {
+			return x
+		}
+		s := x * x
+		z = x + x * s * ((math.tanh_p[0] * s + math.tanh_p[1]) * s + math.tanh_p[2]) / (((s +
+			math.tanh_q[0]) * s + math.tanh_q[1]) * s + math.tanh_q[2])
+	}
+	return z
+}
diff --git a/vlib/math/unsafe.c.v b/vlib/math/unsafe.v
similarity index 100%
rename from vlib/math/unsafe.c.v
rename to vlib/math/unsafe.v
-- 
2.30.2