Add missing builtins

CHANGES.rst

@@ -8,6 +8,31 @@ New Builtins
 * ``CheckAbort``
 * ``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