Add missing builtins (#1133)

I added a bunch of WMA builtins that aren't in Mathics yet:

- NumberDigit
- BellB
- EulerE
- JacobiSymbol
- KroneckerSymbol
- LucasL
- PolygonalNumber
- SquaresR
- LinearRecurrence
- ~~RootSum~~ I think this one needs some more work for converting to
sympy
- SeriesCoefficient
- DivisorSigma
- DivisorSum
- IntegerPart
- IntegerPartitions
- MersennePrimeExponent
- MoebiusMu
- PowersRepresentations
- HypergeometricU
- LambertW
- Subfactorial
- PolyLog
- ReverseSort

I also expanded the functionality of a few existing builtins to match
the behaviour of WMA:

- handle non-integer bounds in Sum
- fixed a bug in RealDigits where it would report the wrong number of
digits for powers of 10 (or whatever base is specified)
- support Fibonacci polynomials
- support symbolic bounds in Table (and other IterationFunctions)
- support providing SeriesData to CoefficientList
- added the listable attribute to a few functions
- allow supplying a precision when calling `N[]` on a complex number

I can split some of these out into separate PRs if they need more work
to merge.

I've started setting up some automated tests to pull all the Mathematica
programs listed in OEIS and check that Mathics can evaluate them to the
correct sequence values - these are all just the low hanging fruit that
fell out of doing that for the first couple of hundred sequences.

---------

Co-authored-by: Juan Mauricio Matera <matera@fisica.unlp.edu.ar>

Author: davidar

CHANGES.rst

@@ -9,6 +9,31 @@ New Builtins
 * ``FileNameDrop``
 * ``SetEnvironment``
 
+By `@davidar <https://github.com/davidar>`_:
+
+* ``BellB``
+* ``DivisorSigma``
+* ``DivisorSum``
+* ``EulerE``
+* ``HypergeometricU``
+* ``IntegerPart``
+* ``IntegerPartitions``
+* ``JacobiSymbol``
+* ``KroneckerSymbol``
+* ``LambertW``
+* ``LinearRecurrence``
+* ``LucasL``
+* ``MersennePrimeExponent``
+* ``MoebiusMu``
+* ``NumberDigit``
+* ``PolygonalNumber``
+* ``PolyLog``
+* ``PowersRepresentations``
+* ``ReverseSort``
+* ``SeriesCoefficient``
+* ``SquaresR``
+* ``Subfactorial``
+
 ``mathics`` command line
 ++++++++++++++++++++++++
 

mathics/builtin/arithmetic.py

@@ -959,6 +959,14 @@ class Sum(IterationFunction, SympyFunction):
     Verify algebraic identities:
     >> Sum[x ^ 2, {x, 1, y}] - y * (y + 1) * (2 * y + 1) / 6
      = 0
+
+    Non-integer bounds:
+    >> Sum[i, {i, 1, 2.5}]
+     = 3
+    >> Sum[i, {i, 1.1, 2.5}]
+     = 3.2
+    >> Sum[k, {k, I, I + 1.5}]
+     = 1 + 2 I
     """
 
     summary_text = "discrete sum"
@@ -1020,7 +1028,11 @@ def to_sympy(self, expr, **kwargs) -> Optional[SympyExpression]:
                 # If we have integer bounds, we'll use Mathics's iterator Sum
                 # (which is Plus)
 
-                if all(hasattr(i, "is_integer") and i.is_integer for i in bounds[1:]):
+                if all(
+                    (hasattr(i, "is_integer") and i.is_integer)
+                    or (hasattr(i, "is_finite") and i.is_finite and i.is_constant())
+                    for i in bounds[1:]
+                ):
                     # When we have integer bounds, it is better to not use Sympy but
                     # use Mathics evaluation. We turn:
                     # Sum[f[x], {<limits>}] into
@@ -1029,9 +1041,10 @@ def to_sympy(self, expr, **kwargs) -> Optional[SympyExpression]:
                     values = Expression(SymbolTable, *expr.elements).evaluate(
                         evaluation
                     )
-                    ret = self.get_result(values.elements).evaluate(evaluation)
-                    # Make sure to convert the result back to sympy.
-                    return ret.to_sympy()
+                    if values.get_head_name() != SymbolTable.get_name():
+                        ret = self.get_result(values.elements).evaluate(evaluation)
+                        # Make sure to convert the result back to sympy.
+                        return ret.to_sympy()
 
             if None not in bounds:
                 return sympy.summation(f_sympy, bounds)

mathics/builtin/atomic/numbers.py

@@ -30,7 +30,7 @@
     MachineReal,
     Rational,
 )
-from mathics.core.attributes import A_LISTABLE, A_PROTECTED
+from mathics.core.attributes import A_LISTABLE, A_PROTECTED, A_READ_PROTECTED
 from mathics.core.builtin import Builtin, Predefined
 from mathics.core.convert.python import from_python
 from mathics.core.expression import Expression
@@ -60,7 +60,9 @@
 
 @lru_cache()
 def log_n_b(py_n, py_b) -> int:
-    return int(mpmath.ceil(mpmath.log(py_n, py_b))) if py_n != 0 and py_n != 1 else 1
+    return (
+        int(mpmath.floor(mpmath.log(py_n, py_b))) + 1 if py_n != 0 and py_n != 1 else 1
+    )
 
 
 def check_finite_decimal(denominator):
@@ -376,6 +378,33 @@ def eval(self, n, b, evaluation):
         return Integer(j)
 
 
+class NumberDigit(Builtin):
+    """
+    <url>:WMA link:
+    https://reference.wolfram.com/language/ref/NumberDigit.html</url>
+
+    <dl>
+      <dt>'NumberDigit[$x$, $n$, $b$]'
+      <dd>returns the coefficient of $b^n$ in the base-$b$ representation of $x$.
+    </dl>
+
+    >> NumberDigit[123456, 2]
+     = 4
+    >> NumberDigit[12.3456, -1]
+     = 3
+
+    """
+
+    attributes = A_PROTECTED | A_READ_PROTECTED
+
+    summary_text = "digits of a real number"
+
+    rules = {
+        "NumberDigit[x_, n_Integer]": "NumberDigit[x, n, 10]",
+        "NumberDigit[x_, n_Integer, b_Integer]": "RealDigits[x, b, 1, n][[1]][[1]]",
+    }
+
+
 class RealDigits(Builtin):
     """
     <url>:WMA link:
@@ -420,6 +449,9 @@ class RealDigits(Builtin):
     Return 25 digits of in base 10:
     >> RealDigits[Pi, 10, 25]
      = {{3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3, 8, 4, 6, 2, 6, 4, 3}, 1}
+
+    >> RealDigits[10]
+     = {{1, 0}, 2}
     """
 
     attributes = A_LISTABLE | A_PROTECTED
@@ -445,7 +477,7 @@ def eval_rational_with_base(self, n, b, evaluation):
         if check_finite_decimal(n.denominator().get_int_value()) and not py_b % 2:
             return self.eval_with_base(n, b, evaluation)
         else:
-            exp = int(mpmath.ceil(mpmath.log(py_n, py_b)))
+            exp = log_n_b(py_n, py_b)
             (head, tails) = convert_repeating_decimal(
                 py_n.as_numer_denom()[0], py_n.as_numer_denom()[1], py_b
             )

mathics/builtin/intfns/recurrence.py

@@ -29,6 +29,8 @@ class Fibonacci(MPMathFunction):
     <dl>
       <dt>'Fibonacci[$n$]'
       <dd>computes the $n$th Fibonacci number.
+      <dt>'Fibonacci[$n$, $x$]'
+      <dd>computes the Fibonacci polynomial $F$_$n$($x$).
     </dl>
 
     >> Fibonacci[0]
@@ -39,6 +41,8 @@ class Fibonacci(MPMathFunction):
      = 55
     >> Fibonacci[200]
      = 280571172992510140037611932413038677189525
+    >> Fibonacci[7, x]
+     = 1 + 6 x ^ 2 + 5 x ^ 4 + x ^ 6
     """
 
     nargs = {1}
@@ -47,6 +51,11 @@ class Fibonacci(MPMathFunction):
     mpmath_name = "fibonacci"
     summary_text = "Fibonacci's numbers"
 
+    rules = {
+        "Fibonacci[0, x_]": "0",
+        "Fibonacci[n_Integer?Negative, x_]": "Fibonacci[-n, x]",
+    }
+
 
 class HarmonicNumber(MPMathFunction):
     """
@@ -73,6 +82,31 @@ class HarmonicNumber(MPMathFunction):
     sympy_name = "harmonic"
 
 
+class LinearRecurrence(Builtin):
+    """
+    <url>:WMA link:https://reference.wolfram.com/language/ref/LinearRecurrence.html</url>
+
+    <dl>
+      <dt>'LinearRecurrence[$ker$, $init$, $n$]'
+      <dd>computes $n$ terms of the linear recurrence with kernel $ker$ and intial values $init$
+    </dl>
+
+    >> LinearRecurrence[{1, 1}, {1, 1}, 10]
+     = {1, 1, 2, 3, 5, 8, 13, 21, 34, 55}
+    >> LinearRecurrence[{1, 1}, {1, 1}, {5, 5}]
+     = {5}
+    """
+
+    attributes = A_PROTECTED | A_READ_PROTECTED
+    summary_text = "linear recurrence"
+
+    rules = {
+        "LinearRecurrence[ker_List, init_List, n_Integer]": "Nest[Append[#, Reverse[ker] . Take[#, -Length[ker]]] &, init, n - Length[init]]",
+        "LinearRecurrence[ker_List, init_List, {n_Integer?Positive}]": "LinearRecurrence[ker, init, n][[n]]",
+        "LinearRecurrence[ker_List, init_List, {nmin_Integer?Positive, nmax_Integer?Positive}]": "LinearRecurrence[ker, init, nmax][[nmin;;nmax]]",
+    }
+
+
 # Note: WL allows StirlingS1[{2, 4, 6}, 2], but we don't (yet).
 class StirlingS1(Builtin):
     """

mathics/builtin/list/constructing.py

@@ -522,6 +522,7 @@ class Table(IterationFunction):
      = {x, x, x}
     >> n = 0; Table[n = n + 1, {5}]
      = {1, 2, 3, 4, 5}
+    #> Clear[n]
     >> Table[i, {i, 4}]
      = {1, 2, 3, 4}
     >> Table[i, {i, 2, 5}]
@@ -536,6 +537,14 @@ class Table(IterationFunction):
     'Table' supports multi-dimensional tables:
     >> Table[{i, j}, {i, {a, b}}, {j, 1, 2}]
      = {{{a, 1}, {a, 2}}, {{b, 1}, {b, 2}}}
+
+    Symbolic bounds:
+    >> Table[x, {x, a, a + 5 n, n}]
+     = {a, a + n, a + 2 n, a + 3 n, a + 4 n, a + 5 n}
+
+    The lower bound is always included even for large step sizes:
+    >> Table[i, {i, 1, 9, Infinity}]
+     = {1}
     """
 
     rules = {

mathics/builtin/numbers/algebra.py

@@ -901,6 +901,8 @@ class CoefficientList(Builtin):
      = {{5, d, 0, b}, {c, 0, 0, 0}, {a, 0, 0, 0}}
     >> CoefficientList[(x - 2 y + 3 z)^3, {x, y, z}]
      = {{{0, 0, 0, 27}, {0, 0, -54, 0}, {0, 36, 0, 0}, {-8, 0, 0, 0}}, {{0, 0, 27, 0}, {0, -36, 0, 0}, {12, 0, 0, 0}, {0, 0, 0, 0}}, {{0, 9, 0, 0}, {-6, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}, {{1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}
+    >> CoefficientList[Series[Log[1-x], {x, 0, 9}], x]
+     = {0, -1, -1 / 2, -1 / 3, -1 / 4, -1 / 5, -1 / 6, -1 / 7, -1 / 8, -1 / 9}
     """
 
     messages = {
@@ -914,7 +916,7 @@ def eval_noform(self, expr, evaluation):
         "CoefficientList[expr_]"
         evaluation.message("CoefficientList", "argtu")
 
-    def eval(self, expr, form, evaluation):
+    def eval(self, expr: Expression, form: Expression, evaluation: Evaluation):
         "CoefficientList[expr_, form_]"
         vars = [form] if not form.has_form("List", None) else [v for v in form.elements]
 
@@ -937,6 +939,18 @@ def eval(self, expr, form, evaluation):
             return ListExpression(expr)
         elif form.has_form("List", 0):
             return expr
+        elif expr.get_head_name() == "System`SeriesData":
+            coeffs: ListExpression
+            nmin: Integer
+            nmax: Integer
+            x, x0, coeffs, nmin, nmax, den = expr.elements
+            if x == form and x0 == Integer0 and den == Integer1:
+                return ListExpression(
+                    *[
+                        coeffs.elements[i - nmin.value] if i >= nmin.value else Integer0
+                        for i in range(0, nmax.value)
+                    ]
+                )
 
         sympy_expr = expr.to_sympy()
         sympy_vars = [v.to_sympy() for v in vars]
@@ -953,8 +967,7 @@ def eval(self, expr, form, evaluation):
 
             # single & multiple variables cases
             if not form.has_form("List", None):
-                return Expression(
-                    SymbolList,
+                return ListExpression(
                     *[
                         _coefficient(
                             self.__class__.__name__, expr, form, Integer(n), evaluation
@@ -964,8 +977,7 @@ def eval(self, expr, form, evaluation):
                 )
             elif form.has_form("List", 1):
                 form = form.elements[0]
-                return Expression(
-                    SymbolList,
+                return ListExpression(
                     *[
                         _coefficient(
                             self.__class__.__name__, expr, form, Integer(n), evaluation

mathics/builtin/numbers/calculus.py

@@ -63,6 +63,7 @@
     SymbolConditionalExpression,
     SymbolD,
     SymbolDerivative,
+    SymbolIndeterminate,
     SymbolInfinity,
     SymbolInfix,
     SymbolIntegrate,
@@ -1731,6 +1732,50 @@ def eval_multivariate_series(self, f, varspec, evaluation: Evaluation):
         return None
 
 
+class SeriesCoefficient(Builtin):
+    """
+    <url>:WMA link:https://reference.wolfram.com/language/ref/SeriesCoefficient.html</url>
+
+    <dl>
+      <dt>'SeriesCoefficient[$series$, $n$]'
+      <dd>Find the $n$th coefficient in the given $series$
+    </dl>
+
+    >> SeriesCoefficient[Series[Exp[Sin[x]], {x, 0, 10}], 8]
+     = 31 / 5760
+    >> SeriesCoefficient[Exp[-x], {x, 0, 5}]
+     = -1 / 120
+
+    >> SeriesCoefficient[SeriesData[x, c, Table[i^2, {i, 10}], 7, 17, 3], 14/3]
+     = 64
+    >> SeriesCoefficient[SeriesData[x, c, Table[i^2, {i, 10}], 7, 17, 3], 6/3]
+     = 0
+    >> SeriesCoefficient[SeriesData[x, c, Table[i^2, {i, 10}], 7, 17, 3], 17/3]
+     = Indeterminate
+    """
+
+    attributes = A_PROTECTED
+    summary_text = "power series coefficient"
+
+    rules = {
+        "SeriesCoefficient[f_, {x_Symbol, x0_, n_Integer}]": "SeriesCoefficient[Series[f, {x, x0, n}], n]"
+    }
+
+    def eval(self, series: Expression, n: Rational, evaluation: Evaluation):
+        """SeriesCoefficient[series_SeriesData, n_]"""
+        coeffs: ListExpression
+        nmin: Integer
+        nmax: Integer
+        den: Integer
+        coeffs, nmin, nmax, den = series.elements[2:]
+        index = n.value * den.value - nmin.value
+        if index < 0:
+            return Integer0
+        if index >= nmax.value - nmin.value:
+            return SymbolIndeterminate
+        return coeffs[index]
+
+
 class SeriesData(Builtin):
     """
 

mathics/builtin/specialfns/bessel.py

@@ -558,6 +558,28 @@ class HankelH2(_Bessel):
     sympy_name = "hankel2"
 
 
+class HypergeometricU(MPMathFunction):
+    """
+    <url>
+    :Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
+    :mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html</url>)
+    <dl>
+      <dt>'HypergeometricU[$a$, $b$, $z$]'
+      <dd>returns $U$($a$, $b$, $z$).
+    </dl>
+
+    >> HypergeometricU[3, 2, 1.]
+     = 0.105479
+    """
+
+    attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
+    summary_text = "Tricomi confluent hypergeometric function"
+    mpmath_name = "hyperu"
+    sympy_name = ""
+    nargs = {3}
+
+
 # Kelvin Functions