From 07282cf8f225d1c3470b7090242e7855a2093d02 Mon Sep 17 00:00:00 2001 From: David A Roberts Date: Sun, 20 Oct 2024 11:16:17 +1000 Subject: [PATCH 1/2] Add combinatorial functions --- mathics/builtin/intfns/combinatorial.py | 209 ++++++++++++++++++++++++ 1 file changed, 209 insertions(+) diff --git a/mathics/builtin/intfns/combinatorial.py b/mathics/builtin/intfns/combinatorial.py index cc1a30989..b53560065 100644 --- a/mathics/builtin/intfns/combinatorial.py +++ b/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): + """ + + :Bell number: https://en.wikipedia.org/wiki/Bell_number ( + :SymPy: https://docs.sympy.org/latest/modules/functions/combinatorial.html#sympy.functions.combinatorial.numbers.bell, + :WMA: https://reference.wolfram.com/language/ref/BellB.html) +
+
'BellB[$n$]' +
Bell number $B$_$n$. + +
'BellB[$n$, $x$]' +
Bell polynomial $B$_$n$($x$). +
+ + >> 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): + """ + + :Euler numbers: https://en.wikipedia.org/wiki/Euler_numbers ( + :SymPy: https://docs.sympy.org/latest/modules/functions/combinatorial.html#sympy.functions.combinatorial.numbers.euler, + :WMA: https://reference.wolfram.com/language/ref/EulerE.html) +
+
'EulerE[$n$]' +
Euler number $E$_$n$. + +
'EulerE[$n$, $x$]' +
Euler polynomial $E$_$n$($x$). +
+ + >> 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): """ @@ -214,6 +276,93 @@ def _compute(self, n, c_ff, c_ft, c_tf, c_tt): ) +class JacobiSymbol(Builtin): + """ + + :Jacobi symbol: https://en.wikipedia.org/wiki/Jacobi_symbol ( + :WMA: https://reference.wolfram.com/language/ref/JacobiSymbol.html) +
+
'JacobiSymbol[$a$, $n$]' +
returns the Jacobi symbol ($a$/$n$). +
+ + >> Table[JacobiSymbol[n, m], {n, 0, 10}, {m, 1, n, 2}] + = {{}, {1}, {1}, {1, 0}, {1, 1}, {1, -1, 0}, {1, 0, 1}, {1, 1, -1, 0}, {1, -1, -1, 1}, {1, 0, 1, 1, 0}, {1, 1, 0, -1, 1}} + """ + + attributes = A_LISTABLE | A_PROTECTED + summary_text = "Jacobi symbol" + + rules = { + "JacobiSymbol[a_, p_?PrimeQ]": "Which[Mod[a, p] == 0, 0, PowerMod[a, (p - 1)/2, p] == 1, 1, True, -1]", # Legendre symbol + "JacobiSymbol[a_, n_]": "Times @@ (JacobiSymbol[a, #1]^#2 & @@@ FactorInteger[n])", + } + + +class KroneckerSymbol(Builtin): + """ + + :Kronecker symbol: https://en.wikipedia.org/wiki/Kronecker_symbol ( + :WMA: https://reference.wolfram.com/language/ref/KroneckerSymbol.html) +
+
'KroneckerSymbol[$a$, $n$]' +
returns the Kronecker symbol ($a$/$n$). +
+ + >> Table[KroneckerSymbol[n, m], {n, 5}, {m, 5}] + = {{1, 1, 1, 1, 1}, {1, 0, -1, 0, -1}, {1, -1, 0, 1, -1}, {1, 0, 1, 0, 1}, {1, -1, -1, 1, 0}} + """ + + attributes = A_LISTABLE | A_PROTECTED | A_READ_PROTECTED + summary_text = "Kronecker symbol" + + rules = { + "KroneckerSymbol[a_, n_?(Positive[#] && OddQ[#] &)]": "JacobiSymbol[a, n]", + "KroneckerSymbol[a_, 0]": "If[Abs[a] == 1, 1, 0]", + "KroneckerSymbol[a_, -1]": "If[a < 0, -1, 1]", + "KroneckerSymbol[a_, 2]": "Which[EvenQ[a], 0, Mod[a, 8] == 1 || Mod[a, 8] == 7, 1, True, -1]", + "KroneckerSymbol[a_, n_]": "Times @@ (KroneckerSymbol[a, #1]^#2 & @@@ FactorInteger[n])", + } + + +class LucasL(SympyFunction): + """ + + :Lucas Number: + https://en.wikipedia.org/wiki/Lucas_number ( + :SymPy: + https://docs.sympy.org/latest/modules/functions/combinatorial.html#sympy.functions.combinatorial.numbers.lucas, \ + + :WMA: + https://reference.wolfram.com/language/ref/LucasL.html) + +
+
'LucasL[$n$]' +
gives the $n$th Lucas number. +
+ + 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): """ :WMA link:https://reference.wolfram.com/language/ref/MatchingDissimilarity.html @@ -276,6 +425,31 @@ def eval(self, values, evaluation): return Expression(SymbolTimes, *elements) +class PolygonalNumber(Builtin): + """ + + :Polygonal number: https://en.wikipedia.org/wiki/Polygonal_number ( + :WMA: https://reference.wolfram.com/language/ref/PolygonalNumber.html) +
+
'PolygonalNumber[$r$, $n$]' +
gives the $n$th $r$-gonal number. +
+ + >> 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): """ @@ -352,6 +526,41 @@ def _compute(self, n, c_ff, c_ft, c_tf, c_tt): return Expression(SymbolDivide, Integer(r), Integer(c_tt + r)) +class SquaresR(Builtin): + """ + + :Sum of squares function: https://en.wikipedia.org/wiki/Sum_of_squares_function ( + :WMA: https://reference.wolfram.com/language/ref/SquaresR.html) + +
+
'SquaresR[$d$, $n$]' +
returns the number of ways to represent $n$ as a sum of $d$ squares. +
+ + >> Table[SquaresR[2, n], {n, 10}] + = {4, 4, 0, 4, 8, 0, 0, 4, 4, 8} + >> Table[Sum[SquaresR[2, k], {k, 0, n^2}], {n, 5}] + = {5, 13, 29, 49, 81} + >> Table[SquaresR[4, n], {n, 10}] + = {8, 24, 32, 24, 48, 96, 64, 24, 104, 144} + >> Table[SquaresR[6, n], {n, 10}] + = {12, 60, 160, 252, 312, 544, 960, 1020, 876, 1560} + >> Table[SquaresR[8, n], {n, 10}] + = {16, 112, 448, 1136, 2016, 3136, 5504, 9328, 12112, 14112} + """ + + attributes = A_LISTABLE | A_PROTECTED | A_READ_PROTECTED + summary_text = "sum of squares function" + + rules = { + "SquaresR[d_Integer, 0]": "1", + "SquaresR[2, n_Integer?Positive]": "4 Total[(-1)^((# - 1)/2) & /@ Select[Divisors[n], Mod[#, 4] == 1 || Mod[#, 4] == 3 &]]", + "SquaresR[4, n_Integer?Positive]": "8 Total[Select[Divisors[n], Mod[#, 4] != 0 &]]", + "SquaresR[6, n_Integer?Positive]": "4 Total[#^2 * (4 * KroneckerSymbol[-4, n/#] - KroneckerSymbol[-4, #]) & /@ Divisors[n]]", + "SquaresR[8, n_Integer?Positive]": "16 Total[(-1)^(n + #) #^3 & /@ Divisors[n]]", + } + + class Subsets(Builtin): """ From e0e88323a9ab4a1b80a7f55732f83f23903b8abb Mon Sep 17 00:00:00 2001 From: David A Roberts Date: Sun, 20 Oct 2024 12:03:13 +1000 Subject: [PATCH 2/2] Move number theory functions --- mathics/builtin/intfns/combinatorial.py | 84 ------------------------- 1 file changed, 84 deletions(-) diff --git a/mathics/builtin/intfns/combinatorial.py b/mathics/builtin/intfns/combinatorial.py index b53560065..33b7333f2 100644 --- a/mathics/builtin/intfns/combinatorial.py +++ b/mathics/builtin/intfns/combinatorial.py @@ -276,55 +276,6 @@ def _compute(self, n, c_ff, c_ft, c_tf, c_tt): ) -class JacobiSymbol(Builtin): - """ - - :Jacobi symbol: https://en.wikipedia.org/wiki/Jacobi_symbol ( - :WMA: https://reference.wolfram.com/language/ref/JacobiSymbol.html) -
-
'JacobiSymbol[$a$, $n$]' -
returns the Jacobi symbol ($a$/$n$). -
- - >> Table[JacobiSymbol[n, m], {n, 0, 10}, {m, 1, n, 2}] - = {{}, {1}, {1}, {1, 0}, {1, 1}, {1, -1, 0}, {1, 0, 1}, {1, 1, -1, 0}, {1, -1, -1, 1}, {1, 0, 1, 1, 0}, {1, 1, 0, -1, 1}} - """ - - attributes = A_LISTABLE | A_PROTECTED - summary_text = "Jacobi symbol" - - rules = { - "JacobiSymbol[a_, p_?PrimeQ]": "Which[Mod[a, p] == 0, 0, PowerMod[a, (p - 1)/2, p] == 1, 1, True, -1]", # Legendre symbol - "JacobiSymbol[a_, n_]": "Times @@ (JacobiSymbol[a, #1]^#2 & @@@ FactorInteger[n])", - } - - -class KroneckerSymbol(Builtin): - """ - - :Kronecker symbol: https://en.wikipedia.org/wiki/Kronecker_symbol ( - :WMA: https://reference.wolfram.com/language/ref/KroneckerSymbol.html) -
-
'KroneckerSymbol[$a$, $n$]' -
returns the Kronecker symbol ($a$/$n$). -
- - >> Table[KroneckerSymbol[n, m], {n, 5}, {m, 5}] - = {{1, 1, 1, 1, 1}, {1, 0, -1, 0, -1}, {1, -1, 0, 1, -1}, {1, 0, 1, 0, 1}, {1, -1, -1, 1, 0}} - """ - - attributes = A_LISTABLE | A_PROTECTED | A_READ_PROTECTED - summary_text = "Kronecker symbol" - - rules = { - "KroneckerSymbol[a_, n_?(Positive[#] && OddQ[#] &)]": "JacobiSymbol[a, n]", - "KroneckerSymbol[a_, 0]": "If[Abs[a] == 1, 1, 0]", - "KroneckerSymbol[a_, -1]": "If[a < 0, -1, 1]", - "KroneckerSymbol[a_, 2]": "Which[EvenQ[a], 0, Mod[a, 8] == 1 || Mod[a, 8] == 7, 1, True, -1]", - "KroneckerSymbol[a_, n_]": "Times @@ (KroneckerSymbol[a, #1]^#2 & @@@ FactorInteger[n])", - } - - class LucasL(SympyFunction): """ @@ -526,41 +477,6 @@ def _compute(self, n, c_ff, c_ft, c_tf, c_tt): return Expression(SymbolDivide, Integer(r), Integer(c_tt + r)) -class SquaresR(Builtin): - """ - - :Sum of squares function: https://en.wikipedia.org/wiki/Sum_of_squares_function ( - :WMA: https://reference.wolfram.com/language/ref/SquaresR.html) - -
-
'SquaresR[$d$, $n$]' -
returns the number of ways to represent $n$ as a sum of $d$ squares. -
- - >> Table[SquaresR[2, n], {n, 10}] - = {4, 4, 0, 4, 8, 0, 0, 4, 4, 8} - >> Table[Sum[SquaresR[2, k], {k, 0, n^2}], {n, 5}] - = {5, 13, 29, 49, 81} - >> Table[SquaresR[4, n], {n, 10}] - = {8, 24, 32, 24, 48, 96, 64, 24, 104, 144} - >> Table[SquaresR[6, n], {n, 10}] - = {12, 60, 160, 252, 312, 544, 960, 1020, 876, 1560} - >> Table[SquaresR[8, n], {n, 10}] - = {16, 112, 448, 1136, 2016, 3136, 5504, 9328, 12112, 14112} - """ - - attributes = A_LISTABLE | A_PROTECTED | A_READ_PROTECTED - summary_text = "sum of squares function" - - rules = { - "SquaresR[d_Integer, 0]": "1", - "SquaresR[2, n_Integer?Positive]": "4 Total[(-1)^((# - 1)/2) & /@ Select[Divisors[n], Mod[#, 4] == 1 || Mod[#, 4] == 3 &]]", - "SquaresR[4, n_Integer?Positive]": "8 Total[Select[Divisors[n], Mod[#, 4] != 0 &]]", - "SquaresR[6, n_Integer?Positive]": "4 Total[#^2 * (4 * KroneckerSymbol[-4, n/#] - KroneckerSymbol[-4, #]) & /@ Divisors[n]]", - "SquaresR[8, n_Integer?Positive]": "16 Total[(-1)^(n + #) #^3 & /@ Divisors[n]]", - } - - class Subsets(Builtin): """