Implement `MachinePrecision` option for large numbers that fall outside of Python's builtin `float` mantissa

[rocky]

Description

N[integer, MachinePrecision] is not implemented.

Currently, we are forcing N for a large Integer argument (doesn't fit into a float mantissa) to be PrecisionReal.

How to Reproduce and Output Given

In the source code that implements the evaluation function, remove lines 94-95 of mathics.core.eval.nevaluator:

if isinstance(expr, Number):
        return expr.round(d)

Now use N with MachinePrecision specified:

In[1]:= N[265613988875874769338781322035779626829233452653394495974574961739092490901302182994384699044001, MachinePrecision]
Out[1]= 265613988875874769338781322035779626829233452653394495974574961739092490901302182994384699044001

I believe this is showing that numeric conversion to MachinePrecision has not been implemented.

Expected behavior

In[2]:= N[3^200, MachinePrecision]                                                                        
Out[2]= 2.65614 10^95

Additional context

Related to issue #1740

[mmatera]

Actually, the problem occurs only when the first argument is an Integer:

In[1]:= N[2`10^300, MachinePrecision]//Precision
Out[1]= MachinePrecision

In[2]:= N[2.^300, MachinePrecision]//Precision
Out[2]= MachinePrecision

In[3]:= N[2^300, MachinePrecision]//Precision
Out[3]= 15.

[rocky]

Actually, the problem occurs only when the first argument is an Integer:

Good - thanks for the clarification and correction. Title updated.

[mmatera]

This is the implementation of Integer.round, used in N[x_Integer, prec_]:

    def round(self, d: Optional[int] = None) -> Union["MachineReal", "PrecisionReal"]:
        """
        Produce a Real approximation of ``self`` with decimal precision ``d``.
        If ``d`` is  ``None``, and self.value fits in a float,
        returns a ``MachineReal`` number.
        Is the low-level equivalent to ``N[self, d]``.
        """
        if d is None:
            d = self.value.bit_length()
            if d <= FP_MANTISA_BINARY_DIGITS:
                return MachineReal(float(self.value))
            else:
                d = MACHINE_PRECISION_VALUE
        return PrecisionReal(sympy.Float(self.value, d))

For N[3^20, MachinePrecision] we get a MachinePrecision number as expected. But for N[3^200, MachinePrecision], the number just does not fit into a float number, and therefore, it is stored in a PrecisionReal value, with precision 15.

In WMA, a similar situation occurs for larger numbers:

 {N[3^646]//Precision, N[3^647] //Precision} 

What we should check is if FP_MANTISA_BINARY_DIGITS has the right value.
The problem is in the condition. A patch to handle this is coming as soon as I finish the tests.

[rocky]

This is the implementation of Integer.round, used in N[x_Integer, prec_]:

    def round(self, d: Optional[int] = None) -> Union["MachineReal", "PrecisionReal"]:
        """
        Produce a Real approximation of ``self`` with decimal precision ``d``.
        If ``d`` is  ``None``, and self.value fits in a float,
        returns a ``MachineReal`` number.
        Is the low-level equivalent to ``N[self, d]``.
        """
        if d is None:
            d = self.value.bit_length()
            if d <= FP_MANTISA_BINARY_DIGITS:
                return MachineReal(float(self.value))
            else:
                d = MACHINE_PRECISION_VALUE
        return PrecisionReal(sympy.Float(self.value, d))

For N[3^20, MachinePrecision] we get a MachinePrecision number as expected. But for N[3^200, MachinePrecision], the number just does not fit into a float number, and therefore, it is stored in a PrecisionReal value, with precision 15.

In WMA, a similar situation occurs for larger numbers:

 {N[3^646]//Precision, N[3^647] //Precision} 

What we should check is if FP_MANTISA_BINARY_DIGITS has the right value. The problem is in the condition. A patch to handle this is coming as soon as I finish the tests.

Again, very interesting indeed!

Google Gemini informs that here we can either use decimal.Decimal, mpmath.mpf, or numpy.float128.

A problem with numpy.float128, it says, is that it is not supported on Windows. When gmpy2 is installed, I imagine that it may be faster than decimal.Decimal.

(Side note, the Python source says it is an implementation of Mike Cowlishaw/IBM's General Decimal Arithmetic Specification. I remember him from my IBM and REXX language days.)

Type	Internal Representation	Approx. Decimal Digits
Built-in float	64-bit (IEEE 754)	15–17
numpy.float128	128-bit (if supported)	33–36
decimal.Decimal	Arbitrary / User-defined	Unlimited (Memory-bound)
mpmath.mpf	Arbitrary / User-defined	Unlimited (Memory-bound)

[mmatera]

#1743 reproduces what WMA does (up to rounding the reported Precision).

[rocky]

#1743 reproduces what WMA does (up to rounding the reported Precision).

See also #1744