v/vlib/math/sqrt.v

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))
}

// sqrti calculates the integer square-root of the provided value. (i64)
pub fn sqrti(a i64) i64 {
	mut x := a
	mut q, mut r := i64(1), i64(0)
	for ; q <= x; {
		q <<= 2
	}
	for ; q > 1; {
		q >>= 2
		t := x - r - q
		r >>= 1
		if t >= 0 {
			x = t
			r += q
		}
	}
	return r
}