From 1fff67db2cd15aaff9750adca8b4bc728bec49aa Mon Sep 17 00:00:00 2001 From: Noah <33094578+coolreader18@users.noreply.github.com> Date: Wed, 24 Mar 2021 11:11:41 -0500 Subject: [PATCH] Correct float rounding algorithm --- Lib/test/test_float.py | 6 ---- common/src/float_ops.rs | 64 ++++++++++++++++++++++++++++++++++++++++ vm/src/builtins/float.rs | 45 +++++----------------------- 3 files changed, 71 insertions(+), 44 deletions(-) diff --git a/Lib/test/test_float.py b/Lib/test/test_float.py index 451f0bcfc..c3a29b5f6 100644 --- a/Lib/test/test_float.py +++ b/Lib/test/test_float.py @@ -832,8 +832,6 @@ class RoundTestCase(unittest.TestCase): self.assertRaises(TypeError, round, NAN, "ceci n'est pas un integer") self.assertRaises(TypeError, round, -0.0, 1j) - # TODO: RUSTPYTHON - @unittest.expectedFailure def test_large_n(self): for n in [324, 325, 400, 2**31-1, 2**31, 2**32, 2**100]: self.assertEqual(round(123.456, n), 123.456) @@ -846,8 +844,6 @@ class RoundTestCase(unittest.TestCase): self.assertEqual(round(1e150, 309), 1e150) self.assertEqual(round(1.4e-315, 315), 1e-315) - # TODO: RUSTPYTHON - @unittest.expectedFailure def test_small_n(self): for n in [-308, -309, -400, 1-2**31, -2**31, -2**31-1, -2**100]: self.assertEqual(round(123.456, n), 0.0) @@ -855,8 +851,6 @@ class RoundTestCase(unittest.TestCase): self.assertEqual(round(1e300, n), 0.0) self.assertEqual(round(1e-320, n), 0.0) - # TODO: RUSTPYTHON - @unittest.expectedFailure def test_overflow(self): self.assertRaises(OverflowError, round, 1.6e308, -308) self.assertRaises(OverflowError, round, -1.7e308, -308) diff --git a/common/src/float_ops.rs b/common/src/float_ops.rs index 731258d9b..83c12375f 100644 --- a/common/src/float_ops.rs +++ b/common/src/float_ops.rs @@ -1,5 +1,6 @@ use num_bigint::{BigInt, ToBigInt}; use num_traits::{Float, Signed, ToPrimitive, Zero}; +use std::f64; pub fn ufrexp(value: f64) -> (f64, i32) { if 0.0 == value { @@ -366,6 +367,69 @@ pub fn ulp(x: f64) -> f64 { } } +pub fn round_float_digits(x: f64, ndigits: i32) -> Option { + 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 { + // according to cpython: pow1 and pow2 are each safe from overflow, but + // pow1*pow2 ~= pow(10.0, ndigits) might overflow + 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 { + // ndigits >= 0 + (z / pow2) / pow1 + } else { + z * pow1 + }; + + if !z.is_finite() { + // overflow + return None; + } + + z + } + }; + Some(float) +} + #[test] fn test_to_hex() { use rand::Rng; diff --git a/vm/src/builtins/float.rs b/vm/src/builtins/float.rs index 51a88db64..017e1d176 100644 --- a/vm/src/builtins/float.rs +++ b/vm/src/builtins/float.rs @@ -1,8 +1,7 @@ use num_bigint::{BigInt, ToBigInt}; use num_complex::Complex64; use num_rational::Ratio; -use num_traits::{pow, Signed, ToPrimitive, Zero}; -use std::f64; +use num_traits::{Signed, ToPrimitive, Zero}; use super::bytes::PyBytes; use super::int::{self, PyInt, PyIntRef}; @@ -395,43 +394,13 @@ impl PyFloat { let ndigits = ndigits.flatten(); let value = if let Some(ndigits) = ndigits { let ndigits = ndigits.borrow_value(); - let float = if ndigits.is_zero() { - let fract = self.value.fract(); - if (fract.abs() - 0.5).abs() < f64::EPSILON { - if self.value.trunc() % 2.0 == 0.0 { - self.value - fract - } else { - self.value + fract - } - } else { - self.value.round() - } - } else { - let ndigits = match ndigits.to_isize() { - Some(n) => n, - None if ndigits.is_positive() => isize::MAX, - None => isize::MIN, - }; - const NDIGITS_MAX: isize = ((f64::MANTISSA_DIGITS as i32 - f64::MIN_EXP) as f64 - * f64::consts::LOG10_2) as isize; - const NDIGITS_MIN: isize = - -(((f64::MAX_EXP + 1) as f64 * f64::consts::LOG10_2) as isize); - if ndigits > NDIGITS_MAX { - self.value - } else if ndigits > NDIGITS_MIN { - 0.0f64.copysign(self.value) - } else if ndigits >= 0 { - (self.value * pow(10.0, ndigits as usize)).round() / pow(10.0, ndigits as usize) - } else { - let result = (self.value / pow(10.0, (-ndigits) as usize)).round() - * pow(10.0, (-ndigits) as usize); - if result.is_nan() { - 0.0 - } else { - result - } - } + let ndigits = match ndigits.to_i32() { + Some(n) => n, + None if ndigits.is_positive() => i32::MAX, + None => i32::MIN, }; + let float = float_ops::round_float_digits(self.value, ndigits) + .ok_or_else(|| vm.new_overflow_error("overflow ocurred during round".to_owned()))?; vm.ctx.new_float(float) } else { let fract = self.value.fract();