From e8a95f4223e6e5cd60148e648843186736bb2f0c Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 23 Oct 2026 10:41:59 +0000 Subject: [PATCH] Make Spouge's constants from exact ratios and powers of e The factorial of a non-integer is Spouge's approximation, with a = 1.3 P constants at P digits, cached per precision. Each constant c_k = (-1)^(k - 1) (a - k)^(k - 1/2) e^(a - k)/(k - 1)! was made with a power to a non-integer exponent and an exponential, a logarithm and two exponentials at 1.95 P digits for each of the 1.3 P terms. The first factorial at a precision took 0.9 s at 100 digits, 20.5 s at 300 and 218 s at 600. (a - k)^(k - 1)/(k - 1)! is a ratio of exact integers, sqrt(a - k) one square root, and e^(a - k) the running product of e, carried with ten more digits for the a multiplications and rounded back to the constants' precision. The first factorial now takes 0.12 s at 100 digits, 0.57 s at 300, 1.5 s at 600, 5.5 s at 1000 and 30 s at 2000. The values are unchanged. Against mpmath, Gamma(1/3), Gamma(1/3 + i), Gamma(2/7), Gamma(-7/3), Gamma(101/2) and Gamma(-5/2 + 3i) have the same error on both builds at 100 and 300 digits, within 1.3 units of the last digit, and stay within that at 600 and 1000. Once compiled, an evaluation at 100 digits costs the same, 7 to 10 ms, on both. A test of (1/3)! and (1/3 + i)! at 300 digits is added. Part of #2075. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Functions/InternalAMExtensions.cs | 25 ++++++++++++------- .../Core/FactorialFarFromZeroTest.cs | 20 ++++++++++++++- 2 files changed, 35 insertions(+), 10 deletions(-) diff --git a/Sources/AngouriMath/Functions/InternalAMExtensions.cs b/Sources/AngouriMath/Functions/InternalAMExtensions.cs index bdf772f15..b41a2d8f6 100644 --- a/Sources/AngouriMath/Functions/InternalAMExtensions.cs +++ b/Sources/AngouriMath/Functions/InternalAMExtensions.cs @@ -1245,17 +1245,24 @@ internal static EDecimal[] GetSpougeFactorialConstants(int a) constants[0] = EDecimal.PI(mc).Multiply(2, mc).Sqrt(mc); - var negative = false; + // c_k = (-1)^(k - 1) (a - k)^(k - 1)/(k - 1)! sqrt(a - k) e^(a - k), the first factor a ratio of exact + // integers and e^(a - k) a product of e's, with digits to spare for the a multiplications that make it. + // With a power of a - k to k - 1/2 and an exponential for each k, the constants for 300 digits took 20 + // seconds, and for 600 digits 218. + var work = EContext.ForPrecision(a * 15 / 10 + 10); + var e = EDecimal.One.Exp(work); + var powersOfE = new EDecimal[a]; + powersOfE[0] = EDecimal.One; + for (int m = 1; m < a; m++) + powersOfE[m] = powersOfE[m - 1].Multiply(e, work); for (int k = 1; k < a; k++) { - var deltaAK = EDecimal.FromInt32(a - k); - var ck = deltaAK.Pow(EDecimal.FromInt32(k).Subtract(0.5m), mc); - ck = deltaAK.Exp(mc).Multiply(ck, mc); - ck = ck.Divide(Factorial(k - 1), mc); - if (negative) - ck = ck.Negate(); - constants[k] = ck; - negative = !negative; + var m = a - k; + var ck = EDecimal.FromEInteger(EInteger.FromInt32(m).Pow(k - 1)).Divide(EDecimal.FromEInteger(Factorial(k - 1)), work) + .Multiply(EDecimal.FromInt32(m).Sqrt(work), work) + .Multiply(powersOfE[m], work) + .RoundToPrecision(mc); + constants[k] = k % 2 == 0 ? ck.Negate() : ck; } return constants; diff --git a/Sources/Tests/UnitTests/Core/FactorialFarFromZeroTest.cs b/Sources/Tests/UnitTests/Core/FactorialFarFromZeroTest.cs index a2b094839..a06b06023 100644 --- a/Sources/Tests/UnitTests/Core/FactorialFarFromZeroTest.cs +++ b/Sources/Tests/UnitTests/Core/FactorialFarFromZeroTest.cs @@ -1,4 +1,4 @@ -// +// // Copyright (c) 2019-2026 Angouri. // AngouriMath is licensed under MIT. // Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. @@ -78,6 +78,24 @@ public void ADivisorsModulusNearAWholeNumberIsNotRoundedToIt() Close("-0.002993872340987814525483818825494753213055396187438519990678101458851071", value.ImaginaryPart.EDecimal, 60); } + /// + /// At 300 digits, on the real line and off it. Spouge's 390 constants were made with a power and an exponential + /// each, and the first factorial at that precision took 20 seconds; they are a ratio of integers, a square root and + /// a product of e's now. + /// + [Fact] + public void AtThreeHundredDigits() + { + var context = EContext.ForPrecision(300); + Close("0.89297951156924921121856431365822588137622979265243370031680944253013920338924793984699429634706292980638634973335742111190636152316815741999257112256948992272242994101666773944751677408238564727251848331757626146492648043525527613387077304719882217534457586134019354517131441367733861313846090028508017368", + Evaluate("(1/3)!", context).RealPart.EDecimal, 295); + var value = Evaluate("(1/3 + i)!", context); + Close("0.55310533708451368643420143855148996023577256508846542223376506819123717723035945992089147907370521801219397887942272059738687781730094054546525636372408182990085240891013953705988378166744479672870436455961985418938874599630727253813312731648659256268261654562175265755195218626689594794087234395264053064", + value.RealPart.EDecimal, 295); + Close("0.0097530148049768242524751667053160295250646412881078165005696874705860541880090280173509525760804027158706765197226375317213454843960898145149813764056673011748957471092198259369271410261828558753507851947578513569017247888018679994921424076865568363297838883827120115940210102202849352569014259005743340527", + value.ImaginaryPart.EDecimal, 293); + } + /// The reflection formula off the real line, which the complex factorial now takes left of zero. [Fact] public void LeftOfZeroOffTheRealLine()