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
}