Skip to content

The factorial, gamma and the totient compile under Compile<TIn, TOut> #1620

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Sep 30, 2026
Merged

The factorial, gamma and the totient compile under Compile<TIn, TOut> #1620

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 21 additions & 0 deletions BREAKING-CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2354,6 +2354,27 @@ two of them, the points are samples of a continuum and not the set
| `"floor(x)".Compile("x")`, and `ceil`, `round`, `max`, `min` | `UncompilableNodeException` | compiled |
| `"max(x, 1)".ToEntity().Differentiate("x")` | `derivative(max(x, 1), x)` | `(1 + sgn(x - 1)) / 2 provided not x - 1 = 0` |

### The compiled factorial is exact at whole numbers and finite to `170!`, and the compiled totient has no value at a fraction

Both compilers now evaluate the factorial, the gamma function and Euler's totient the same way,
and `Compile<TIn, TOut>` has a form for them where it had none
([#1607](https://github.com/asc-community/AngouriMath/issues/1607)). The compiled factorial was
Lanczos's approximation everywhere, and its power overflowed before the value did: from just past
`141!`, about `2e243`, it compiled to NaN, up to `170!`. At a whole number it was off from the
integer by a rounding, which a compilation over integers would truncate. It is now the exact
factorial rounded once, up to `170!`; infinite past it; NaN at the negative whole numbers, the
poles; and the same approximation elsewhere, taken in logarithms so that it overflows only where
the value does. The compiled totient truncated a fraction to the whole number below it. There it
now has no value, as the interpreter's has none.

| Input | Was (2.5.0) | Now |
|---|---|---|
| `"x!".Compile("x").Call(4)` | `23.999999999999996` | `24` |
| `"x!".Compile("x").Call(0)` | `0.9999999999999998` | `1` |
| `"x!".Compile("x").Call(142)`, and on to `170` | NaN | `2.695364137888163E+245` |
| `"phi(x)".Compile("x").Call(2.5)` | `1` | NaN |
| `"x!".ToEntity().Compile<double, double>("x")`, and `gamma`, `phi` | `UncompilableNodeException` | compiled |

### An equation the solver cannot invert is left unsolved, not answered with no roots

`x! = 6` was answered `{ }`, a claim that it has no roots, and it has 3. The solver isolates `x`
Expand Down
28 changes: 2 additions & 26 deletions Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -258,24 +258,7 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value
stack.Push(AngouriMath.Numerics.SpecialFunctions.Chi(stack.Pop()));
break;
case InstructionType.CALL_FACTORIAL:
// https://stackoverflow.com/a/15454784/5429648
const int g = 7;
static System.Numerics.Complex Gamma(System.Numerics.Complex z)
{
if (z.Real < 0.5) return System.Math.PI / (System.Numerics.Complex.Sin(System.Math.PI * z) * Gamma(1 - z));
else
{
z -= 1;

System.Numerics.Complex x = gammaCoeffs[0];
for (var i = 1; i < g + 2; i++)
x += gammaCoeffs[i] / (z + i);

var t = z + g + 0.5;
return System.Math.Sqrt(2 * System.Math.PI) * System.Numerics.Complex.Pow(t, z + 0.5) * System.Numerics.Complex.Exp(-t) * x;
}
}
stack.Push(Gamma(stack.Pop() + 1));
stack.Push(AngouriMath.Numerics.SpecialFunctions.Factorial(stack.Pop()));
break;
case InstructionType.CALL_SIGNUM:
stack.Push(stack.Pop().Signum());
Expand All @@ -284,21 +267,14 @@ static System.Numerics.Complex Gamma(System.Numerics.Complex z)
stack.Push(stack.Pop().Abs());
break;
case InstructionType.CALL_PHI:
var n = (long) stack.Pop().Real;
stack.Push(n.Phi());
stack.Push(AngouriMath.Numerics.SpecialFunctions.Phi(stack.Pop()));
break;
}
if (stack.Count != 1)
throw new AngouriBugException($"Unused values remain in the stack: {stack.Count} instead of 1");
return stack.Pop();
}

[ConstantField] static readonly double[] gammaCoeffs = {
0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
};

/// <summary>Might be useful for debug if a function works too slowly</summary>
public override string ToString() => string.Join(" \n| ", instructions);
}
Expand Down
2 changes: 2 additions & 0 deletions Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,8 @@ Sif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Si",
Cif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Ci", 1, newE.Type), newE),
Shif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Shi", 1, newE.Type), newE),
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),
Phif => Expression.Call(GetDef("Phi", 1, e.Type), e),

Notf => Expression.Not(e),

Expand Down
23 changes: 21 additions & 2 deletions Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -288,8 +288,17 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a)
public static int Min(int a, int b) => Math.Min(a, b);
public static BigInteger Min(BigInteger a, BigInteger b) => BigInteger.Min(a, b);

// The special functions, in double precision as the kernels define them.
// https://github.com/asc-community/AngouriMath/issues/1607
// Euler's totient: exact over the integers, and over Complex and double a value at a whole
// number and NaN at any other. https://github.com/asc-community/AngouriMath/issues/1607
public static System.Numerics.Complex Phi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
public static double Phi(double a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
public static float Phi(float a) => (float)AngouriMath.Numerics.SpecialFunctions.Phi(a);
public static long Phi(long a) => a.Phi();
public static int Phi(int a) => (int)((long)a).Phi();
public static BigInteger Phi(BigInteger a) => PeterO.Numbers.EInteger.FromBytes(a.ToByteArray(), littleEndian: true).Phi().ToBigInteger();

// The special functions, in double precision as the kernels define them, and the factorial,
// the gamma function one along. https://github.com/asc-community/AngouriMath/issues/1607
public static System.Numerics.Complex Erf(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erf(a);
public static double Erf(double a) => AngouriMath.Numerics.SpecialFunctions.Erf(a);
public static System.Numerics.Complex Erfc(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erfc(a);
Expand All @@ -308,6 +317,8 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a)
public static double Shi(double a) => AngouriMath.Numerics.SpecialFunctions.Shi(a);
public static System.Numerics.Complex Chi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Chi(a);
public static double Chi(double a) => AngouriMath.Numerics.SpecialFunctions.Chi(a);
public static System.Numerics.Complex Factorial(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Factorial(a);
public static double Factorial(double a) => AngouriMath.Numerics.SpecialFunctions.Factorial(a);

/// <summary>
/// True where <c>System.Numerics.Complex.IsNaN</c> is, spelled out because that overload is
Expand Down Expand Up @@ -459,6 +470,12 @@ public static bool IsNaN(System.Numerics.Complex a)
{ ("Round", 1, typeof(long)), Def<long>(a => Round(a)) },
{ ("Round", 1, typeof(int)), Def<int>(a => Round(a)) },
{ ("Round", 1, typeof(BigInteger)), Def<BigInteger>(a => Round(a)) },
{ ("Phi", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Phi(a)) },
{ ("Phi", 1, typeof(double)), Def<double>(a => Phi(a)) },
{ ("Phi", 1, typeof(float)), Def<float>(a => Phi(a)) },
{ ("Phi", 1, typeof(long)), Def<long>(a => Phi(a)) },
{ ("Phi", 1, typeof(int)), Def<int>(a => Phi(a)) },
{ ("Phi", 1, typeof(BigInteger)), Def<BigInteger>(a => Phi(a)) },
{ ("Erf", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Erf(a)) },
{ ("Erf", 1, typeof(double)), Def<double>(a => Erf(a)) },
{ ("Erfc", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Erfc(a)) },
Expand All @@ -477,6 +494,8 @@ public static bool IsNaN(System.Numerics.Complex a)
{ ("Shi", 1, typeof(double)), Def<double>(a => Shi(a)) },
{ ("Chi", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Chi(a)) },
{ ("Chi", 1, typeof(double)), Def<double>(a => Chi(a)) },
{ ("Factorial", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Factorial(a)) },
{ ("Factorial", 1, typeof(double)), Def<double>(a => Factorial(a)) },
};

// An expression tree is how C# spells "the MethodInfo of this method" without a string:
Expand Down
17 changes: 13 additions & 4 deletions Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,11 @@
// to one and forgotten in the other.
var unary = new[]{ "Sin", "Cos", "Tan", "Asin", "Acos", "Atan",
"Cot", "Sec", "Csc", "Acot", "Asec", "Acsc",
"Abs", "Sgn", "Floor", "Ceil", "Round" };
"Abs", "Sgn", "Floor", "Ceil", "Round", "Phi" };
var binary = new[]{ "Log", "Pow", "Max", "Min" };
// Defined for Complex and double alone: an integer or a float is widened to double first, as
// for the trigonometry.
var special = new[]{ "Erf", "Erfc", "Erfi", "Ei", "Li", "Si", "Ci", "Shi", "Chi" };
var special = new[]{ "Erf", "Erfc", "Erfi", "Ei", "Li", "Si", "Ci", "Shi", "Chi", "Factorial" };
var types = new[]{ "System.Numerics.Complex", "double", "float", "long", "int", "BigInteger" };
#>
//
Expand Down Expand Up @@ -191,8 +191,17 @@ namespace AngouriMath.Core.Compilation.IntoLinq
public static int Min(int a, int b) => Math.Min(a, b);
public static BigInteger Min(BigInteger a, BigInteger b) => BigInteger.Min(a, b);

// The special functions, in double precision as the kernels define them.
// https://github.com/asc-community/AngouriMath/issues/1607
// Euler's totient: exact over the integers, and over Complex and double a value at a whole
// number and NaN at any other. https://github.com/asc-community/AngouriMath/issues/1607
public static System.Numerics.Complex Phi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
public static double Phi(double a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
public static float Phi(float a) => (float)AngouriMath.Numerics.SpecialFunctions.Phi(a);
public static long Phi(long a) => a.Phi();
public static int Phi(int a) => (int)((long)a).Phi();
public static BigInteger Phi(BigInteger a) => PeterO.Numbers.EInteger.FromBytes(a.ToByteArray(), littleEndian: true).Phi().ToBigInteger();

// The special functions, in double precision as the kernels define them, and the factorial,
// the gamma function one along. https://github.com/asc-community/AngouriMath/issues/1607
<# foreach (var name in special) { #>
public static System.Numerics.Complex <#= name #>(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.<#= name #>(a);
public static double <#= name #>(double a) => AngouriMath.Numerics.SpecialFunctions.<#= name #>(a);
Expand Down
114 changes: 114 additions & 0 deletions Sources/AngouriMath/Numerics/SpecialFunctions.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -380,6 +380,120 @@ private static NumericsComplex HalfOfTheExponentialSeries(NumericsComplex z, boo
return sum;
}

/// <summary>
/// <c>x!</c>, the gamma function one along. At a whole number up to 170 it is the exact
/// factorial rounded once, and past 170, where a double cannot hold it, it is infinite. At a
/// negative whole number, one of the gamma function's poles, it is NaN, as the interpreter's
/// is. Elsewhere it is Lanczos's approximation with <c>g = 7</c> and nine terms, the one the
/// compiler has always used. Its <c>t^(z + 1/2) e^(-t)</c> is taken in logarithms, so that
/// it overflows only where the value does: the power alone passed the largest double just past
/// 141, where the factorial is about <c>2e243</c>, and the compiled factorial was NaN from
/// there to 170.
/// https://github.com/asc-community/AngouriMath/issues/1607
/// </summary>
public static NumericsComplex Factorial(NumericsComplex z)
{
if (z.Imaginary == 0)
return Factorial(z.Real);
if (IsNaN(z) || IsInfinite(z))
return NaN;
return Gamma(new NumericsComplex(z.Real + 1, z.Imaginary));
}

public static double Factorial(double x)
{
if (double.IsNaN(x) || double.IsNegativeInfinity(x))
return double.NaN;
if (double.IsPositiveInfinity(x))
return double.PositiveInfinity;
if (x == Math.Floor(x))
return x < 0 ? double.NaN : x <= 170 ? factorials[(int)x] : double.PositiveInfinity;
return Gamma(x + 1);
}

/// <summary><c>0!</c> to <c>170!</c>, each the exact integer rounded once to a double.</summary>
[AngouriMath.Core.ConstantField] private static readonly double[] factorials = ExactFactorials();

private static double[] ExactFactorials()
{
var table = new double[171];
var factorial = PeterO.Numbers.EInteger.One;
for (var n = 0; n <= 170; n++)
{
if (n > 0)
factorial *= n;
table[n] = PeterO.Numbers.EFloat.FromEInteger(factorial).ToDouble();
}
return table;
}

[AngouriMath.Core.ConstantField] private static readonly double[] lanczos = {
0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 };

private static double Gamma(double x)
{
// The reflection formula below 1/2.
if (x < 0.5)
return Math.PI / (Math.Sin(Math.PI * x) * Gamma(1 - x));
x -= 1;
var sum = lanczos[0];
for (var i = 1; i < 9; i++)
sum += lanczos[i] / (x + i);
var t = x + 7.5;
return Math.Exp((x + 0.5) * Math.Log(t) - t + Math.Log(Math.Sqrt(2 * Math.PI) * sum));
}

/// <summary>
/// The gamma function off the real line, as the exponential of its logarithm, so that it
/// overflows and underflows only where the value does. The reflection formula
/// <c>pi/(sin(pi z) gamma(1 - z))</c> is taken in logarithms too: far from the real line
/// <c>sin(pi z)</c> overflows where the value it divides is tiny, and it made a NaN.
/// </summary>
private static NumericsComplex Gamma(NumericsComplex z)
=> z.Real < 0.5
? ExpTimes(Math.Log(Math.PI) - LogSine(Math.PI * z) - LogGamma(1 - z), 1)
: ExpTimes(LogGamma(z), 1);

/// <summary>A logarithm of the gamma function for <c>Re z &gt;= 1/2</c>, by Lanczos's approximation.</summary>
private static NumericsComplex LogGamma(NumericsComplex z)
{
z -= 1;
NumericsComplex sum = lanczos[0];
for (var i = 1; i < 9; i++)
sum += lanczos[i] / (z + i);
var t = z + 7.5;
return (z + 0.5) * NumericsComplex.Log(t) - t + NumericsComplex.Log(Math.Sqrt(2 * Math.PI) * sum);
}

/// <summary>
/// A logarithm of <c>sin w</c>. Past <c>|Im w| = 20</c> one of <c>e^(i w)</c> and
/// <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
/// other. Its logarithm is then written out, where <c>sin w</c> itself would overflow.
/// </summary>
private static NumericsComplex LogSine(NumericsComplex w)
{
if (Math.Abs(w.Imaginary) < 20)
return NumericsComplex.Log(NumericsComplex.Sin(w));
// (i/2) e^(-i w) above the real line, and -(i/2) e^(i w) below it.
return w.Imaginary > 0
? new NumericsComplex(w.Imaginary - Math.Log(2), Math.PI / 2 - w.Real)
: new NumericsComplex(-w.Imaginary - Math.Log(2), w.Real - Math.PI / 2);
}

/// <summary>
/// Euler's totient at a whole number, 0 at one that is not positive, as the interpreter has
/// it, and NaN at any other number, where the interpreter's is NaN too. The compiler
/// truncated a fraction to the whole number below it, so <c>phi(2.5)</c> was <c>phi(2)</c>.
/// Past <c>long</c>'s range it is NaN: the interpreter's totient cannot go there either.
/// https://github.com/asc-community/AngouriMath/issues/1607
/// </summary>
public static NumericsComplex Phi(NumericsComplex z) => z.Imaginary == 0 ? Phi(z.Real) : NaN;

public static double Phi(double x)
=> x == Math.Floor(x) && x >= long.MinValue && x < long.MaxValue ? ((long)x).Phi() : double.NaN;

/// <summary>
/// Half of <paramref name="value"/>, part by part: dividing a complex number by 2 multiplies it
/// by the complex number 2, and an infinite part times the other's zero makes a NaN.
Expand Down
Loading
Loading