use malachite_bigint::{BigInt, ToBigInt};
use num_traits::{Float, Signed, ToPrimitive, Zero};
use std::f64;
pub fn ufrexp(value: f64) -> (f64, i32) {
if 0.0 == value {
(0.0, 0i32)
} else {
let bits = value.to_bits();
let exponent: i32 = ((bits >> 52) & 0x7ff) as i32 - 1022;
let mantissa_bits = bits & (0x000f_ffff_ffff_ffff) | (1022 << 52);
(f64::from_bits(mantissa_bits), exponent)
}
}
pub fn eq_int(value: f64, other: &BigInt) -> bool {
if let (Some(self_int), Some(other_float)) = (value.to_bigint(), other.to_f64()) {
value == other_float && self_int == *other
} else {
false
}
}
pub fn lt_int(value: f64, other_int: &BigInt) -> bool {
match (value.to_bigint(), other_int.to_f64()) {
(Some(self_int), Some(other_float)) => value < other_float || self_int < *other_int,
(Some(_), None) => other_int.is_positive(),
_ if value.is_infinite() => value.is_sign_negative(),
_ => false,
}
}
pub fn gt_int(value: f64, other_int: &BigInt) -> bool {
match (value.to_bigint(), other_int.to_f64()) {
(Some(self_int), Some(other_float)) => value > other_float || self_int > *other_int,
(Some(_), None) => other_int.is_negative(),
_ if value.is_infinite() => value.is_sign_positive(),
_ => false,
}
}
pub fn div(v1: f64, v2: f64) -> Option<f64> {
if v2 != 0.0 {
Some(v1 / v2)
} else {
None
}
}
pub fn mod_(v1: f64, v2: f64) -> Option<f64> {
if v2 != 0.0 {
let val = v1 % v2;
match (v1.signum() as i32, v2.signum() as i32) {
(1, 1) | (-1, -1) => Some(val),
_ if (v1 == 0.0) || (v1.abs() == v2.abs()) => Some(val.copysign(v2)),
_ => Some((val + v2).copysign(v2)),
}
} else {
None
}
}
pub fn floordiv(v1: f64, v2: f64) -> Option<f64> {
if v2 != 0.0 {
Some((v1 / v2).floor())
} else {
None
}
}
pub fn divmod(v1: f64, v2: f64) -> Option<(f64, f64)> {
if v2 != 0.0 {
let mut m = v1 % v2;
let mut d = (v1 - m) / v2;
if v2.is_sign_negative() != m.is_sign_negative() {
m += v2;
d -= 1.0;
}
Some((d, m))
} else {
None
}
}
#[allow(clippy::float_cmp)]
pub fn nextafter(x: f64, y: f64) -> f64 {
if x == y {
y
} else if x.is_nan() || y.is_nan() {
f64::NAN
} else if x >= f64::INFINITY {
f64::MAX
} else if x <= f64::NEG_INFINITY {
f64::MIN
} else if x == 0.0 {
f64::from_bits(1).copysign(y)
} else {
let b = x.to_bits();
let bits = if (y > x) == (x > 0.0) { b + 1 } else { b - 1 };
let ret = f64::from_bits(bits);
if ret == 0.0 {
ret.copysign(x)
} else {
ret
}
}
}
pub fn ulp(x: f64) -> f64 {
if x.is_nan() {
return x;
}
let x = x.abs();
let x2 = nextafter(x, f64::INFINITY);
if x2.is_infinite() {
let x2 = nextafter(x, f64::NEG_INFINITY);
x - x2
} else {
x2 - x
}
}
pub fn round_float_digits(x: f64, ndigits: i32) -> Option<f64> {
let float = if ndigits.is_zero() {
let fract = x.fract();
if (fract.abs() - 0.5).abs() < f64::EPSILON {
if x.trunc() % 2.0 == 0.0 {
x - fract
} else {
x + fract
}
} else {
x.round()
}
} else {
const NDIGITS_MAX: i32 =
((f64::MANTISSA_DIGITS as i32 - f64::MIN_EXP) as f64 * f64::consts::LOG10_2) as i32;
const NDIGITS_MIN: i32 = -(((f64::MAX_EXP + 1) as f64 * f64::consts::LOG10_2) as i32);
if ndigits > NDIGITS_MAX {
x
} else if ndigits < NDIGITS_MIN {
0.0f64.copysign(x)
} else {
let (y, pow1, pow2) = if ndigits >= 0 {
let (pow1, pow2) = if ndigits > 22 {
(10.0.powf((ndigits - 22) as f64), 1e22)
} else {
(10.0.powf(ndigits as f64), 1.0)
};
let y = (x * pow1) * pow2;
if !y.is_finite() {
return Some(x);
}
(y, pow1, Some(pow2))
} else {
let pow1 = 10.0.powf((-ndigits) as f64);
(x / pow1, pow1, None)
};
let z = y.round();
#[allow(clippy::float_cmp)]
let z = if (y - z).abs() == 0.5 {
2.0 * (y / 2.0).round()
} else {
z
};
let z = if let Some(pow2) = pow2 {
(z / pow2) / pow1
} else {
z * pow1
};
if !z.is_finite() {
return None;
}
z
}
};
Some(float)
}