diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 423f1e962..dfe88344c 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -714,7 +714,8 @@ the name, so `sum(eulergamma, eulergamma, 1, 3)` is `6` **Not silent: `catalannumber(n)` is the Catalan number, where it was a product.** It is `binomial(2n, n) / (n + 1)` for a whole `n >= 0`, and `4^n Γ(n + 1/2) / (√π Γ(n + 2))` elsewhere, as SymPy's `catalan` is: `-1/2` at `-1`, `0` at every whole number below it, and no value at a negative half-integer -([#2182](https://github.com/asc-community/AngouriMath/issues/2182)). +([#2182](https://github.com/asc-community/AngouriMath/issues/2182)). It compiles, on the same terms, in both compilers +([#2195](https://github.com/asc-community/AngouriMath/issues/2195)). | Input | Was (2.6.0) | Now | |---|---|---| diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs index 2b56e8204..18e5c3970 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs @@ -557,6 +557,17 @@ private protected override void CompileNode(Compiler compiler) // The binomial coefficient, the factorial powers and the double factorial, on the // interpreter's terms. https://github.com/asc-community/AngouriMath/issues/2195 + public partial record CatalanNumberf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_CATALAN_NUMBER)); + } + } + public partial record DoubleFactorialf { private protected override bool CompilesItself => true; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs index ba6963ad6..3f54d16c6 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs @@ -56,6 +56,7 @@ internal enum InstructionType CALL_IM, CALL_NOT, CALL_DOUBLE_FACTORIAL, + CALL_CATALAN_NUMBER, // 2-arg functions CALL_SUM = 100, @@ -367,6 +368,9 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value case InstructionType.CALL_HARMONIC_NUMBER: stack.Push(AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(stack.Pop(), stack.Pop())); break; + case InstructionType.CALL_CATALAN_NUMBER: + stack.Push(AngouriMath.Numerics.CombinatorialFunctions.CatalanNumber(stack.Pop())); + break; case InstructionType.CALL_DOUBLE_FACTORIAL: stack.Push(AngouriMath.Numerics.CombinatorialFunctions.DoubleFactorial(stack.Pop())); break; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs index 93bfbeac6..e7d4c8f89 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs @@ -157,6 +157,7 @@ Shif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Shi", Chif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Chi", 1, newE.Type), newE), Factorialf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Factorial", 1, newE.Type), newE), DoubleFactorialf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("DoubleFactorial", 1, newE.Type), newE), + CatalanNumberf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("CatalanNumber", 1, newE.Type), newE), Phif => Expression.Call(GetDef("Phi", 1, e.Type), e), // The conjugate, the argument and the parts, widened to double as the special diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs index bfa9ad7eb..5fd183372 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs @@ -335,6 +335,8 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a) // terms. https://github.com/asc-community/AngouriMath/issues/2195 public static System.Numerics.Complex DoubleFactorial(System.Numerics.Complex a) => AngouriMath.Numerics.CombinatorialFunctions.DoubleFactorial(a); public static double DoubleFactorial(double a) => AngouriMath.Numerics.CombinatorialFunctions.DoubleFactorial(a); + public static System.Numerics.Complex CatalanNumber(System.Numerics.Complex a) => AngouriMath.Numerics.CombinatorialFunctions.CatalanNumber(a); + public static double CatalanNumber(double a) => AngouriMath.Numerics.CombinatorialFunctions.CatalanNumber(a); public static System.Numerics.Complex Binomial(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.CombinatorialFunctions.Binomial(a, b); public static double Binomial(double a, double b) => AngouriMath.Numerics.CombinatorialFunctions.Binomial(a, b); public static System.Numerics.Complex FallingFactorial(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.CombinatorialFunctions.FallingFactorial(a, b); @@ -540,6 +542,8 @@ public static bool IsNaN(System.Numerics.Complex a) { ("Im", 1, typeof(double)), Def(a => Im(a)) }, { ("DoubleFactorial", 1, typeof(System.Numerics.Complex)), Def(a => DoubleFactorial(a)) }, { ("DoubleFactorial", 1, typeof(double)), Def(a => DoubleFactorial(a)) }, + { ("CatalanNumber", 1, typeof(System.Numerics.Complex)), Def(a => CatalanNumber(a)) }, + { ("CatalanNumber", 1, typeof(double)), Def(a => CatalanNumber(a)) }, { ("Binomial", 2, typeof(System.Numerics.Complex)), Def((a, b) => Binomial(a, b)) }, { ("Binomial", 2, typeof(double)), Def((a, b) => Binomial(a, b)) }, { ("FallingFactorial", 2, typeof(System.Numerics.Complex)), Def((a, b) => FallingFactorial(a, b)) }, diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt index d8f28e101..f6e87ece9 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt @@ -17,7 +17,7 @@ var parts = new[]{ "Conj", "Arg", "Re", "Im" }; // The combinatorial functions, for Complex and double alone in the same way: one argument, and two. // https://github.com/asc-community/AngouriMath/issues/2195 - var combinatorial = new[]{ "DoubleFactorial" }; + var combinatorial = new[]{ "DoubleFactorial", "CatalanNumber" }; var combinatorialBinary = new[]{ "Binomial", "FallingFactorial", "RisingFactorial" }; var types = new[]{ "System.Numerics.Complex", "double", "float", "long", "int", "BigInteger" }; #> diff --git a/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs b/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs index 5b4448da1..849ac540a 100644 --- a/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs +++ b/Sources/AngouriMath/Numerics/CombinatorialFunctions.cs @@ -119,6 +119,39 @@ public static NumericsComplex DoubleFactorial(NumericsComplex n) public static double DoubleFactorial(double n) => Real(DoubleFactorial(new NumericsComplex(n, 0))); + /// + /// The Catalan number: (2n)!/(n! (n + 1)!) for a whole n >= 0, -1/2 at -1, 0 at + /// a whole number below it, no value at a negative half-integer, and 4^n Γ(n + 1/2)/(√π Γ(n + 2)) + /// elsewhere. + /// + public static NumericsComplex CatalanNumber(NumericsComplex n) + { + if (IsWhole(n, out var whole)) + { + if (whole == -1) + return -0.5; + if (whole < -1) + return 0; + // Exactly, by C_(k + 1) = C_k 2 (2k + 1)/(k + 2), whose division leaves no remainder, and rounded + // once, since a product of the ratios in doubles rounds at each step and C_5 would be + // 41.99999999999999. Past the index the interpreter computes exactly, the gamma ratio. + if (whole <= 1000) + { + var catalan = System.Numerics.BigInteger.One; + for (var k = 0; k < whole; k++) + catalan = catalan * (2 * (2 * k + 1)) / (k + 2); + return (double)catalan; + } + } + if (IsWhole(n + 0.5, out var half) && half <= 0) + return NaN; + var top = SpecialFunctions.Factorial(n - 0.5); + var bottom = SpecialFunctions.Factorial(n + 1); + return NumericsComplex.Pow(4, n) * top / (Math.Sqrt(Math.PI) * bottom); + } + + public static double CatalanNumber(double n) => Real(CatalanNumber(new NumericsComplex(n, 0))); + // x (x - 1) ... (x - n + 1) falling and x (x + 1) ... (x + n - 1) rising for n > 0, and for // n < 0 the reciprocal of (x + 1) ... (x - n) falling and of (x - 1) ... (x + n) rising. private static NumericsComplex OfAWholeOrder(NumericsComplex x, long n, bool falling) diff --git a/Sources/Tests/AotSmokeTest/Program.cs b/Sources/Tests/AotSmokeTest/Program.cs index f1dc291f8..65000789b 100644 --- a/Sources/Tests/AotSmokeTest/Program.cs +++ b/Sources/Tests/AotSmokeTest/Program.cs @@ -313,6 +313,8 @@ private static void Compilation() CheckClose("compiled lambertw(-1, x) at x = -0.3", compiledW(-0.3), -1.7813370234216276); var compiledPolygamma = "polygamma(1, x) + harmonicnumber(x, 2)".Compile(x); CheckClose("compiled polygamma(1, x) + harmonicnumber(x, 2) at x = 1/2", compiledPolygamma(0.5), 4.934802200544679 + 0.7101318663035472); + var compiledCatalan = "catalannumber(x)".Compile(x); + Check("compiled catalannumber at 5 and -1", compiledCatalan.Call(5) == 42 && compiledCatalan.Call(-1) == -0.5, $"{compiledCatalan.Call(5)}, {compiledCatalan.Call(-1)}"); // FastExpression: a statement is its truth value, and a piecewise the case that holds. var fastPiecewise = "piecewise(x provided x > 0, -x)".Compile(x); var fastStatement = "x > 1 or x < 0".Compile(x); diff --git a/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs b/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs index c32662e32..146c02ca7 100644 --- a/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs +++ b/Sources/Tests/UnitTests/Common/CombinatorialCompilationTest.cs @@ -59,6 +59,34 @@ public void ATwoArgumentFunctionIsWhatEvaluationSays(string expression) Assert.True(compared > values.Length * values.Length * 3 / 4, $"only {compared} of {values.Length * values.Length} evaluated"); } + /// + /// The Catalan numbers: exact on the whole numbers, -1/2 at -1, 0 below it, no value at a + /// negative half-integer, and the gamma ratio elsewhere. + /// + [Fact] + public void TheCatalanNumbersAreWhatEvaluationSays() + { + var fast = "catalannumber(x)".ToEntity().Compile("x"); + var linq = "catalannumber(x)".ToEntity().Compile("x"); + var real = "catalannumber(x)".ToEntity().Compile("x"); + foreach (var a in values.Concat(new[] { "-4", "10", "30", "7/3", "-7/3" })) + { + var expected = Assert.IsAssignableFrom($"catalannumber({a})".ToEntity().Evaled).ToNumerics(); + var x = Numeric(a); + Assert.True(Agree(expected, fast.Call(x)), $"catalannumber({a}): {expected} evaluated, {fast.Call(x)} by FastExpression"); + Assert.True(Agree(expected, linq(x)), $"catalannumber({a}): {expected} evaluated, {linq(x)} by LINQ"); + if (x.Imaginary == 0 && expected.Imaginary == 0) + Assert.True(Agree(expected, real(x.Real)), $"catalannumber({a}): {expected} evaluated, {real(x.Real)} over double"); + } + } + + /// A Catalan number at a whole index is the integer, exactly, as far as a double holds it. + [Theory] + [InlineData(5, 42)] + [InlineData(30, 3814986502092304)] + public void ACatalanNumberAtAWholeIndexIsExact(int n, long expected) + => Assert.Equal((double)expected, "catalannumber(x)".ToEntity().Compile("x").Call(n).Real); + [Fact] public void TheDoubleFactorialIsWhatEvaluationSays() {