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"));