diff --git a/mathics/builtin/numeric.py b/mathics/builtin/numeric.py index c35f298ff..c6ef8162b 100644 --- a/mathics/builtin/numeric.py +++ b/mathics/builtin/numeric.py @@ -42,15 +42,12 @@ PrecisionValueError, get_closest_precision, is_inexact, - is_numeric, ) from mathics.core.symbols import ( Symbol, - SymbolDivide, SymbolFalse, SymbolList, SymbolMachinePrecision, - SymbolTimes, SymbolTrue, ) from mathics.core.systemsymbols import SymbolPiecewise @@ -64,6 +61,7 @@ eval_Abs, eval_negate_number, eval_RealSign, + eval_Round, eval_Sign, eval_UnitStep, eval_UnitStep_multidimensional, @@ -748,6 +746,12 @@ class Round(Builtin): >> Round[-1.4] = -1 + Large exact numbers are rounded exactly: + >> Round[10^20 + 1/3] + = 100000000000000000000 + >> Round[Pi 10^20] + = 314159265358979323846 + Expressions other than numbers remain unevaluated: >> Round[x] = Round[x] @@ -779,21 +783,7 @@ def eval_one_arg(self, expr, evaluation: Evaluation): def eval_two_args(self, expr, k, evaluation: Evaluation): "Round[expr_, k_]" - if not is_numeric(expr, evaluation): - # We can't evaluate, so keep the symbolic representation. - return - - n = Expression(SymbolDivide, expr, k).round_to_float( - evaluation, permit_complex=True - ) - if n is None: - return - elif isinstance(n, complex): - n = round(n.real) - else: - n = round(n) - n = int(n) - return Expression(SymbolTimes, Integer(n), k) + return eval_Round(expr, k, evaluation) class Sign(SympyFunction): diff --git a/mathics/eval/numeric.py b/mathics/eval/numeric.py index a8a9b08f2..69dee5ff3 100644 --- a/mathics/eval/numeric.py +++ b/mathics/eval/numeric.py @@ -8,6 +8,7 @@ """ import math +from fractions import Fraction from typing import Callable, Optional import mpmath @@ -33,6 +34,7 @@ from mathics.core.convert.mpmath import from_mpmath from mathics.core.convert.sympy import from_sympy from mathics.core.element import BaseElement +from mathics.core.evaluation import Evaluation from mathics.core.expression import Expression from mathics.core.number import ( FP_MANTISA_BINARY_DIGITS, @@ -40,6 +42,7 @@ SpecialValueError, is_inexact, is_machine_precision, + is_numeric, is_zero, min_prec, prec as number_prec, @@ -49,6 +52,7 @@ SymbolAbs, SymbolComplexInfinity, SymbolDirectedInfinity, + SymbolDivide, SymbolExp, SymbolI, SymbolIndeterminate, @@ -272,6 +276,51 @@ def eval_RealSign(expr: BaseElement) -> Optional[Integer]: return eval_RealSign(try_add) +def eval_Round( + expr: BaseElement, k: BaseElement, evaluation: Evaluation +) -> Optional[BaseElement]: + """ + Round[expr, k]: the multiple of `k` closest to `expr`, ties to even. + Returns None when `expr` is not numeric or cannot be rounded. + """ + if not is_numeric(expr, evaluation): + # We can't evaluate, so keep the symbolic representation. + return None + + quotient = Expression(SymbolDivide, expr, k) + n = round_real_exactly(quotient.evaluate(evaluation)) + if n is None: + n = quotient.round_to_float(evaluation, permit_complex=True) + if n is None: + return None + elif isinstance(n, complex): + n = round(n.real) + else: + n = round(n) + n = int(n) + return Expression(SymbolTimes, Integer(n), k) + + +def round_real_exactly(x: BaseElement) -> Optional[int]: + """ + The integer closest to the real number `x`, ties to even, computed + without going through a machine float: a machine float has 53 bits, so + Round[10^20 + 1] would otherwise give 10^20. Returns None when `x` is + not known to be real. + """ + value = x.to_sympy() + if value is None or not value.is_number or value.is_real is not True: + return None + if value.is_Rational or value.is_Float: + # A Float is a binary fraction; Rational() gives it exactly. + return round(Fraction(*map(int, sympy.Rational(value).as_numer_denom()))) + # An exact irrational number such as Pi 10^20 is never halfway + # between two integers; SymPy finds its floor to whatever precision + # that needs. + floor = sympy.floor(value + sympy.Rational(1, 2)) + return int(floor) if floor.is_Integer else None + + def eval_Sign(expr: Number) -> Optional[BaseElement]: """ if expr is a number, return its sign. diff --git a/test/builtin/test_numeric.py b/test/builtin/test_numeric.py index 9734288b8..e2e707746 100644 --- a/test/builtin/test_numeric.py +++ b/test/builtin/test_numeric.py @@ -138,6 +138,17 @@ def test_realvalued(): ), ("Sign[1 - 4*I] == (1/17 - 4 I/17) Sqrt[17]", None, "True", None), ('Sign["20"]', None, "Sign[20]", None), + # Round must not go through a machine float (53 bits). + ("Round[10^20 + 1]", None, "100000000000000000001", None), + ("Round[-10^20 - 1]", None, "-100000000000000000001", None), + ("Round[123456789012345678901/10]", None, "12345678901234567890", None), + ("Round[(10^20 + 1)/2]", None, "50000000000000000000", None), + ("Round[Pi 10^20]", None, "314159265358979323846", None), + ("Round[Sqrt[2] 10^30]", None, "1414213562373095048801688724210", None), + ("Round[1.5`30 10^25]", None, "15000000000000000000000000", None), + ("Round[10^20 + 7, 10]", None, "100000000000000000010", None), + ("Round[1/2 + 10^-40]", None, "1", None), + ("Round[{5/2, 7/2, -5/2}]", None, "{2, 4, -2}", None), ], ) def test_numeric(str_expr, msgs, str_expected, fail_msg):