1 | module math |
2 | |
3 | // Floating-point mod function. |
4 | // mod returns the floating-point remainder of x/y. |
5 | // The magnitude of the result is less than y and its |
6 | // sign agrees with that of x. |
7 | // |
8 | // special cases are: |
9 | // mod(±inf, y) = nan |
10 | // mod(nan, y) = nan |
11 | // mod(x, 0) = nan |
12 | // mod(x, ±inf) = x |
13 | // mod(x, nan) = nan |
14 | pub fn mod(x f64, y f64) f64 { |
15 | return fmod(x, y) |
16 | } |
17 | |
18 | // fmod returns the floating-point remainder of number / denom (rounded towards zero) |
19 | pub fn fmod(x f64, y f64) f64 { |
20 | if y == 0 || is_inf(x, 0) || is_nan(x) || is_nan(y) { |
21 | return nan() |
22 | } |
23 | abs_y := abs(y) |
24 | abs_y_fr, abs_y_exp := frexp(abs_y) |
25 | mut r := x |
26 | if x < 0 { |
27 | r = -x |
28 | } |
29 | for r >= abs_y { |
30 | rfr, mut rexp := frexp(r) |
31 | if rfr < abs_y_fr { |
32 | rexp = rexp - 1 |
33 | } |
34 | r = r - ldexp(abs_y, rexp - abs_y_exp) |
35 | } |
36 | if x < 0 { |
37 | r = -r |
38 | } |
39 | return r |
40 | } |
41 | |
42 | // gcd calculates greatest common (positive) divisor (or zero if a and b are both zero). |
43 | pub fn gcd(a_ i64, b_ i64) i64 { |
44 | mut a := a_ |
45 | mut b := b_ |
46 | if a < 0 { |
47 | a = -a |
48 | } |
49 | if b < 0 { |
50 | b = -b |
51 | } |
52 | for b != 0 { |
53 | a %= b |
54 | if a == 0 { |
55 | return b |
56 | } |
57 | b %= a |
58 | } |
59 | return a |
60 | } |
61 | |
62 | // egcd returns (gcd(a, b), x, y) such that |a*x + b*y| = gcd(a, b) |
63 | pub fn egcd(a i64, b i64) (i64, i64, i64) { |
64 | mut old_r, mut r := a, b |
65 | mut old_s, mut s := i64(1), i64(0) |
66 | mut old_t, mut t := i64(0), i64(1) |
67 | |
68 | for r != 0 { |
69 | quot := old_r / r |
70 | old_r, r = r, old_r % r |
71 | old_s, s = s, old_s - quot * s |
72 | old_t, t = t, old_t - quot * t |
73 | } |
74 | return if old_r < 0 { -old_r } else { old_r }, old_s, old_t |
75 | } |
76 | |
77 | // lcm calculates least common (non-negative) multiple. |
78 | pub fn lcm(a i64, b i64) i64 { |
79 | if a == 0 { |
80 | return a |
81 | } |
82 | res := a * (b / gcd(b, a)) |
83 | if res < 0 { |
84 | return -res |
85 | } |
86 | return res |
87 | } |