Skip to content

Commit 58a75c7

Browse files
Rafael-SOWNetclaude
Rafael-SOWNet
and
claude
authored
The factorial, gamma and the totient compile under Compile<TIn, TOut> (#1620)
Compile<TIn, TOut> had no form for x!, gamma or phi, and Compile's own had defects. Both compilers now call one routine in SpecialFunctions. The factorial is exact at a whole number up to 170, rounded once to a double, and infinite past it; NaN at a negative whole number, one of the gamma function's poles, as in the interpreter; and elsewhere the Lanczos approximation Compile always used (g = 7, nine terms), taken in logarithms, the reflection formula included. Compile's overflowed its own power before the value did, so from just past 141!, about 2e243, it compiled to NaN, up to 170!. At a whole number it was a rounding off the integer, 23.999999999999996 at 4, which a compilation over integers truncates to 23. Off the real line the reflection's sin(pi z) overflowed where the value is tiny; past |Im(pi z)| = 20 it is a logarithm written out now. gamma(x) is (x - 1)!. The totient is the totient at a whole number, 0 at one that is not positive, and NaN at any other number, as the interpreter's is. Compile truncated a fraction to the whole number below it, so phi(2.5) was phi(2). Under Compile<TIn, TOut>, long, int and BigInteger take it exactly. Against mpmath the compiled factorial is within 1.1e-13 of the value on a grid of [-14, 14]^2, 2.4e-13 up the imaginary axis to 10^4.5 and 3.2e-13 on circles out to radius 300. It is further out only where the conditioning forces it: next to a pole, and in the subnormal range past -177. SpecialFunctionCompilationTest gains the factorial and gamma against the interpreter, through both compilers and over doubles; whole numbers exact, including under Compile<long, long>; 171 infinite and the poles NaN; and the totient at whole numbers over every integer type, NaN elsewhere. UserInvalidExceptions keeps binomial, gcd, lcm and prime as the nodes with no compiled form. Closes #1607. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
1 parent 8cb650e commit 58a75c7

8 files changed

Lines changed: 271 additions & 38 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2357,6 +2357,27 @@ two of them, the points are samples of a continuum and not the set
23572357
| `"floor(x)".Compile("x")`, and `ceil`, `round`, `max`, `min` | `UncompilableNodeException` | compiled |
23582358
| `"max(x, 1)".ToEntity().Differentiate("x")` | `derivative(max(x, 1), x)` | `(1 + sgn(x - 1)) / 2 provided not x - 1 = 0` |
23592359

2360+
### The compiled factorial is exact at whole numbers and finite to `170!`, and the compiled totient has no value at a fraction
2361+
2362+
Both compilers now evaluate the factorial, the gamma function and Euler's totient the same way,
2363+
and `Compile<TIn, TOut>` has a form for them where it had none
2364+
([#1607](https://github.com/asc-community/AngouriMath/issues/1607)). The compiled factorial was
2365+
Lanczos's approximation everywhere, and its power overflowed before the value did: from just past
2366+
`141!`, about `2e243`, it compiled to NaN, up to `170!`. At a whole number it was off from the
2367+
integer by a rounding, which a compilation over integers would truncate. It is now the exact
2368+
factorial rounded once, up to `170!`; infinite past it; NaN at the negative whole numbers, the
2369+
poles; and the same approximation elsewhere, taken in logarithms so that it overflows only where
2370+
the value does. The compiled totient truncated a fraction to the whole number below it. There it
2371+
now has no value, as the interpreter's has none.
2372+
2373+
| Input | Was (2.5.0) | Now |
2374+
|---|---|---|
2375+
| `"x!".Compile("x").Call(4)` | `23.999999999999996` | `24` |
2376+
| `"x!".Compile("x").Call(0)` | `0.9999999999999998` | `1` |
2377+
| `"x!".Compile("x").Call(142)`, and on to `170` | NaN | `2.695364137888163E+245` |
2378+
| `"phi(x)".Compile("x").Call(2.5)` | `1` | NaN |
2379+
| `"x!".ToEntity().Compile<double, double>("x")`, and `gamma`, `phi` | `UncompilableNodeException` | compiled |
2380+
23602381
### An equation the solver cannot invert is left unsolved, not answered with no roots
23612382

23622383
`x! = 6` was answered `{ }`, a claim that it has no roots, and it has 3. The solver isolates `x`

‎Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs‎

Lines changed: 2 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -258,24 +258,7 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value
258258
stack.Push(AngouriMath.Numerics.SpecialFunctions.Chi(stack.Pop()));
259259
break;
260260
case InstructionType.CALL_FACTORIAL:
261-
// https://stackoverflow.com/a/15454784/5429648
262-
const int g = 7;
263-
static System.Numerics.Complex Gamma(System.Numerics.Complex z)
264-
{
265-
if (z.Real < 0.5) return System.Math.PI / (System.Numerics.Complex.Sin(System.Math.PI * z) * Gamma(1 - z));
266-
else
267-
{
268-
z -= 1;
269-
270-
System.Numerics.Complex x = gammaCoeffs[0];
271-
for (var i = 1; i < g + 2; i++)
272-
x += gammaCoeffs[i] / (z + i);
273-
274-
var t = z + g + 0.5;
275-
return System.Math.Sqrt(2 * System.Math.PI) * System.Numerics.Complex.Pow(t, z + 0.5) * System.Numerics.Complex.Exp(-t) * x;
276-
}
277-
}
278-
stack.Push(Gamma(stack.Pop() + 1));
261+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Factorial(stack.Pop()));
279262
break;
280263
case InstructionType.CALL_SIGNUM:
281264
stack.Push(stack.Pop().Signum());
@@ -284,21 +267,14 @@ static System.Numerics.Complex Gamma(System.Numerics.Complex z)
284267
stack.Push(stack.Pop().Abs());
285268
break;
286269
case InstructionType.CALL_PHI:
287-
var n = (long) stack.Pop().Real;
288-
stack.Push(n.Phi());
270+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Phi(stack.Pop()));
289271
break;
290272
}
291273
if (stack.Count != 1)
292274
throw new AngouriBugException($"Unused values remain in the stack: {stack.Count} instead of 1");
293275
return stack.Pop();
294276
}
295277

296-
[ConstantField] static readonly double[] gammaCoeffs = {
297-
0.99999999999980993, 676.5203681218851, -1259.1392167224028,
298-
771.32342877765313, -176.61502916214059, 12.507343278686905,
299-
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
300-
};
301-
302278
/// <summary>Might be useful for debug if a function works too slowly</summary>
303279
public override string ToString() => string.Join(" \n| ", instructions);
304280
}

‎Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -149,6 +149,8 @@ Sif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Si",
149149
Cif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Ci", 1, newE.Type), newE),
150150
Shif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Shi", 1, newE.Type), newE),
151151
Chif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Chi", 1, newE.Type), newE),
152+
Factorialf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Factorial", 1, newE.Type), newE),
153+
Phif => Expression.Call(GetDef("Phi", 1, e.Type), e),
152154

153155
Notf => Expression.Not(e),
154156

‎Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs‎

Lines changed: 21 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -288,8 +288,17 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a)
288288
public static int Min(int a, int b) => Math.Min(a, b);
289289
public static BigInteger Min(BigInteger a, BigInteger b) => BigInteger.Min(a, b);
290290

291-
// The special functions, in double precision as the kernels define them.
292-
// https://github.com/asc-community/AngouriMath/issues/1607
291+
// Euler's totient: exact over the integers, and over Complex and double a value at a whole
292+
// number and NaN at any other. https://github.com/asc-community/AngouriMath/issues/1607
293+
public static System.Numerics.Complex Phi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
294+
public static double Phi(double a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
295+
public static float Phi(float a) => (float)AngouriMath.Numerics.SpecialFunctions.Phi(a);
296+
public static long Phi(long a) => a.Phi();
297+
public static int Phi(int a) => (int)((long)a).Phi();
298+
public static BigInteger Phi(BigInteger a) => PeterO.Numbers.EInteger.FromBytes(a.ToByteArray(), littleEndian: true).Phi().ToBigInteger();
299+
300+
// The special functions, in double precision as the kernels define them, and the factorial,
301+
// the gamma function one along. https://github.com/asc-community/AngouriMath/issues/1607
293302
public static System.Numerics.Complex Erf(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erf(a);
294303
public static double Erf(double a) => AngouriMath.Numerics.SpecialFunctions.Erf(a);
295304
public static System.Numerics.Complex Erfc(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erfc(a);
@@ -308,6 +317,8 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a)
308317
public static double Shi(double a) => AngouriMath.Numerics.SpecialFunctions.Shi(a);
309318
public static System.Numerics.Complex Chi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Chi(a);
310319
public static double Chi(double a) => AngouriMath.Numerics.SpecialFunctions.Chi(a);
320+
public static System.Numerics.Complex Factorial(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Factorial(a);
321+
public static double Factorial(double a) => AngouriMath.Numerics.SpecialFunctions.Factorial(a);
311322

312323
/// <summary>
313324
/// True where <c>System.Numerics.Complex.IsNaN</c> is, spelled out because that overload is
@@ -459,6 +470,12 @@ public static bool IsNaN(System.Numerics.Complex a)
459470
{ ("Round", 1, typeof(long)), Def<long>(a => Round(a)) },
460471
{ ("Round", 1, typeof(int)), Def<int>(a => Round(a)) },
461472
{ ("Round", 1, typeof(BigInteger)), Def<BigInteger>(a => Round(a)) },
473+
{ ("Phi", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Phi(a)) },
474+
{ ("Phi", 1, typeof(double)), Def<double>(a => Phi(a)) },
475+
{ ("Phi", 1, typeof(float)), Def<float>(a => Phi(a)) },
476+
{ ("Phi", 1, typeof(long)), Def<long>(a => Phi(a)) },
477+
{ ("Phi", 1, typeof(int)), Def<int>(a => Phi(a)) },
478+
{ ("Phi", 1, typeof(BigInteger)), Def<BigInteger>(a => Phi(a)) },
462479
{ ("Erf", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Erf(a)) },
463480
{ ("Erf", 1, typeof(double)), Def<double>(a => Erf(a)) },
464481
{ ("Erfc", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Erfc(a)) },
@@ -477,6 +494,8 @@ public static bool IsNaN(System.Numerics.Complex a)
477494
{ ("Shi", 1, typeof(double)), Def<double>(a => Shi(a)) },
478495
{ ("Chi", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Chi(a)) },
479496
{ ("Chi", 1, typeof(double)), Def<double>(a => Chi(a)) },
497+
{ ("Factorial", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Factorial(a)) },
498+
{ ("Factorial", 1, typeof(double)), Def<double>(a => Factorial(a)) },
480499
};
481500

482501
// An expression tree is how C# spells "the MethodInfo of this method" without a string:

‎Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt‎

Lines changed: 13 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -7,11 +7,11 @@
77
// to one and forgotten in the other.
88
var unary = new[]{ "Sin", "Cos", "Tan", "Asin", "Acos", "Atan",
99
"Cot", "Sec", "Csc", "Acot", "Asec", "Acsc",
10-
"Abs", "Sgn", "Floor", "Ceil", "Round" };
10+
"Abs", "Sgn", "Floor", "Ceil", "Round", "Phi" };
1111
var binary = new[]{ "Log", "Pow", "Max", "Min" };
1212
// Defined for Complex and double alone: an integer or a float is widened to double first, as
1313
// for the trigonometry.
14-
var special = new[]{ "Erf", "Erfc", "Erfi", "Ei", "Li", "Si", "Ci", "Shi", "Chi" };
14+
var special = new[]{ "Erf", "Erfc", "Erfi", "Ei", "Li", "Si", "Ci", "Shi", "Chi", "Factorial" };
1515
var types = new[]{ "System.Numerics.Complex", "double", "float", "long", "int", "BigInteger" };
1616
#>
1717
//
@@ -191,8 +191,17 @@ namespace AngouriMath.Core.Compilation.IntoLinq
191191
public static int Min(int a, int b) => Math.Min(a, b);
192192
public static BigInteger Min(BigInteger a, BigInteger b) => BigInteger.Min(a, b);
193193

194-
// The special functions, in double precision as the kernels define them.
195-
// https://github.com/asc-community/AngouriMath/issues/1607
194+
// Euler's totient: exact over the integers, and over Complex and double a value at a whole
195+
// number and NaN at any other. https://github.com/asc-community/AngouriMath/issues/1607
196+
public static System.Numerics.Complex Phi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
197+
public static double Phi(double a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
198+
public static float Phi(float a) => (float)AngouriMath.Numerics.SpecialFunctions.Phi(a);
199+
public static long Phi(long a) => a.Phi();
200+
public static int Phi(int a) => (int)((long)a).Phi();
201+
public static BigInteger Phi(BigInteger a) => PeterO.Numbers.EInteger.FromBytes(a.ToByteArray(), littleEndian: true).Phi().ToBigInteger();
202+
203+
// The special functions, in double precision as the kernels define them, and the factorial,
204+
// the gamma function one along. https://github.com/asc-community/AngouriMath/issues/1607
196205
<# foreach (var name in special) { #>
197206
public static System.Numerics.Complex <#= name #>(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.<#= name #>(a);
198207
public static double <#= name #>(double a) => AngouriMath.Numerics.SpecialFunctions.<#= name #>(a);

‎Sources/AngouriMath/Numerics/SpecialFunctions.cs‎

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -380,6 +380,120 @@ private static NumericsComplex HalfOfTheExponentialSeries(NumericsComplex z, boo
380380
return sum;
381381
}
382382

383+
/// <summary>
384+
/// <c>x!</c>, the gamma function one along. At a whole number up to 170 it is the exact
385+
/// factorial rounded once, and past 170, where a double cannot hold it, it is infinite. At a
386+
/// negative whole number, one of the gamma function's poles, it is NaN, as the interpreter's
387+
/// is. Elsewhere it is Lanczos's approximation with <c>g = 7</c> and nine terms, the one the
388+
/// compiler has always used. Its <c>t^(z + 1/2) e^(-t)</c> is taken in logarithms, so that
389+
/// it overflows only where the value does: the power alone passed the largest double just past
390+
/// 141, where the factorial is about <c>2e243</c>, and the compiled factorial was NaN from
391+
/// there to 170.
392+
/// https://github.com/asc-community/AngouriMath/issues/1607
393+
/// </summary>
394+
public static NumericsComplex Factorial(NumericsComplex z)
395+
{
396+
if (z.Imaginary == 0)
397+
return Factorial(z.Real);
398+
if (IsNaN(z) || IsInfinite(z))
399+
return NaN;
400+
return Gamma(new NumericsComplex(z.Real + 1, z.Imaginary));
401+
}
402+
403+
public static double Factorial(double x)
404+
{
405+
if (double.IsNaN(x) || double.IsNegativeInfinity(x))
406+
return double.NaN;
407+
if (double.IsPositiveInfinity(x))
408+
return double.PositiveInfinity;
409+
if (x == Math.Floor(x))
410+
return x < 0 ? double.NaN : x <= 170 ? factorials[(int)x] : double.PositiveInfinity;
411+
return Gamma(x + 1);
412+
}
413+
414+
/// <summary><c>0!</c> to <c>170!</c>, each the exact integer rounded once to a double.</summary>
415+
[AngouriMath.Core.ConstantField] private static readonly double[] factorials = ExactFactorials();
416+
417+
private static double[] ExactFactorials()
418+
{
419+
var table = new double[171];
420+
var factorial = PeterO.Numbers.EInteger.One;
421+
for (var n = 0; n <= 170; n++)
422+
{
423+
if (n > 0)
424+
factorial *= n;
425+
table[n] = PeterO.Numbers.EFloat.FromEInteger(factorial).ToDouble();
426+
}
427+
return table;
428+
}
429+
430+
[AngouriMath.Core.ConstantField] private static readonly double[] lanczos = {
431+
0.99999999999980993, 676.5203681218851, -1259.1392167224028,
432+
771.32342877765313, -176.61502916214059, 12.507343278686905,
433+
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 };
434+
435+
private static double Gamma(double x)
436+
{
437+
// The reflection formula below 1/2.
438+
if (x < 0.5)
439+
return Math.PI / (Math.Sin(Math.PI * x) * Gamma(1 - x));
440+
x -= 1;
441+
var sum = lanczos[0];
442+
for (var i = 1; i < 9; i++)
443+
sum += lanczos[i] / (x + i);
444+
var t = x + 7.5;
445+
return Math.Exp((x + 0.5) * Math.Log(t) - t + Math.Log(Math.Sqrt(2 * Math.PI) * sum));
446+
}
447+
448+
/// <summary>
449+
/// The gamma function off the real line, as the exponential of its logarithm, so that it
450+
/// overflows and underflows only where the value does. The reflection formula
451+
/// <c>pi/(sin(pi z) gamma(1 - z))</c> is taken in logarithms too: far from the real line
452+
/// <c>sin(pi z)</c> overflows where the value it divides is tiny, and it made a NaN.
453+
/// </summary>
454+
private static NumericsComplex Gamma(NumericsComplex z)
455+
=> z.Real < 0.5
456+
? ExpTimes(Math.Log(Math.PI) - LogSine(Math.PI * z) - LogGamma(1 - z), 1)
457+
: ExpTimes(LogGamma(z), 1);
458+
459+
/// <summary>A logarithm of the gamma function for <c>Re z &gt;= 1/2</c>, by Lanczos's approximation.</summary>
460+
private static NumericsComplex LogGamma(NumericsComplex z)
461+
{
462+
z -= 1;
463+
NumericsComplex sum = lanczos[0];
464+
for (var i = 1; i < 9; i++)
465+
sum += lanczos[i] / (z + i);
466+
var t = z + 7.5;
467+
return (z + 0.5) * NumericsComplex.Log(t) - t + NumericsComplex.Log(Math.Sqrt(2 * Math.PI) * sum);
468+
}
469+
470+
/// <summary>
471+
/// A logarithm of <c>sin w</c>. Past <c>|Im w| = 20</c> one of <c>e^(i w)</c> and
472+
/// <c>e^(-i w)</c> in <c>sin w = (e^(i w) - e^(-i w))/(2 i)</c> is below <c>1e-17</c> of the
473+
/// other. Its logarithm is then written out, where <c>sin w</c> itself would overflow.
474+
/// </summary>
475+
private static NumericsComplex LogSine(NumericsComplex w)
476+
{
477+
if (Math.Abs(w.Imaginary) < 20)
478+
return NumericsComplex.Log(NumericsComplex.Sin(w));
479+
// (i/2) e^(-i w) above the real line, and -(i/2) e^(i w) below it.
480+
return w.Imaginary > 0
481+
? new NumericsComplex(w.Imaginary - Math.Log(2), Math.PI / 2 - w.Real)
482+
: new NumericsComplex(-w.Imaginary - Math.Log(2), w.Real - Math.PI / 2);
483+
}
484+
485+
/// <summary>
486+
/// Euler's totient at a whole number, 0 at one that is not positive, as the interpreter has
487+
/// it, and NaN at any other number, where the interpreter's is NaN too. The compiler
488+
/// truncated a fraction to the whole number below it, so <c>phi(2.5)</c> was <c>phi(2)</c>.
489+
/// Past <c>long</c>'s range it is NaN: the interpreter's totient cannot go there either.
490+
/// https://github.com/asc-community/AngouriMath/issues/1607
491+
/// </summary>
492+
public static NumericsComplex Phi(NumericsComplex z) => z.Imaginary == 0 ? Phi(z.Real) : NaN;
493+
494+
public static double Phi(double x)
495+
=> x == Math.Floor(x) && x >= long.MinValue && x < long.MaxValue ? ((long)x).Phi() : double.NaN;
496+
383497
/// <summary>
384498
/// Half of <paramref name="value"/>, part by part: dividing a complex number by 2 multiplies it
385499
/// by the complex number 2, and an infinite part times the other's zero makes a NaN.

0 commit comments

Comments
 (0)