1 | module math |
2 | |
3 | // special cases are: |
4 | // sqrt(+inf) = +inf |
5 | // sqrt(±0) = ±0 |
6 | // sqrt(x < 0) = nan |
7 | // sqrt(nan) = nan |
8 | [inline] |
9 | pub fn sqrt(a f64) f64 { |
10 | mut x := a |
11 | if x == 0.0 || is_nan(x) || is_inf(x, 1) { |
12 | return x |
13 | } |
14 | if x < 0.0 { |
15 | return nan() |
16 | } |
17 | z, ex := frexp(x) |
18 | w := x |
19 | // approximate square root of number between 0.5 and 1 |
20 | // relative error of approximation = 7.47e-3 |
21 | x = 4.173075996388649989089e-1 + 5.9016206709064458299663e-1 * z // adjust for odd powers of 2 |
22 | if (ex & 1) != 0 { |
23 | x *= sqrt2 |
24 | } |
25 | x = ldexp(x, ex >> 1) |
26 | // newton iterations |
27 | x = 0.5 * (x + w / x) |
28 | x = 0.5 * (x + w / x) |
29 | x = 0.5 * (x + w / x) |
30 | return x |
31 | } |
32 | |
33 | // sqrtf calculates square-root of the provided value. (float32) |
34 | [inline] |
35 | pub fn sqrtf(a f32) f32 { |
36 | return f32(sqrt(a)) |
37 | } |
38 | |
39 | // sqrti calculates the integer square-root of the provided value. (i64) |
40 | pub fn sqrti(a i64) i64 { |
41 | mut x := a |
42 | mut q, mut r := i64(1), i64(0) |
43 | for ; q <= x; { |
44 | q <<= 2 |
45 | } |
46 | for ; q > 1; { |
47 | q >>= 2 |
48 | t := x - r - q |
49 | r >>= 1 |
50 | if t >= 0 { |
51 | x = t |
52 | r += q |
53 | } |
54 | } |
55 | return r |
56 | } |