Round: round exact and arbitrary-precision numbers exactly

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):

mathics/eval/numeric.py

@@ -8,6 +8,7 @@
 """
 
 import math
+from fractions import Fraction
 from typing import Callable, Optional
 
 import mpmath
@@ -33,13 +34,15 @@
 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,
     MIN_STORED_BITS,
     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.

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):