The Catalan numbers compile

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 |
 |---|---|---|

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;

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;

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

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<double>(a => Im(a)) },
             { ("DoubleFactorial", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => DoubleFactorial(a)) },
             { ("DoubleFactorial", 1, typeof(double)), Def<double>(a => DoubleFactorial(a)) },
+            { ("CatalanNumber", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => CatalanNumber(a)) },
+            { ("CatalanNumber", 1, typeof(double)), Def<double>(a => CatalanNumber(a)) },
             { ("Binomial", 2, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>((a, b) => Binomial(a, b)) },
             { ("Binomial", 2, typeof(double)), Def<double>((a, b) => Binomial(a, b)) },
             { ("FallingFactorial", 2, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>((a, b) => FallingFactorial(a, 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" };
 #>

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)));
 
+        /// <summary>
+        /// The Catalan number: <c>(2n)!/(n! (n + 1)!)</c> for a whole <c>n &gt;= 0</c>, <c>-1/2</c> at <c>-1</c>, <c>0</c> at
+        /// a whole number below it, no value at a negative half-integer, and <c>4^n Γ(n + 1/2)/(√π Γ(n + 2))</c>
+        /// elsewhere.
+        /// </summary>
+        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)

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<double, double>(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);

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");
         }
 
+        /// <summary>
+        /// The Catalan numbers: exact on the whole numbers, <c>-1/2</c> at <c>-1</c>, <c>0</c> below it, no value at a
+        /// negative half-integer, and the gamma ratio elsewhere.
+        /// </summary>
+        [Fact]
+        public void TheCatalanNumbersAreWhatEvaluationSays()
+        {
+            var fast = "catalannumber(x)".ToEntity().Compile("x");
+            var linq = "catalannumber(x)".ToEntity().Compile<Complex, Complex>("x");
+            var real = "catalannumber(x)".ToEntity().Compile<double, double>("x");
+            foreach (var a in values.Concat(new[] { "-4", "10", "30", "7/3", "-7/3" }))
+            {
+                var expected = Assert.IsAssignableFrom<Entity.Number.Complex>($"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");
+            }
+        }
+
+        /// <summary>A Catalan number at a whole index is the integer, exactly, as far as a double holds it.</summary>
+        [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()
         {