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