Add missing builtins: combinatorial functions

mathics/builtin/intfns/combinatorial.py

@@ -17,12 +17,14 @@
 from mathics.core.atoms import Integer
 from mathics.core.attributes import (
     A_LISTABLE,
+    A_N_HOLD_FIRST,
     A_NUMERIC_FUNCTION,
     A_ORDERLESS,
     A_PROTECTED,
     A_READ_PROTECTED,
 )
 from mathics.core.builtin import Builtin, MPMathFunction, SympyFunction
+from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
 from mathics.core.symbols import (
@@ -39,6 +41,36 @@
 SymbolSubsets = Symbol("Subsets")
 
 
+class BellB(SympyFunction):
+    """
+    <url>
+    :Bell number: https://en.wikipedia.org/wiki/Bell_number</url> (<url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/combinatorial.html#sympy.functions.combinatorial.numbers.bell</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/BellB.html</url>)
+    <dl>
+      <dt>'BellB[$n$]'
+      <dd>Bell number $B$_$n$.
+
+      <dt>'BellB[$n$, $x$]'
+      <dd>Bell polynomial $B$_$n$($x$).
+    </dl>
+
+    >> BellB[10]
+     = 115975
+
+    >> BellB[5, x]
+     = x + 15 x ^ 2 + 25 x ^ 3 + 10 x ^ 4 + x ^ 5
+    """
+
+    attributes = A_LISTABLE | A_N_HOLD_FIRST | A_PROTECTED | A_READ_PROTECTED
+    summary_text = "Bell numbers"
+    sympy_name = "bell"
+
+    def eval(self, z, evaluation: Evaluation):
+        "%(name)s[z__]"
+        return super().eval(z, evaluation)
+
+
 class _BooleanDissimilarity(Builtin):
     @staticmethod
     def _to_bool_vector(u):
@@ -184,6 +216,36 @@ def _compute(self, n, c_ff, c_ft, c_tf, c_tt):
         )
 
 
+class EulerE(SympyFunction):
+    """
+    <url>
+    :Euler numbers: https://en.wikipedia.org/wiki/Euler_numbers</url> (<url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/combinatorial.html#sympy.functions.combinatorial.numbers.euler</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/EulerE.html</url>)
+    <dl>
+      <dt>'EulerE[$n$]'
+      <dd>Euler number $E$_$n$.
+
+      <dt>'EulerE[$n$, $x$]'
+      <dd>Euler polynomial $E$_$n$($x$).
+    </dl>
+
+    >> Table[EulerE[k], {k, 0, 10}]
+     = {1, 0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521}
+
+    >> EulerE[5, z]
+     = -1 / 2 + 5 z ^ 2 / 2 - 5 z ^ 4 / 2 + z ^ 5
+    """
+
+    attributes = A_LISTABLE | A_PROTECTED
+    summary_text = "Euler numbers"
+    sympy_name = "euler"
+
+    def eval(self, z, evaluation: Evaluation):
+        "%(name)s[z__]"
+        return super().eval(z, evaluation)
+
+
 class JaccardDissimilarity(_BooleanDissimilarity):
     """
     <url>
@@ -214,6 +276,44 @@ def _compute(self, n, c_ff, c_ft, c_tf, c_tt):
         )
 
 
+class LucasL(SympyFunction):
+    """
+    <url>
+    :Lucas Number:
+    https://en.wikipedia.org/wiki/Lucas_number</url> (<url>
+    :SymPy:
+    https://docs.sympy.org/latest/modules/functions/combinatorial.html#sympy.functions.combinatorial.numbers.lucas</url>, \
+    <url>
+    :WMA:
+    https://reference.wolfram.com/language/ref/LucasL.html</url>)
+
+    <dl>
+      <dt>'LucasL[$n$]'
+      <dd>gives the $n$th Lucas number.
+    </dl>
+
+    A list of the first five Lucas numbers:
+    >> Table[LucasL[n], {n, 1, 5}]
+     = {1, 3, 4, 7, 11}
+    >> Series[LucasL[1/2, x], {x, 0, 5}]
+     = 1 + 1 / 4 x + 1 / 32 x ^ 2 + (-1 / 128) x ^ 3 + (-5 / 2048) x ^ 4 + 7 / 8192 x ^ 5 + O[x] ^ 6
+    """
+
+    attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
+
+    summary_text = "lucas number"
+    sympy_name = "lucas"
+
+    rules = {
+        "LucasL[n_, 1]": "LucasL[n]",
+        "LucasL[n_, x_]": "(x/2 + Sqrt[1 + x^2 / 4])^n + Cos[n Pi] / (x/2 + Sqrt[1 + x^2 / 4])^n // Simplify",
+    }
+
+    def eval_integer(self, n: Integer, evaluation):
+        "LucasL[n_Integer]"
+        return self.eval(n, evaluation)
+
+
 class MatchingDissimilarity(_BooleanDissimilarity):
     """
     <url>:WMA link:https://reference.wolfram.com/language/ref/MatchingDissimilarity.html</url>
@@ -276,6 +376,31 @@ def eval(self, values, evaluation):
         return Expression(SymbolTimes, *elements)
 
 
+class PolygonalNumber(Builtin):
+    """
+    <url>
+    :Polygonal number: https://en.wikipedia.org/wiki/Polygonal_number</url> (<url>
+    :WMA: https://reference.wolfram.com/language/ref/PolygonalNumber.html</url>)
+    <dl>
+      <dt>'PolygonalNumber[$r$, $n$]'
+      <dd>gives the $n$th $r$-gonal number.
+    </dl>
+
+    >> Table[PolygonalNumber[n], {n, 10}]
+     = {1, 3, 6, 10, 15, 21, 28, 36, 45, 55}
+    >> Table[PolygonalNumber[r, 10], {r, 3, 10}]
+     = {55, 100, 145, 190, 235, 280, 325, 370}
+    """
+
+    attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
+    summary_text = "polygonal number"
+
+    rules = {
+        "PolygonalNumber[n_Integer]": "PolygonalNumber[3, n]",
+        "PolygonalNumber[r_Integer, n_Integer]": "(1/2) n (n (r - 2) - r + 4)",
+    }
+
+
 class RogersTanimotoDissimilarity(_BooleanDissimilarity):
     """
     <url>