From c32ee42fa2c16282b29d22a024f41a42e9e26e45 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 30 Oct 2026 02:07:29 +0000 Subject: [PATCH] The number-theoretic functions compile mobius, divisorsigma, primenu, primeomega, valuation, prime and primepi compile in FastExpression and in the LINQ compiler, over Complex and over double. prime and valuation were refused with UncompilableNodeException; the other five are new since 2.6.0. Each is computed by the routines the interpreter uses: Primes.Factorization for the functions of a factorisation, IsPrime for the valuation's prime, GetPrime for the n-th prime and CountUpTo for the prime count. So a compiled value is exact at a whole argument and NaN at any other, as evaluation gives, and a number too large to factorise or a prime past the millionth is NaN where the interpreter leaves it as written. The test compares both compilers with evaluation at 68 arguments for each function of a factorisation, including 1099511627776 = 2^40, 600851475143 and the non-whole 5/2 and 1 + i; the divisor sigma at six orders, among them -1, 1/2 and i; the valuation at seven numbers and six bases, two of them not prime; the n-th prime up to the 10000th; and the prime count up to 100000. prime(x) at 10: UncompilableNodeException (2.6.0), 29 now. valuation(x, 2) at 12: UncompilableNodeException (2.6.0), 2 now. Part of #2195. Co-Authored-By: Claude Opus 5.5 (1M context) --- BREAKING-CHANGES.md | 15 +++ Sources/.editorconfig | 6 + .../Functions/Compilation/IntoFE/Compiler.cs | 81 +++++++++++++ .../Compilation/IntoFE/FastExpression.cs | 28 +++++ .../IntoLinq/CompilationProtocol.cs | 11 ++ .../Compilation/IntoLinq/MathAllMethods.cs | 31 +++++ .../Compilation/IntoLinq/MathAllMethods.tt | 24 ++++ .../Numerics/NumberTheoreticFunctions.cs | 111 ++++++++++++++++++ Sources/Tests/AotSmokeTest/Program.cs | 5 + .../Common/NumberTheoreticCompilationTest.cs | 83 +++++++++++++ .../UnitTests/Core/UserInvalidExceptions.cs | 3 +- 11 files changed, 396 insertions(+), 2 deletions(-) create mode 100644 Sources/AngouriMath/Numerics/NumberTheoreticFunctions.cs create mode 100644 Sources/Tests/UnitTests/Common/NumberTheoreticCompilationTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index dfe88344c..35f5e6f2e 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -337,6 +337,21 @@ factorial has no value off the whole numbers or at a negative even one | `"binomial(x, 2)".ToEntity().Compile("x").Call(5)` | `UncompilableNodeException` | `(10, 0)` | | `"binomial(x, 1/2)".ToEntity().Compile("x")(3)` | `UncompilableNodeException` | `2.037183271576257` | +### The number-theoretic functions compile + +**Not silent: what was refused is compiled.** `prime` and `valuation` were refused by both compilers with +`UncompilableNodeException`, and `mobius`, `divisorsigma`, `primenu`, `primeomega` and `primepi`, new since 2.6.0, +compile from the start. Each compiles by the routines the interpreter uses, exactly at a whole argument and to `NaN` at +any other, as evaluation gives. A number too large to factorise, and a prime past the millionth, are `NaN` where the +interpreter leaves them as written ([#2195](https://github.com/asc-community/AngouriMath/issues/2195)). + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"prime(x)".ToEntity().Compile("x").Call(10)` | `UncompilableNodeException` | `(29, 0)` | +| `"prime(x)".ToEntity().Compile("x")(2.5)` | `UncompilableNodeException` | `NaN` | +| `"valuation(x, 2)".ToEntity().Compile("x").Call(12)` | `UncompilableNodeException` | `(2, 0)` | +| `"valuation(x, 2)".ToEntity().Compile("x")(0)` | `UncompilableNodeException` | `Infinity` | + ### `mod` compiles over `Complex` **Not silent: what was refused is compiled.** `Compile` refused `mod` with `UncompilableNodeException`, diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 14650bf62..c04f68364 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -375,6 +375,12 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Common/PolygammaCompilationTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[AngouriMath/Numerics/NumberTheoreticFunctions.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[Tests/UnitTests/Common/NumberTheoreticCompilationTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [Tests/UnitTests/Common/CompileMatrixTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs index 18e5c3970..ac2f221f0 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs @@ -653,6 +653,87 @@ private protected override void CompileNode(Compiler compiler) } } + // The number-theoretic functions, exact at a whole argument and NaN at any other. + // https://github.com/asc-community/AngouriMath/issues/2195 + public partial record Mobiusf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_MOBIUS)); + } + } + + public partial record PrimeNuf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_PRIME_NU)); + } + } + + public partial record PrimeOmegaf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_PRIME_OMEGA)); + } + } + + public partial record Primef + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_PRIME)); + } + } + + public partial record PrimePif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_PRIME_PI)); + } + } + + public partial record DivisorSigmaf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + Order.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_DIVISOR_SIGMA)); + } + } + + public partial record Valuationf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Prime.InnerCompile(compiler); + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_VALUATION)); + } + } + public partial record Maxf { private protected override bool CompilesItself => true; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs index 3f54d16c6..cab923c8f 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs @@ -57,6 +57,11 @@ internal enum InstructionType CALL_NOT, CALL_DOUBLE_FACTORIAL, CALL_CATALAN_NUMBER, + CALL_MOBIUS, + CALL_PRIME_NU, + CALL_PRIME_OMEGA, + CALL_PRIME, + CALL_PRIME_PI, // 2-arg functions CALL_SUM = 100, @@ -85,6 +90,8 @@ internal enum InstructionType CALL_LAMBERTW, CALL_POLYGAMMA, CALL_HARMONIC_NUMBER, + CALL_DIVISOR_SIGMA, + CALL_VALUATION, // 3-arg: a value where its predicate holds, and the one below it elsewhere. CALL_SELECT = 150, @@ -368,6 +375,27 @@ 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_MOBIUS: + stack.Push(AngouriMath.Numerics.NumberTheoreticFunctions.Mobius(stack.Pop())); + break; + case InstructionType.CALL_PRIME_NU: + stack.Push(AngouriMath.Numerics.NumberTheoreticFunctions.PrimeNu(stack.Pop())); + break; + case InstructionType.CALL_PRIME_OMEGA: + stack.Push(AngouriMath.Numerics.NumberTheoreticFunctions.PrimeOmega(stack.Pop())); + break; + case InstructionType.CALL_PRIME: + stack.Push(AngouriMath.Numerics.NumberTheoreticFunctions.Prime(stack.Pop())); + break; + case InstructionType.CALL_PRIME_PI: + stack.Push(AngouriMath.Numerics.NumberTheoreticFunctions.PrimePi(stack.Pop())); + break; + case InstructionType.CALL_DIVISOR_SIGMA: + stack.Push(AngouriMath.Numerics.NumberTheoreticFunctions.DivisorSigma(stack.Pop(), stack.Pop())); + break; + case InstructionType.CALL_VALUATION: + stack.Push(AngouriMath.Numerics.NumberTheoreticFunctions.Valuation(stack.Pop(), stack.Pop())); + break; case InstructionType.CALL_CATALAN_NUMBER: stack.Push(AngouriMath.Numerics.CombinatorialFunctions.CatalanNumber(stack.Pop())); break; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs index e7d4c8f89..f5dbb737c 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs @@ -158,6 +158,13 @@ Chif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Chi", 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), + // The number-theoretic functions, exact at a whole argument and NaN at any other. + // https://github.com/asc-community/AngouriMath/issues/2195 + Mobiusf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Mobius", 1, newE.Type), newE), + PrimeNuf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("PrimeNu", 1, newE.Type), newE), + PrimeOmegaf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("PrimeOmega", 1, newE.Type), newE), + Primef when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Prime", 1, newE.Type), newE), + PrimePif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("PrimePi", 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 @@ -249,6 +256,10 @@ Polygammaf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDou => Expression.Call(GetDef("Polygamma", 2, newLeft.Type), newLeft, newRight), HarmonicNumberf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(right) is var newRight => Expression.Call(GetDef("HarmonicNumber", 2, newLeft.Type), newLeft, newRight), + DivisorSigmaf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(right) is var newRight + => Expression.Call(GetDef("DivisorSigma", 2, newLeft.Type), newLeft, newRight), + Valuationf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(right) is var newRight + => Expression.Call(GetDef("Valuation", 2, newLeft.Type), newLeft, newRight), Maxf => Expression.Call(GetDef("Max", 2, left.Type), left, right), Minf => Expression.Call(GetDef("Min", 2, left.Type), left, right), diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs index 5fd183372..4af291eac 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs @@ -356,6 +356,23 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a) public static System.Numerics.Complex HarmonicNumber(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); public static double HarmonicNumber(double a, double b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); + // The functions of a factorisation, the valuation, the n-th prime and the prime-counting function, by the + // interpreter's routines. https://github.com/asc-community/AngouriMath/issues/2195 + public static System.Numerics.Complex Mobius(System.Numerics.Complex a) => AngouriMath.Numerics.NumberTheoreticFunctions.Mobius(a); + public static double Mobius(double a) => AngouriMath.Numerics.NumberTheoreticFunctions.Mobius(a); + public static System.Numerics.Complex PrimeNu(System.Numerics.Complex a) => AngouriMath.Numerics.NumberTheoreticFunctions.PrimeNu(a); + public static double PrimeNu(double a) => AngouriMath.Numerics.NumberTheoreticFunctions.PrimeNu(a); + public static System.Numerics.Complex PrimeOmega(System.Numerics.Complex a) => AngouriMath.Numerics.NumberTheoreticFunctions.PrimeOmega(a); + public static double PrimeOmega(double a) => AngouriMath.Numerics.NumberTheoreticFunctions.PrimeOmega(a); + public static System.Numerics.Complex Prime(System.Numerics.Complex a) => AngouriMath.Numerics.NumberTheoreticFunctions.Prime(a); + public static double Prime(double a) => AngouriMath.Numerics.NumberTheoreticFunctions.Prime(a); + public static System.Numerics.Complex PrimePi(System.Numerics.Complex a) => AngouriMath.Numerics.NumberTheoreticFunctions.PrimePi(a); + public static double PrimePi(double a) => AngouriMath.Numerics.NumberTheoreticFunctions.PrimePi(a); + public static System.Numerics.Complex DivisorSigma(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.NumberTheoreticFunctions.DivisorSigma(a, b); + public static double DivisorSigma(double a, double b) => AngouriMath.Numerics.NumberTheoreticFunctions.DivisorSigma(a, b); + public static System.Numerics.Complex Valuation(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.NumberTheoreticFunctions.Valuation(a, b); + public static double Valuation(double a, double b) => AngouriMath.Numerics.NumberTheoreticFunctions.Valuation(a, b); + /// /// True where System.Numerics.Complex.IsNaN is, spelled out because that overload is /// .NET 7 and later and this assembly also targets netstandard2.0. A complex number is NaN @@ -556,6 +573,20 @@ public static bool IsNaN(System.Numerics.Complex a) { ("LambertW", 2, typeof(double)), Def((a, b) => LambertW(a, b)) }, { ("Polygamma", 2, typeof(double)), Def((a, b) => Polygamma(a, b)) }, { ("HarmonicNumber", 2, typeof(double)), Def((a, b) => HarmonicNumber(a, b)) }, + { ("Mobius", 1, typeof(System.Numerics.Complex)), Def(a => Mobius(a)) }, + { ("Mobius", 1, typeof(double)), Def(a => Mobius(a)) }, + { ("PrimeNu", 1, typeof(System.Numerics.Complex)), Def(a => PrimeNu(a)) }, + { ("PrimeNu", 1, typeof(double)), Def(a => PrimeNu(a)) }, + { ("PrimeOmega", 1, typeof(System.Numerics.Complex)), Def(a => PrimeOmega(a)) }, + { ("PrimeOmega", 1, typeof(double)), Def(a => PrimeOmega(a)) }, + { ("Prime", 1, typeof(System.Numerics.Complex)), Def(a => Prime(a)) }, + { ("Prime", 1, typeof(double)), Def(a => Prime(a)) }, + { ("PrimePi", 1, typeof(System.Numerics.Complex)), Def(a => PrimePi(a)) }, + { ("PrimePi", 1, typeof(double)), Def(a => PrimePi(a)) }, + { ("DivisorSigma", 2, typeof(System.Numerics.Complex)), Def((a, b) => DivisorSigma(a, b)) }, + { ("DivisorSigma", 2, typeof(double)), Def((a, b) => DivisorSigma(a, b)) }, + { ("Valuation", 2, typeof(System.Numerics.Complex)), Def((a, b) => Valuation(a, b)) }, + { ("Valuation", 2, typeof(double)), Def((a, b) => Valuation(a, b)) }, }; // An expression tree is how C# spells "the MethodInfo of this method" without a string: diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt index f6e87ece9..d0022b23c 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt @@ -19,6 +19,9 @@ // https://github.com/asc-community/AngouriMath/issues/2195 var combinatorial = new[]{ "DoubleFactorial", "CatalanNumber" }; var combinatorialBinary = new[]{ "Binomial", "FallingFactorial", "RisingFactorial" }; + // The number-theoretic functions, exact at a whole argument and NaN at any other. + var numberTheoretic = new[]{ "Mobius", "PrimeNu", "PrimeOmega", "Prime", "PrimePi" }; + var numberTheoreticBinary = new[]{ "DivisorSigma", "Valuation" }; var types = new[]{ "System.Numerics.Complex", "double", "float", "long", "int", "BigInteger" }; #> // @@ -244,6 +247,17 @@ namespace AngouriMath.Core.Compilation.IntoLinq public static System.Numerics.Complex HarmonicNumber(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); public static double HarmonicNumber(double a, double b) => AngouriMath.Numerics.PolygammaFunction.HarmonicNumber(a, b); + // The functions of a factorisation, the valuation, the n-th prime and the prime-counting function, by the + // interpreter's routines. https://github.com/asc-community/AngouriMath/issues/2195 +<# foreach (var name in numberTheoretic) { #> + public static System.Numerics.Complex <#= name #>(System.Numerics.Complex a) => AngouriMath.Numerics.NumberTheoreticFunctions.<#= name #>(a); + public static double <#= name #>(double a) => AngouriMath.Numerics.NumberTheoreticFunctions.<#= name #>(a); +<# } #> +<# foreach (var name in numberTheoreticBinary) { #> + public static System.Numerics.Complex <#= name #>(System.Numerics.Complex a, System.Numerics.Complex b) => AngouriMath.Numerics.NumberTheoreticFunctions.<#= name #>(a, b); + public static double <#= name #>(double a, double b) => AngouriMath.Numerics.NumberTheoreticFunctions.<#= name #>(a, b); +<# } #> + /// /// True where System.Numerics.Complex.IsNaN is, spelled out because that overload is /// .NET 7 and later and this assembly also targets netstandard2.0. A complex number is NaN @@ -302,6 +316,16 @@ namespace AngouriMath.Core.Compilation.IntoLinq { ("LambertW", 2, typeof(<#= type #>)), Def<<#= type #>>((a, b) => LambertW(a, b)) }, { ("Polygamma", 2, typeof(<#= type #>)), Def<<#= type #>>((a, b) => Polygamma(a, b)) }, { ("HarmonicNumber", 2, typeof(<#= type #>)), Def<<#= type #>>((a, b) => HarmonicNumber(a, b)) }, +<# } #> +<# foreach (var name in numberTheoretic) { #> +<# foreach (var type in new[]{ "System.Numerics.Complex", "double" }) { #> + { ("<#= name #>", 1, typeof(<#= type #>)), Def<<#= type #>>(a => <#= name #>(a)) }, +<# } #> +<# } #> +<# foreach (var name in numberTheoreticBinary) { #> +<# foreach (var type in new[]{ "System.Numerics.Complex", "double" }) { #> + { ("<#= name #>", 2, typeof(<#= type #>)), Def<<#= type #>>((a, b) => <#= name #>(a, b)) }, +<# } #> <# } #> }; diff --git a/Sources/AngouriMath/Numerics/NumberTheoreticFunctions.cs b/Sources/AngouriMath/Numerics/NumberTheoreticFunctions.cs new file mode 100644 index 000000000..7b659ebc0 --- /dev/null +++ b/Sources/AngouriMath/Numerics/NumberTheoreticFunctions.cs @@ -0,0 +1,111 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Functions; +using PeterO.Numbers; +using NumericsComplex = System.Numerics.Complex; + +namespace AngouriMath.Numerics +{ + /// + /// The functions of a factorisation, the p-adic valuation, the n-th prime and the prime-counting + /// function, for the compilers, by the routines the interpreter uses: exact at a whole argument, and + /// NaN at any other, where none of them has a value, and where the interpreter leaves one as + /// written, a number too large to factorise or a prime too far to build towards. + /// https://github.com/asc-community/AngouriMath/issues/2195 + /// + internal static class NumberTheoreticFunctions + { + [AngouriMath.Core.ConstantField] private static readonly NumericsComplex NaN = new(double.NaN, 0); + + /// Whether is a whole number a holds exactly. + private static bool IsWhole(NumericsComplex z, out long whole) + { + whole = 0; + if (z.Imaginary != 0 || double.IsNaN(z.Real) || double.IsInfinity(z.Real) + || Math.Floor(z.Real) != z.Real || Math.Abs(z.Real) > 9e15) + return false; + whole = (long)z.Real; + return true; + } + + // The factorisation of a positive whole number, or null where there is none or it is not found. + private static System.Collections.Generic.List<(EInteger Prime, int Power)>? Factors(NumericsComplex n) + => IsWhole(n, out var whole) && whole > 0 ? Primes.Factorization(EInteger.FromInt64(whole)) : null; + + /// The Möbius function: 0 where a square divides n, and (-1)^k for k distinct primes. + public static NumericsComplex Mobius(NumericsComplex n) + => Factors(n) is { } factors ? (factors.Exists(factor => factor.Power > 1) ? 0 : factors.Count % 2 == 0 ? 1 : -1) : NaN; + + /// The number of distinct primes dividing n. + public static NumericsComplex PrimeNu(NumericsComplex n) => Factors(n) is { } factors ? factors.Count : NaN; + + /// The number of primes dividing n, with multiplicity. + public static NumericsComplex PrimeOmega(NumericsComplex n) + { + if (Factors(n) is not { } factors) + return NaN; + var count = 0; + foreach (var factor in factors) + count += factor.Power; + return count; + } + + /// The sum of the k-th powers of the divisors of n, by its product over the prime powers. + public static NumericsComplex DivisorSigma(NumericsComplex k, NumericsComplex n) + { + if (Factors(n) is not { } factors || double.IsNaN(k.Real) || double.IsNaN(k.Imaginary)) + return NaN; + NumericsComplex product = 1; + foreach (var (prime, power) in factors) + { + NumericsComplex sum = 0; + var p = (double)prime.ToInt64Checked(); + for (var j = 0; j <= power; j++) + sum += NumericsComplex.Pow(Math.Pow(p, j), k); + product *= sum; + } + return product; + } + + /// The p-adic valuation of a whole n at a prime p: +oo at 0, and no value for a p not prime. + public static NumericsComplex Valuation(NumericsComplex n, NumericsComplex p) + { + if (!IsWhole(n, out var whole) || !IsWhole(p, out var prime) || Primes.IsPrime(EInteger.FromInt64(prime)) is not true) + return NaN; + if (whole == 0) + return double.PositiveInfinity; + var count = 0; + for (var rest = Math.Abs(whole); rest % prime == 0; rest /= prime) + count++; + return count; + } + + /// The n-th prime, for a whole n from 1 to a million, as far as the table is built. + public static NumericsComplex Prime(NumericsComplex n) + => IsWhole(n, out var index) && index >= 1 && index <= 1000000 ? Primes.GetPrime((int)index - 1).EInteger.ToInt64Checked() : NaN; + + /// The number of primes up to a real x, as far as the sieve counts. + public static NumericsComplex PrimePi(NumericsComplex x) + { + if (x.Imaginary != 0 || double.IsNaN(x.Real) || double.IsInfinity(x.Real) || Math.Abs(x.Real) > 9e15) + return NaN; + return Primes.CountUpTo(EInteger.FromInt64((long)Math.Floor(x.Real))) is { } count ? count : NaN; + } + + public static double Mobius(double n) => Real(Mobius(new NumericsComplex(n, 0))); + public static double PrimeNu(double n) => Real(PrimeNu(new NumericsComplex(n, 0))); + public static double PrimeOmega(double n) => Real(PrimeOmega(new NumericsComplex(n, 0))); + public static double DivisorSigma(double k, double n) => Real(DivisorSigma(new NumericsComplex(k, 0), new NumericsComplex(n, 0))); + public static double Valuation(double n, double p) => Real(Valuation(new NumericsComplex(n, 0), new NumericsComplex(p, 0))); + public static double Prime(double n) => Real(Prime(new NumericsComplex(n, 0))); + public static double PrimePi(double x) => Real(PrimePi(new NumericsComplex(x, 0))); + + private static double Real(NumericsComplex z) => z.Imaginary == 0 ? z.Real : double.NaN; + } +} diff --git a/Sources/Tests/AotSmokeTest/Program.cs b/Sources/Tests/AotSmokeTest/Program.cs index 65000789b..e1df0093f 100644 --- a/Sources/Tests/AotSmokeTest/Program.cs +++ b/Sources/Tests/AotSmokeTest/Program.cs @@ -315,6 +315,11 @@ private static void Compilation() 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)}"); + // mobius(30) = -1, prime(30) = 113, primepi(30) = 10; divisorsigma(1, 12) = 28, valuation(12, 2) = 2. + var fastNumberTheory = "mobius(x) + prime(x) + primepi(x)".Compile(x); + var linqNumberTheory = "divisorsigma(1, x) + valuation(x, 2)".Compile(x); + Check("compiled mobius(x) + prime(x) + primepi(x) at 30, divisorsigma(1, x) + valuation(x, 2) at 12", + fastNumberTheory.Call(30) == 122 && linqNumberTheory(12) == 30, $"{fastNumberTheory.Call(30)}, {linqNumberTheory(12)}"); // 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); diff --git a/Sources/Tests/UnitTests/Common/NumberTheoreticCompilationTest.cs b/Sources/Tests/UnitTests/Common/NumberTheoreticCompilationTest.cs new file mode 100644 index 000000000..b90f5e9af --- /dev/null +++ b/Sources/Tests/UnitTests/Common/NumberTheoreticCompilationTest.cs @@ -0,0 +1,83 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System.Linq; +using System.Numerics; +using AngouriMath; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Common +{ + /// + /// The functions of a factorisation, the valuation, the n-th prime and the prime-counting function + /// compile, in FastExpression and in the LINQ compiler, to what evaluation gives: exactly at a whole + /// argument, and NaN at any other. Each was refused. + /// #2195 + /// + [Trait("Area", "Common")] + public sealed class NumberTheoreticCompilationTest + { + private static readonly string[] wholes = + new[] { "-3", "0", "97", "1000003", "1099511627776", "600851475143", "5/2", "1 + i" } + .Concat(Enumerable.Range(1, 60).Select(n => n.ToString())).ToArray(); + + private static Complex Numeric(string value) => ((Entity.Number.Complex)value.ToEntity().Evaled).ToNumerics(); + + // Both NaN, or equal to a part in 1e12. + private static bool Agree(Complex expected, Complex actual) + => double.IsNaN(expected.Real) || double.IsNaN(expected.Imaginary) + ? double.IsNaN(actual.Real) || double.IsNaN(actual.Imaginary) + : double.IsInfinity(expected.Real) ? expected == actual + : Complex.Abs(expected - actual) <= 1e-12 * System.Math.Max(1, Complex.Abs(expected)); + + private static void Compare(string expression, string[] xs, string[] ys) + { + var fast = expression.ToEntity().Compile("x", "y"); + var linq = expression.ToEntity().Compile("x", "y"); + var real = expression.ToEntity().Compile("x", "y"); + var compared = 0; + foreach (var a in xs) + foreach (var b in ys) + { + if (expression.ToEntity().Substitute("x", a).Substitute("y", b).Evaled is not Entity.Number.Complex evaluated) + continue; + var expected = evaluated.ToNumerics(); + var (x, y) = (Numeric(a), Numeric(b)); + Assert.True(Agree(expected, fast.Call(x, y)), $"{expression} at ({a}, {b}): {expected} evaluated, {fast.Call(x, y)} by FastExpression"); + Assert.True(Agree(expected, linq(x, y)), $"{expression} at ({a}, {b}): {expected} evaluated, {linq(x, y)} by LINQ"); + if (x.Imaginary == 0 && y.Imaginary == 0 && expected.Imaginary == 0) + Assert.True(Agree(expected, real(x.Real, y.Real)), $"{expression} at ({a}, {b}): {expected} evaluated, {real(x.Real, y.Real)} over double"); + compared++; + } + Assert.True(compared > 0, $"nothing of {expression} evaluated"); + } + + [Theory] + [InlineData("mobius(x) + 0 y")] + [InlineData("primenu(x) + 0 y")] + [InlineData("primeomega(x) + 0 y")] + public void AFunctionOfAFactorisationIsWhatEvaluationSays(string expression) + => Compare(expression, wholes, new[] { "0" }); + + [Fact] + public void TheDivisorSigmaIsWhatEvaluationSays() + => Compare("divisorsigma(y, x)", wholes, new[] { "0", "1", "2", "-1", "1/2", "i" }); + + [Fact] + public void TheValuationIsWhatEvaluationSays() + => Compare("valuation(x, y)", new[] { "0", "12", "-12", "1024", "7", "5/2", "600851475143" }, new[] { "2", "3", "4", "7", "1/2", "71" }); + + [Fact] + public void TheNthPrimeIsWhatEvaluationSays() + => Compare("prime(x) + 0 y", new[] { "1", "2", "10", "1000", "10000", "0", "-2", "5/2", "i" }, new[] { "0" }); + + [Fact] + public void ThePrimeCountIsWhatEvaluationSays() + => Compare("primepi(x) + 0 y", new[] { "-3", "0", "2", "5/2", "100", "100000", "i" }, new[] { "0" }); + } +} diff --git a/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs b/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs index 687c1f163..e06e22ac3 100644 --- a/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs +++ b/Sources/Tests/UnitTests/Core/UserInvalidExceptions.cs @@ -57,13 +57,12 @@ [Fact] public void UncompilableNode4() => // of these were nodes 2.0 added and never taught the compiler; floor, ceil and round have // been taught since (https://github.com/asc-community/AngouriMath/issues/1603), phi, // gamma and the factorial (https://github.com/asc-community/AngouriMath/issues/1607), and - // the binomial coefficient and the polygamma functions + // the binomial coefficient, the polygamma functions and the n-th prime // (https://github.com/asc-community/AngouriMath/issues/2195). [Theory] [InlineData("gamma(2, x)")] [InlineData("gcd(x, 6)")] [InlineData("lcm(x, 6)")] - [InlineData("prime(x)")] public void ANodeWithNoCompiledFormSaysSoRatherThanAskingForABugReport(string expression) => Assert.Throws(() => expression.ToEntity().Compile("x"));