From 7290956911a88f4c9f97d43b34d239cbf2aaf508 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 30 Sep 2026 07:14:06 +0000 Subject: [PATCH 1/2] The special functions compile, in double precision erf, erfc, erfi, Ei, li, Si, Ci, Shi and Chi had no compiled form. So Compile and Compile threw UncompilableNodeException on them, and Solve's numerical search, which compiles, declined any equation with one in it. The interpreter's kernels work in arbitrary precision and take 0.7 to 5.5 ms a call at 20 digits. Newton's method makes about 3,000 calls an equation, so it cannot use them. SpecialFunctions gives each a double-precision routine, and both compilers call it: through new opcodes under Compile, and through definitions for Complex and double under Compile, where an integer or a float is widened to double first, as for the trigonometry. - The error functions come from the Faddeeva function, by Weideman's rational approximation with 40 terms, and from erf's Taylor series near the origin. - Ei is summed as its series where the terms cancel little, is E1's continued fraction elsewhere out to |z| = 40, and is its asymptotic series beyond. - li, Si, Ci, Shi and Chi follow from Ei and from their own series, by the relations the kernels state. Each keeps the kernel's branch cuts and its values on the axes. Over doubles, a value off the real line is NaN, as Math.Sqrt of a negative number is. Measured against the kernels at 30 digits on about ten thousand points, through the compiled delegates: the integrals agree to within 6e-14 of the value, and the error functions to within 3e-13, which is the rounding of z^2 in e^(-z^2) at large |z|. They take 0.5 to 2.4 us a call. Solve answers erf(x) = x/2 + 1/7 with its three real roots, each checked against mpmath, and more off the real line. An equation that isolates a special function, erf(x) = 1/2, is still left unsolved: the solver finds no inverse to write, and so never reaches the search. Part of #1607. Co-Authored-By: Claude Opus 5.5 --- BREAKING-CHANGES.md | 15 +- Sources/.editorconfig | 3 + .../Functions/Compilation/IntoFE/Compiler.cs | 101 +++++ .../Compilation/IntoFE/FastExpression.cs | 38 ++ .../IntoLinq/CompilationProtocol.cs | 12 + .../Compilation/IntoLinq/MathAllMethods.cs | 39 ++ .../Compilation/IntoLinq/MathAllMethods.tt | 15 + .../Functions/Compilation/SpecialFunctions.cs | 405 ++++++++++++++++++ .../Algebra/SolveTest/NewtonFallbackTest.cs | 36 +- .../SpecialFunctionCompilationTest.cs | 118 +++++ 10 files changed, 772 insertions(+), 10 deletions(-) create mode 100644 Sources/AngouriMath/Functions/Compilation/SpecialFunctions.cs create mode 100644 Sources/Tests/UnitTests/Convenience/SpecialFunctionCompilationTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index bce7cee9a..88a5535cc 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -2316,9 +2316,9 @@ resort, Newton's method, compiles the expression and its derivative. They compil not ordered; - `max` and `min` differentiate where their arguments are real, as `|f|` does. -The search still declines what has no compiled form, so `erf(x) = x/2 + 1/7` is left unsolved. It -now also declines a set of roots that are not isolated: where the equation holds midway between two -of them, the points are samples of a continuum and not the set +The search still declines what has no compiled form, so `binomial(x, 3/2) = x/2` is left unsolved. +It now also declines a set of roots that are not isolated: where the equation holds midway between +two of them, the points are samples of a continuum and not the set ([#1603](https://github.com/asc-community/AngouriMath/issues/1603), [#1420](https://github.com/asc-community/AngouriMath/issues/1420)). @@ -2387,7 +2387,10 @@ e^(-t^2) dt`, its complement `erfc(z) = 1 - erf(z)` and the imaginary error func `erfi(z) = -i erf(i z)` are nodes (`MathS.Erf`, `MathS.Erfc`, `MathS.Erfi`), the first of #1501's special functions. Each is differentiated, printed in LaTeX as `\operatorname{erf}` and so on, and evaluated to the working precision anywhere in the complex plane (`Entity.Number.Erf`, `Erfc`, -`Erfi`). An equation in one of them is left unsolved rather than answered with no roots. `erf` +`Erfi`). They compile, through `Compile` and `Compile`, in double precision +([#1607](https://github.com/asc-community/AngouriMath/issues/1607)). So an equation the solver's +numerical search reaches, `erf(x) = x/2 + 1/7`, is answered with the roots it finds, while one that +isolates `erf`, `erf(x) = 1/2`, is left unsolved rather than answered with no roots. `erf` was refused by name, but `erfc` and `erfi` were not, so each was read as a product with an undeclared variable. @@ -2398,6 +2401,8 @@ undeclared variable. | `"erfi(x)"` | `erfi * x` | `erfi(x)` | | `"erfc(1)".Simplify()` | `erfc` | `erfc(1)`, and `0.157299207050285…` evaluated | | `"erf(x)".Differentiate("x")` | `UnrecognizedFunctionParseException` | `2 / sqrt(pi) * e ^ (-x ^ 2)` | +| `"erf(x)".Compile("x")`, and `erfc`, `erfi` | `UnrecognizedFunctionParseException`; for `erfc` and `erfi`, `UncompilableNodeException`: not among the variables | compiled | +| `"erf(x) = x/2 + 1/7".ToEntity().Solve("x")` | `UnrecognizedFunctionParseException` | its roots near `-2.283`, `0.235` and `1.679`, and more off the real line | ### `Ei` and `li` are functions @@ -2417,6 +2422,7 @@ with an undeclared variable | `"li(2)".Simplify()` | `2 * li` | `li(2)`, and `1.045163780117492…` evaluated | | `"Ei(1)".Simplify()` | `Ei` | `Ei(1)`, and `1.895117816355936…` evaluated | | `"Ei(x)".Differentiate("x")` | `Ei` | `e ^ x / x` | +| `"Ei(x)".Compile("x")`, and `li` | `UncompilableNodeException`: `Ei` is not among the variables | compiled ([#1607](https://github.com/asc-community/AngouriMath/issues/1607)) | ### The exponential and logarithmic integrals are integrated to @@ -2455,6 +2461,7 @@ was refused, so each was read as a product with an undeclared variable | `"Ci(2)".Simplify()` | `2 * Ci` | `Ci(2)`, and `0.4229808287748649…` evaluated | | `"Chi(1)".Simplify()` | `Chi` | `Chi(1)`, and `0.8378669409802082…` evaluated | | `"Si(x)".Differentiate("x")` | `Si` | `sin(x) / x` | +| `"Si(x)".Compile("x")`, and `Ci`, `Shi`, `Chi` | `UncompilableNodeException`: `Si` is not among the variables | compiled ([#1607](https://github.com/asc-community/AngouriMath/issues/1607)) | ### A quadratic below a sine or a cosine is split over its two roots diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 2d7db79f0..7667f5a13 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -342,6 +342,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [AngouriMath/Core/Entity/Continuous/Number/ErrorFunctions.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/Functions/Compilation/SpecialFunctions.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/Core/Entity/Continuous/Entity.Continuous.ErrorFunctions.Classes.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 5fb13510d..38db33971 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs @@ -442,6 +442,107 @@ private protected override void CompileNode(Compiler compiler) compiler.Instructions.Add(new(InstructionType.CALL_PHI)); } } + + // The special functions, in double precision as the kernels define them. + // https://github.com/asc-community/AngouriMath/issues/1607 + public partial record Erff + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_ERF)); + } + } + + public partial record Erfcf + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_ERFC)); + } + } + + public partial record Erfif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_ERFI)); + } + } + + public partial record Eif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_EI)); + } + } + + public partial record Lif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_LI)); + } + } + + public partial record Sif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_SI)); + } + } + + public partial record Cif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_CI)); + } + } + + public partial record Shif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_SHI)); + } + } + + public partial record Chif + { + private protected override bool CompilesItself => true; + + private protected override void CompileNode(Compiler compiler) + { + Argument.InnerCompile(compiler); + compiler.Instructions.Add(new(InstructionType.CALL_CHI)); + } + } } } diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs index 32b4a0cde..28381dc36 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs @@ -40,6 +40,15 @@ internal enum InstructionType CALL_FLOOR, CALL_CEIL, CALL_ROUND, + CALL_ERF, + CALL_ERFC, + CALL_ERFI, + CALL_EI, + CALL_LI, + CALL_SI, + CALL_CI, + CALL_SHI, + CALL_CHI, // 2-arg functions CALL_SUM = 100, @@ -219,6 +228,35 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value case InstructionType.CALL_ARCCOSECANT: stack.Push(Core.Compilation.ComplexBranches.Arccosecant(stack.Pop())); break; + // In double precision, as the kernels define them. + // https://github.com/asc-community/AngouriMath/issues/1607 + case InstructionType.CALL_ERF: + stack.Push(Core.Compilation.SpecialFunctions.Erf(stack.Pop())); + break; + case InstructionType.CALL_ERFC: + stack.Push(Core.Compilation.SpecialFunctions.Erfc(stack.Pop())); + break; + case InstructionType.CALL_ERFI: + stack.Push(Core.Compilation.SpecialFunctions.Erfi(stack.Pop())); + break; + case InstructionType.CALL_EI: + stack.Push(Core.Compilation.SpecialFunctions.Ei(stack.Pop())); + break; + case InstructionType.CALL_LI: + stack.Push(Core.Compilation.SpecialFunctions.Li(stack.Pop())); + break; + case InstructionType.CALL_SI: + stack.Push(Core.Compilation.SpecialFunctions.Si(stack.Pop())); + break; + case InstructionType.CALL_CI: + stack.Push(Core.Compilation.SpecialFunctions.Ci(stack.Pop())); + break; + case InstructionType.CALL_SHI: + stack.Push(Core.Compilation.SpecialFunctions.Shi(stack.Pop())); + break; + case InstructionType.CALL_CHI: + stack.Push(Core.Compilation.SpecialFunctions.Chi(stack.Pop())); + break; case InstructionType.CALL_FACTORIAL: // https://stackoverflow.com/a/15454784/5429648 const int g = 7; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs index 2568988b4..6db3fa3a4 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs @@ -137,6 +137,18 @@ Arccosecantf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef Floorf => Expression.Call(GetDef("Floor", 1, e.Type), e), Ceilf => Expression.Call(GetDef("Ceil", 1, e.Type), e), Roundf => Expression.Call(GetDef("Round", 1, e.Type), e), + + // The special functions, widened to double first as the trigonometry is. + // https://github.com/asc-community/AngouriMath/issues/1607 + Erff when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Erf", 1, newE.Type), newE), + Erfcf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Erfc", 1, newE.Type), newE), + Erfif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Erfi", 1, newE.Type), newE), + Eif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Ei", 1, newE.Type), newE), + Lif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Li", 1, newE.Type), newE), + Sif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Si", 1, newE.Type), newE), + 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), Notf => Expression.Not(e), diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs index 0b687f76d..e5d2b0594 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs @@ -288,6 +288,27 @@ 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 + public static System.Numerics.Complex Erf(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Erf(a); + public static double Erf(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Erf(a); + public static System.Numerics.Complex Erfc(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfc(a); + public static double Erfc(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfc(a); + public static System.Numerics.Complex Erfi(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfi(a); + public static double Erfi(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfi(a); + public static System.Numerics.Complex Ei(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Ei(a); + public static double Ei(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Ei(a); + public static System.Numerics.Complex Li(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Li(a); + public static double Li(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Li(a); + public static System.Numerics.Complex Si(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Si(a); + public static double Si(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Si(a); + public static System.Numerics.Complex Ci(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Ci(a); + public static double Ci(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Ci(a); + public static System.Numerics.Complex Shi(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Shi(a); + public static double Shi(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Shi(a); + public static System.Numerics.Complex Chi(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Chi(a); + public static double Chi(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Chi(a); + /// /// 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 @@ -438,6 +459,24 @@ public static bool IsNaN(System.Numerics.Complex a) { ("Round", 1, typeof(long)), Def(a => Round(a)) }, { ("Round", 1, typeof(int)), Def(a => Round(a)) }, { ("Round", 1, typeof(BigInteger)), Def(a => Round(a)) }, + { ("Erf", 1, typeof(System.Numerics.Complex)), Def(a => Erf(a)) }, + { ("Erf", 1, typeof(double)), Def(a => Erf(a)) }, + { ("Erfc", 1, typeof(System.Numerics.Complex)), Def(a => Erfc(a)) }, + { ("Erfc", 1, typeof(double)), Def(a => Erfc(a)) }, + { ("Erfi", 1, typeof(System.Numerics.Complex)), Def(a => Erfi(a)) }, + { ("Erfi", 1, typeof(double)), Def(a => Erfi(a)) }, + { ("Ei", 1, typeof(System.Numerics.Complex)), Def(a => Ei(a)) }, + { ("Ei", 1, typeof(double)), Def(a => Ei(a)) }, + { ("Li", 1, typeof(System.Numerics.Complex)), Def(a => Li(a)) }, + { ("Li", 1, typeof(double)), Def(a => Li(a)) }, + { ("Si", 1, typeof(System.Numerics.Complex)), Def(a => Si(a)) }, + { ("Si", 1, typeof(double)), Def(a => Si(a)) }, + { ("Ci", 1, typeof(System.Numerics.Complex)), Def(a => Ci(a)) }, + { ("Ci", 1, typeof(double)), Def(a => Ci(a)) }, + { ("Shi", 1, typeof(System.Numerics.Complex)), Def(a => Shi(a)) }, + { ("Shi", 1, typeof(double)), Def(a => Shi(a)) }, + { ("Chi", 1, typeof(System.Numerics.Complex)), Def(a => Chi(a)) }, + { ("Chi", 1, typeof(double)), Def(a => Chi(a)) }, }; // 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 f89bb3b67..f539284af 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt @@ -9,6 +9,9 @@ "Cot", "Sec", "Csc", "Acot", "Asec", "Acsc", "Abs", "Sgn", "Floor", "Ceil", "Round" }; 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 types = new[]{ "System.Numerics.Complex", "double", "float", "long", "int", "BigInteger" }; #> // @@ -188,6 +191,13 @@ 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 +<# foreach (var name in special) { #> + public static System.Numerics.Complex <#= name #>(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.<#= name #>(a); + public static double <#= name #>(double a) => AngouriMath.Core.Compilation.SpecialFunctions.<#= name #>(a); +<# } #> + /// /// 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 @@ -221,6 +231,11 @@ namespace AngouriMath.Core.Compilation.IntoLinq <# foreach (var type in types) { #> { ("<#= name #>", 1, typeof(<#= type #>)), Def<<#= type #>>(a => <#= name #>(a)) }, <# } #> +<# } #> +<# foreach (var name in special) { #> +<# foreach (var type in new[]{ "System.Numerics.Complex", "double" }) { #> + { ("<#= name #>", 1, typeof(<#= type #>)), Def<<#= type #>>(a => <#= name #>(a)) }, +<# } #> <# } #> }; diff --git a/Sources/AngouriMath/Functions/Compilation/SpecialFunctions.cs b/Sources/AngouriMath/Functions/Compilation/SpecialFunctions.cs new file mode 100644 index 000000000..26bc9023c --- /dev/null +++ b/Sources/AngouriMath/Functions/Compilation/SpecialFunctions.cs @@ -0,0 +1,405 @@ +// +// 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 NumericsComplex = System.Numerics.Complex; + +namespace AngouriMath.Core.Compilation +{ + /// + /// The special functions in double precision, for both compilers: erf, erfc, + /// erfi, Ei, li, Si, Ci, Shi and Chi. + /// The interpreter's kernels, and the + /// others, work in arbitrary precision. At 20 digits they take 0.7 to 5.5 ms a call. Newton's + /// method makes thousands of calls. These take a microsecond or two. + /// https://github.com/asc-community/AngouriMath/issues/1607 + /// + /// + /// + /// Each is the kernel's function, with the kernel's branch cuts and the kernel's values on the + /// axes. They were measured against the kernels at 30 digits on some ten thousand points: a + /// grid of [-14, 14]^2, both axes out to 40 and down to 1e-4, either side of the + /// negative real axis, and circles out to 300. The integrals agree to within 6e-14 of the + /// value. The error functions agree to within 3e-13, which is the rounding of + /// z^2 in their e^(-z^2) where |z| is large. An imaginary part of + /// either sign of zero is read as the axis, as the kernels read it. + /// + /// + /// The error functions come from the Faddeeva function + /// w(z) = e^(-z^2) erfc(-i z), by Weideman's rational approximation with 40 terms + /// (J. A. C. Weideman, Computation of the complex error function, SIAM J. Numer. Anal. + /// 31, 1994), and near the origin from the Taylor series of erf. + /// + /// + /// Ei is summed as its series where the terms cancel little, which is where + /// |z| - Re z is at most 3 or |z| at most 2. Out to |z| = 40 it is + /// -E1(-z) otherwise, by E1's continued fraction, plus the branch term. From + /// 40 on it is the asymptotic series. li is Ei of the logarithm. Shi and + /// Chi are summed as their own series by the same rule with |Re z|, and + /// otherwise come from Ei(z) and Ei(-z). Si and Ci are those at + /// i z. Every relation is one the kernels state. + /// + /// + internal static class SpecialFunctions + { + private const double SqrtPi = 1.7724538509055160273; + private const double EulerGamma = 0.57721566490153286061; + + /// The real part where the value is real, and NaN where it is not, as of a negative number is. + private static double RealOrNaN(NumericsComplex value) => value.Imaginary == 0 ? value.Real : double.NaN; + + private static bool IsNaN(NumericsComplex z) => double.IsNaN(z.Real) || double.IsNaN(z.Imaginary); + private static bool IsInfinite(NumericsComplex z) => double.IsInfinity(z.Real) || double.IsInfinity(z.Imaginary); + [ConstantField] private static readonly NumericsComplex NaN = new(double.NaN, double.NaN); + + public static NumericsComplex Erf(NumericsComplex z) + { + if (IsNaN(z)) return NaN; + if (IsInfinite(z)) + return z.Imaginary == 0 ? (z.Real < 0 ? -1 : 1) : NaN; + NumericsComplex value; + if (z.Real * z.Real + z.Imaginary * z.Imaginary < 1) + value = ErfByTheSeries(z); + else if (z.Real >= 0) + value = 1 - ExpTimes(-z * z, Faddeeva(new NumericsComplex(-z.Imaginary, z.Real))); + else + value = ExpTimes(-z * z, Faddeeva(new NumericsComplex(z.Imaginary, -z.Real))) - 1; + // Real on the real line and imaginary on the imaginary axis. + if (z.Imaginary == 0) return new(value.Real, 0); + if (z.Real == 0) return new(0, value.Imaginary); + return value; + } + + public static NumericsComplex Erfc(NumericsComplex z) + { + if (IsNaN(z)) return NaN; + if (IsInfinite(z)) + return z.Imaginary == 0 ? (z.Real < 0 ? 2 : 0) : NaN; + NumericsComplex value; + if (z.Real * z.Real + z.Imaginary * z.Imaginary < 1) + value = 1 - ErfByTheSeries(z); + else if (z.Real >= 0) + value = ExpTimes(-z * z, Faddeeva(new NumericsComplex(-z.Imaginary, z.Real))); + else + value = 2 - ExpTimes(-z * z, Faddeeva(new NumericsComplex(z.Imaginary, -z.Real))); + // Real on the real line, and 1 - i erfi(y) on the imaginary axis. + if (z.Imaginary == 0) return new(value.Real, 0); + if (z.Real == 0) return new(1, value.Imaginary); + return value; + } + + /// -i erf(i z). + public static NumericsComplex Erfi(NumericsComplex z) + { + var erf = Erf(new NumericsComplex(-z.Imaginary, z.Real)); + // -i (a + i b) = b - i a + var value = new NumericsComplex(erf.Imaginary, -erf.Real); + if (z.Imaginary == 0) return new(value.Real, 0); + if (z.Real == 0) return new(0, value.Imaginary); + return value; + } + + public static double Erf(double x) => RealOrNaN(Erf(new NumericsComplex(x, 0))); + public static double Erfc(double x) => RealOrNaN(Erfc(new NumericsComplex(x, 0))); + public static double Erfi(double x) => RealOrNaN(Erfi(new NumericsComplex(x, 0))); + + /// 2/sqrt(pi) sum (-1)^n z^(2n + 1)/(n! (2n + 1)), for |z| < 1. + private static NumericsComplex ErfByTheSeries(NumericsComplex z) + { + var minusZSquared = -z * z; + NumericsComplex term = z, sum = z; + for (var n = 1; n < 100; n++) + { + term *= minusZSquared / n; + var next = term / (2 * n + 1); + sum += next; + if (next.Magnitude <= 1e-17 * sum.Magnitude) + break; + } + return 2 / SqrtPi * sum; + } + + /// + /// w(z) for Im z >= 0, Weideman's + /// 2 p(Z)/(L - i z)^2 + 1/(sqrt(pi) (L - i z)) with Z = (L + i z)/(L - i z). + /// + private static NumericsComplex Faddeeva(NumericsComplex z) + { + // L - i z and L + i z, for z = x + i y + var below = new NumericsComplex(WeidemanL + z.Imaginary, -z.Real); + var above = new NumericsComplex(WeidemanL - z.Imaginary, z.Real); + var zz = above / below; + NumericsComplex p = weideman[WeidemanTerms]; + for (var n = WeidemanTerms - 1; n >= 1; n--) + p = p * zz + weideman[n]; + return 2 * p / (below * below) + 1 / (SqrtPi * below); + } + + private const int WeidemanTerms = 40; + [ConstantField] private static readonly double WeidemanL = Math.Sqrt(WeidemanTerms / Math.Sqrt(2)); + + /// + /// The polynomial's coefficients, a_n = 1/(4N) sum_k f(t_k) cos(pi n k/(2N)) over + /// |k| < 2N, with t_k = L tan(pi k/(4N)) and f(t) = e^(-t^2) (L^2 + t^2): + /// the discrete Fourier transform Weideman takes, written out, since f is even. + /// Index 0 is unused. + /// + [ConstantField] private static readonly double[] weideman = WeidemanCoefficients(); + + private static double[] WeidemanCoefficients() + { + var m = 2 * WeidemanTerms; + var coefficients = new double[WeidemanTerms + 1]; + for (var n = 1; n <= WeidemanTerms; n++) + { + var sum = 0.0; + for (var k = -m + 1; k <= m - 1; k++) + { + var t = WeidemanL * Math.Tan(k * Math.PI / (2 * m)); + sum += Math.Exp(-t * t) * (WeidemanL * WeidemanL + t * t) * Math.Cos(Math.PI * n * k / m); + } + coefficients[n] = sum / (4 * WeidemanTerms); + } + return coefficients; + } + + /// + /// e^a f. Where e^a overflows, each part of the product is scaled in logarithms, + /// so that an infinite factor and a finite one make an infinity and not a NaN, and a part with + /// a small cosine or sine stays finite if it is. + /// + private static NumericsComplex ExpTimes(NumericsComplex a, NumericsComplex f) + { + if (a.Real <= 700) + return NumericsComplex.Exp(a) * f; + var s = a + NumericsComplex.Log(f); + return new(Scaled(s.Real, Math.Cos(s.Imaginary)), Scaled(s.Real, Math.Sin(s.Imaginary))); + } + + private static double Scaled(double exponent, double factor) + => factor == 0 ? 0 : Math.Sign(factor) * Math.Exp(exponent + Math.Log(Math.Abs(factor))); + + /// The principal logarithm, with a zero imaginary part of either sign read as the axis. + private static NumericsComplex Ln(NumericsComplex z) + => new(Math.Log(z.Magnitude), z.Imaginary == 0 ? (z.Real < 0 ? Math.PI : 0) : Math.Atan2(z.Imaginary, z.Real)); + + /// i pi sgn(Im z), the term the logarithm's branch adds off the real axis. + private static NumericsComplex Branch(NumericsComplex z) + => z.Imaginary == 0 ? NumericsComplex.Zero : new(0, Math.Sign(z.Imaginary) * Math.PI); + + /// + /// Ei(z) = gamma + (ln z - ln(1/z))/2 + sum_(k >= 1) z^k/(k k!), real on the whole real + /// line, as defines it. + /// + public static NumericsComplex Ei(NumericsComplex z) + { + if (IsNaN(z)) return NaN; + if (IsInfinite(z)) + return z.Imaginary == 0 ? (z.Real < 0 ? 0 : double.PositiveInfinity) : NaN; + if (z == NumericsComplex.Zero) + return double.NegativeInfinity; + var modulus = z.Magnitude; + NumericsComplex value; + if (modulus >= 40) + value = EiByTheAsymptoticSeries(z) + Branch(z); + else if (modulus <= 2 || modulus - z.Real <= 3) + // (ln z - ln(1/z))/2: ln |z| on the negative axis. + value = EulerGamma + (z.Imaginary == 0 && z.Real < 0 ? Math.Log(-z.Real) : Ln(z)) + ExponentialSeries(z); + else + value = Branch(z) - E1ByTheContinuedFraction(-z); + return z.Imaginary == 0 ? new(value.Real, 0) : value; + } + + /// li(z) = Ei(ln z), with the principal logarithm; li(0) = 0 and li(1) = -oo. + public static NumericsComplex Li(NumericsComplex z) + { + if (IsNaN(z)) return NaN; + if (IsInfinite(z)) + return z.Imaginary == 0 && z.Real > 0 ? double.PositiveInfinity : NaN; + if (z == NumericsComplex.Zero) + return NumericsComplex.Zero; + if (z == NumericsComplex.One) + return double.NegativeInfinity; + var value = Ei(Ln(z)); + return z.Imaginary == 0 && z.Real > 0 ? new(value.Real, 0) : value; + } + + public static double Ei(double x) => RealOrNaN(Ei(new NumericsComplex(x, 0))); + public static double Li(double x) => RealOrNaN(Li(new NumericsComplex(x, 0))); + + /// sum_(k >= 1) z^k/(k k!). + private static NumericsComplex ExponentialSeries(NumericsComplex z) + { + NumericsComplex power = 1, sum = 0; + var modulus = z.Magnitude; + for (var k = 1; k < 1000; k++) + { + power = power * z / k; + var term = power / k; + sum += term; + if (k > modulus && term.Magnitude <= 1e-17 * (sum.Magnitude + 1)) + break; + } + return sum; + } + + /// + /// E1(w) = e^(-w)/(w + 1 - 1/(w + 3 - 4/(w + 5 - ...))), by the modified Lentz method. + /// Used where w is at least 2 out and away from its cut, the negative real axis, and + /// there it takes at most a couple of hundred steps. + /// + private static NumericsComplex E1ByTheContinuedFraction(NumericsComplex w) + { + const double tiny = 1e-300; + var b = w + 1; + NumericsComplex c = 1 / tiny, d = 1 / b, h = d; + for (var i = 1; i < 5000; i++) + { + var a = -(double)i * i; + b += 2; + d = a * d + b; + if (d == NumericsComplex.Zero) d = tiny; + c = b + a / c; + if (c == NumericsComplex.Zero) c = tiny; + d = 1 / d; + var delta = c * d; + h *= delta; + if ((delta - 1).Magnitude < 1e-16) + break; + } + return ExpTimes(-w, h); + } + + /// + /// e^z/z sum_k k!/z^k, cut at its smallest term, which from |z| = 40 on is below + /// 1e-16 of the sum. + /// + private static NumericsComplex EiByTheAsymptoticSeries(NumericsComplex z) + { + NumericsComplex term = 1, sum = 1; + var previous = 1.0; + for (var k = 1; k < 200; k++) + { + var next = term * k / z; + var size = next.Magnitude; + if (size >= previous || size < 1e-17) + break; + sum += next; + (term, previous) = (next, size); + } + return ExpTimes(z, sum / z); + } + + public static NumericsComplex Si(NumericsComplex z) + { + if (IsNaN(z)) return NaN; + if (IsInfinite(z)) + return z.Imaginary == 0 ? (z.Real < 0 ? -Math.PI / 2 : Math.PI / 2) : NaN; + // -i Shi(i z) + var shi = HyperbolicIntegral(new NumericsComplex(-z.Imaginary, z.Real), odd: true); + return OnTheAxes(z, new NumericsComplex(shi.Imaginary, -shi.Real), odd: true); + } + + public static NumericsComplex Ci(NumericsComplex z) + { + if (IsNaN(z)) return NaN; + if (z == NumericsComplex.Zero) return double.NegativeInfinity; + if (IsInfinite(z)) + return z.Imaginary == 0 ? (z.Real > 0 ? NumericsComplex.Zero : new NumericsComplex(0, Math.PI)) : NaN; + // Chi(i z) + ln z - ln(i z), the last two -i pi/2 where arg z <= pi/2 and 3 i pi/2 past it. + var chi = HyperbolicIntegral(new NumericsComplex(-z.Imaginary, z.Real), odd: false); + var argument = z.Imaginary == 0 ? (z.Real < 0 ? Math.PI : 0) : Math.Atan2(z.Imaginary, z.Real); + return OnTheAxes(z, chi + new NumericsComplex(0, argument <= Math.PI / 2 ? -Math.PI / 2 : 3 * Math.PI / 2), odd: false); + } + + public static NumericsComplex Shi(NumericsComplex z) => HyperbolicIntegral(z, odd: true); + public static NumericsComplex Chi(NumericsComplex z) => HyperbolicIntegral(z, odd: false); + + public static double Si(double x) => RealOrNaN(Si(new NumericsComplex(x, 0))); + public static double Ci(double x) => RealOrNaN(Ci(new NumericsComplex(x, 0))); + public static double Shi(double x) => RealOrNaN(Shi(new NumericsComplex(x, 0))); + public static double Chi(double x) => RealOrNaN(Chi(new NumericsComplex(x, 0))); + + /// Shi where , and Chi otherwise. + private static NumericsComplex HyperbolicIntegral(NumericsComplex z, bool odd) + { + if (IsNaN(z)) return NaN; + if (z == NumericsComplex.Zero) + return odd ? NumericsComplex.Zero : double.NegativeInfinity; + if (IsInfinite(z)) + { + if (z.Imaginary != 0) return NaN; + return odd ? (z.Real < 0 ? double.NegativeInfinity : double.PositiveInfinity) + : (z.Real < 0 ? NaN : double.PositiveInfinity); + } + var modulus = z.Magnitude; + NumericsComplex value; + if (modulus < 40 && (modulus <= 2 || modulus - Math.Abs(z.Real) <= 3)) + value = odd ? HalfOfTheExponentialSeries(z, odd: true) : EulerGamma + Ln(z) + HalfOfTheExponentialSeries(z, odd: false); + else + { + // Shi(z) = (Ei(z) - Ei(-z))/2 - i pi/2 sgn(Im z), Chi(z) = (Ei(z) + Ei(-z))/2 + i pi/2 + // sgn(Im z), and i pi more for Chi on the negative real axis. + var plus = Ei(z); + var minus = Ei(-z); + value = Half(odd ? plus - minus : plus + minus); + if (z.Imaginary != 0) + { + var branch = new NumericsComplex(0, Math.Sign(z.Imaginary) * Math.PI / 2); + value = odd ? value - branch : value + branch; + } + else if (!odd && z.Real < 0) + value += new NumericsComplex(0, Math.PI); + } + return OnTheAxes(z, value, odd); + } + + /// + /// sum z^(2k + 1)/((2k + 1) (2k + 1)!) where , and + /// sum_(k >= 1) z^(2k)/(2k (2k)!) otherwise. + /// + private static NumericsComplex HalfOfTheExponentialSeries(NumericsComplex z, bool odd) + { + NumericsComplex power = 1, sum = 0; + var modulus = z.Magnitude; + for (var n = 1; n < 2000; n++) + { + power = power * z / n; + if ((n % 2 == 1) != odd) + continue; + var term = power / n; + sum += term; + if (n > modulus && term.Magnitude <= 1e-17 * (sum.Magnitude + 1)) + break; + } + return sum; + } + + /// + /// Half of , 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. + /// + private static NumericsComplex Half(NumericsComplex value) => new(value.Real / 2, value.Imaginary / 2); + + /// + /// What the kernels give exactly on the axes: the odd ones real on the real line and imaginary on + /// the imaginary axis; the even ones real on the positive real axis, i pi more on the + /// negative, and i pi/2 sgn(y) more than a real value on the imaginary axis. That holds + /// however large the value, so an infinite one is given its exact part as well. + /// + private static NumericsComplex OnTheAxes(NumericsComplex z, NumericsComplex value, bool odd) + { + if (double.IsNaN(value.Real) && double.IsNaN(value.Imaginary)) + return value; + if (z.Imaginary == 0) + return new(value.Real, odd || z.Real >= 0 ? 0 : Math.PI); + if (z.Real == 0) + return odd ? new NumericsComplex(0, value.Imaginary) : new NumericsComplex(value.Real, Math.Sign(z.Imaginary) * Math.PI / 2); + return value; + } + } +} diff --git a/Sources/Tests/UnitTests/Algebra/SolveTest/NewtonFallbackTest.cs b/Sources/Tests/UnitTests/Algebra/SolveTest/NewtonFallbackTest.cs index d2631fac7..b85dce666 100644 --- a/Sources/Tests/UnitTests/Algebra/SolveTest/NewtonFallbackTest.cs +++ b/Sources/Tests/UnitTests/Algebra/SolveTest/NewtonFallbackTest.cs @@ -18,9 +18,10 @@ namespace AngouriMath.Tests.Algebra /// or the unevaluated derivative(max(x, 1), x) that max differentiated to -- the /// solver threw out of the public method, in 2.5.0 as on /// master. The floors, the rounding and the extremes compile now, and max and min - /// differentiate; what still has no compiled form is declined, and so is a set of roots that - /// are not isolated. + /// differentiate, and so do the special functions; what still has no compiled form is declined, + /// and so is a set of roots that are not isolated. /// #1603 + /// #1607 /// [Trait("Area", "Algebra")] public sealed class NewtonFallbackTest @@ -44,13 +45,36 @@ public void TheFloorsAndTheExtremesAreSearched(string equation, double[] expecte Assert.Contains(expected, e => System.Math.Abs(e - root.EvalNumerical().RealPart.EDecimal.ToDouble()) < 1e-9); } - /// No analytical route answers these, and the compiler still has no form for them. + /// + /// No analytical route answers these, and the compiler has no form for the binomial + /// coefficient or for gcd, nor for their unevaluated derivatives. + /// [Theory] - [InlineData("erf(x) = x/2 + 1/7")] - [InlineData("Si(x) = x/2 + 1/7")] + [InlineData("binomial(x, 3/2) = x/2")] + [InlineData("gcd(x, 6) = x/2")] public void AnEquationTheCompilerHasNoFormForIsLeftUnsolved(string equation) => Assert.IsType(equation.ToEntity().Solve("x")); + /// + /// The special functions compile, so the search reaches these. Among the roots it finds off + /// the real line are the three real ones, each checked against mpmath. Every root it gives + /// holds, to the precision of the double it was found in. + /// #1607 + /// + [Theory] + [InlineData("erf(x) - x/2 - 1/7", new[] { -2.2832295840505545496, 0.23498157773996424592, 1.6791581585243440031 })] + [InlineData("Si(x) - x/2 - 1/7", new[] { -3.8544741217072080362, 0.28837215854663464706, 3.3983595626337253252 })] + public void TheSpecialFunctionsAreSearched(string expression, double[] realRoots) + { + var roots = Assert.IsType((expression + " = 0").ToEntity().Solve("x")); + foreach (var expected in realRoots) + Assert.Contains(roots, root => root.EvalNumerical() is var value + && System.Math.Abs(value.ImaginaryPart.EDecimal.ToDouble()) < 1e-12 + && System.Math.Abs(value.RealPart.EDecimal.ToDouble() - expected) < 1e-9); + foreach (var root in roots) + Assert.True(expression.ToEntity().Substitute("x", root).EvalNumerical().Abs() < 1e-9, $"{root} is not a root"); + } + /// /// Each holds on the whole of [0, +oo), and the search's grid converged to fifty-one points /// of [0, 10] as though they were the set. @@ -74,7 +98,7 @@ public void AnEquationTheCompilerHasAFormForIsStillSearched() /// [Fact] public void NewtonsMethodAskedByNameStillSaysWhy() - => Assert.Throws(() => "erf(x) - x/2 - 1/7".ToEntity().SolveNt("x")); + => Assert.Throws(() => "binomial(x, 3/2) - x/2".ToEntity().SolveNt("x")); /// /// The compiled floors, rounding and extremes against the interpreter, off the real line diff --git a/Sources/Tests/UnitTests/Convenience/SpecialFunctionCompilationTest.cs b/Sources/Tests/UnitTests/Convenience/SpecialFunctionCompilationTest.cs new file mode 100644 index 000000000..0732104c1 --- /dev/null +++ b/Sources/Tests/UnitTests/Convenience/SpecialFunctionCompilationTest.cs @@ -0,0 +1,118 @@ +// +// 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 AngouriMath.Extensions; +using Xunit; +using NumericsComplex = System.Numerics.Complex; + +namespace AngouriMath.Tests.Convenience +{ + /// + /// The special functions compile, through Compile and Compile<TIn, TOut> + /// alike, in double precision and in agreement with the interpreter: its branch cuts, and its + /// values on the axes. + /// #1607 + /// + [Trait("Area", "Convenience")] + public sealed class SpecialFunctionCompilationTest + { + public static readonly TheoryData Points = new() + { + // Near the origin, where erf is its Taylor series, and the others their own series. + { "erf(x)", 0.3, 0.2 }, { "erf(x)", 2.5, 0 }, { "erf(x)", -1.7, 0 }, { "erf(x)", 1.5, -2.5 }, + { "erf(x)", -3, 4 }, { "erf(x)", 0, 3 }, { "erf(x)", 20, 30 }, + { "erfc(x)", 0.3, 0.2 }, { "erfc(x)", 2.5, 0 }, { "erfc(x)", -1.7, 0 }, { "erfc(x)", 1.5, -2.5 }, + { "erfc(x)", 5, 1 }, { "erfc(x)", 0, 3 }, + { "erfi(x)", 0.3, 0.2 }, { "erfi(x)", 2.5, 0 }, { "erfi(x)", -1.7, 0 }, { "erfi(x)", 1.5, -2.5 }, + // Ei, by its series, by E1's continued fraction, by the asymptotic series, and either side + // of its cut, where it is i pi more or less than on the axis. + { "Ei(x)", 0.3, 0.2 }, { "Ei(x)", 2.5, 0 }, { "Ei(x)", -1.7, 0 }, { "Ei(x)", -3, 4 }, + { "Ei(x)", 12, 0 }, { "Ei(x)", -5, 1e-9 }, { "Ei(x)", -5, -1e-9 }, { "Ei(x)", 45, 10 }, + { "Ei(x)", -50, 0 }, { "Ei(x)", 0, 3 }, + { "li(x)", 2.5, 0 }, { "li(x)", 0.4, 0 }, { "li(x)", -2, 0 }, { "li(x)", 1.5, -2.5 }, { "li(x)", 1e6, 0 }, + { "Si(x)", 0.3, 0.2 }, { "Si(x)", 2.5, 0 }, { "Si(x)", -1.7, 0 }, { "Si(x)", 1.5, -2.5 }, + { "Si(x)", 20, 3 }, { "Si(x)", 0, 3 }, { "Si(x)", 50, 0 }, + { "Ci(x)", 0.3, 0.2 }, { "Ci(x)", 2.5, 0 }, { "Ci(x)", -1.7, 0 }, { "Ci(x)", -5, 1e-9 }, + { "Ci(x)", -5, -1e-9 }, { "Ci(x)", 0, 3 }, { "Ci(x)", 0, -3 }, { "Ci(x)", 20, 3 }, + { "Shi(x)", 0.3, 0.2 }, { "Shi(x)", 2.5, 0 }, { "Shi(x)", -1.7, 0 }, { "Shi(x)", 3, 20 }, + { "Chi(x)", 0.3, 0.2 }, { "Chi(x)", 2.5, 0 }, { "Chi(x)", -1.7, 0 }, { "Chi(x)", -5, 1e-9 }, + { "Chi(x)", 3, 20 }, { "Chi(x)", 0, 3 }, + }; + + private static NumericsComplex Interpreted(string expression, double re, double im) + { + // With the downcasting off, so that a small value is not read as zero. + using var _ = MathS.Settings.DowncastingEnabled.Set(false); + return ((Entity.Number.Complex)expression.ToEntity().Substitute("x", Entity.Number.Complex.Create(re, im)).EvalNumerical()).ToNumerics(); + } + + private static void Close(NumericsComplex expected, NumericsComplex actual) + => Assert.True((actual - expected).Magnitude <= 1e-12 * expected.Magnitude, $"expected {expected}, got {actual}"); + + [Theory] + [MemberData(nameof(Points))] + public void TheCompiledFormAgreesWithTheInterpreter(string expression, double re, double im) + => Close(Interpreted(expression, re, im), expression.ToEntity().Compile("x").Call(new NumericsComplex(re, im))); + + [Theory] + [MemberData(nameof(Points))] + public void TheExpressionTreeAgreesWithTheInterpreter(string expression, double re, double im) + => Close(Interpreted(expression, re, im), + expression.ToEntity().Compile("x")(new NumericsComplex(re, im))); + + /// + /// Over doubles, a value where the function is real, and NaN where it is not, as + /// of a negative number is: li, Ci and Chi + /// are i pi more than a real value left of 0. + /// + [Theory] + [InlineData("erf(x)", 2.5)] + [InlineData("erfc(x)", -1.7)] + [InlineData("erfi(x)", 2.5)] + [InlineData("Ei(x)", -1.7)] + [InlineData("li(x)", 0.4)] + [InlineData("Si(x)", -1.7)] + [InlineData("Ci(x)", 2.5)] + [InlineData("Shi(x)", -1.7)] + [InlineData("Chi(x)", 2.5)] + public void OverDoublesARealValueCompiles(string expression, double at) + { + var expected = Interpreted(expression, at, 0); + Assert.Equal(0, expected.Imaginary); + var compiled = expression.ToEntity().Compile("x")(at); + Assert.True(System.Math.Abs(compiled - expected.Real) <= 1e-12 * System.Math.Abs(expected.Real), $"expected {expected.Real}, got {compiled}"); + } + + /// + /// A value too large for a double is infinite in its large part and keeps the part that is + /// exact on the axis, as the interpreter's does. Halving Ei(z) - Ei(-z) as a complex + /// number made that part NaN, since an infinity times the other part's zero is one. + /// + [Theory] + [InlineData("Shi(x)", 1000, 0, double.PositiveInfinity, 0)] + [InlineData("Shi(x)", -1000, 0, double.NegativeInfinity, 0)] + [InlineData("Chi(x)", -1000, 0, double.PositiveInfinity, System.Math.PI)] + [InlineData("Si(x)", 0, 1000, 0, double.PositiveInfinity)] + [InlineData("Ci(x)", 0, 1000, double.PositiveInfinity, System.Math.PI / 2)] + public void AnInfiniteValueKeepsItsExactPart(string expression, double re, double im, double expectedRe, double expectedIm) + { + var compiled = expression.ToEntity().Compile("x").Call(new NumericsComplex(re, im)); + Assert.Equal(expectedRe, compiled.Real); + Assert.Equal(expectedIm, compiled.Imaginary); + var interpreted = Interpreted(expression, re, im); + Assert.Equal(expectedRe, interpreted.Real); + Assert.Equal(expectedIm, interpreted.Imaginary, 15); + } + + [Theory] + [InlineData("li(x)", -2)] + [InlineData("Ci(x)", -1.7)] + [InlineData("Chi(x)", -1.7)] + public void OverDoublesAValueOffTheRealLineIsNaN(string expression, double at) + => Assert.True(double.IsNaN(expression.ToEntity().Compile("x")(at))); + } +} From 19678c932cc7b94cd8b2a0f49bc4dda6c3a1242a Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 30 Sep 2026 08:30:31 +0000 Subject: [PATCH 2/2] SpecialFunctions moves to Numerics, beside the interval evaluation The routines are the compilers' today, and the numerical evaluation's fast tier is where they are to be shared next, as asked on #1607: the interval evaluation reads no special function, so an integrator check with erf or Ei in it falls back to the decimal evaluation. Nothing but the namespace and the file's place changes. Part of #1607. Co-Authored-By: Claude Opus 5.5 --- Sources/.editorconfig | 2 +- .../Compilation/IntoFE/FastExpression.cs | 18 +++++----- .../Compilation/IntoLinq/MathAllMethods.cs | 36 +++++++++---------- .../Compilation/IntoLinq/MathAllMethods.tt | 4 +-- .../SpecialFunctions.cs | 13 +++---- 5 files changed, 37 insertions(+), 36 deletions(-) rename Sources/AngouriMath/{Functions/Compilation => Numerics}/SpecialFunctions.cs (96%) diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 7667f5a13..9a89de0e1 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -342,7 +342,7 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [AngouriMath/Core/Entity/Continuous/Number/ErrorFunctions.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/Functions/Compilation/SpecialFunctions.cs] +[AngouriMath/Numerics/SpecialFunctions.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/Core/Entity/Continuous/Entity.Continuous.ErrorFunctions.Classes.cs] diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs index 28381dc36..3d08fa302 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs @@ -231,31 +231,31 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value // In double precision, as the kernels define them. // https://github.com/asc-community/AngouriMath/issues/1607 case InstructionType.CALL_ERF: - stack.Push(Core.Compilation.SpecialFunctions.Erf(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Erf(stack.Pop())); break; case InstructionType.CALL_ERFC: - stack.Push(Core.Compilation.SpecialFunctions.Erfc(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Erfc(stack.Pop())); break; case InstructionType.CALL_ERFI: - stack.Push(Core.Compilation.SpecialFunctions.Erfi(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Erfi(stack.Pop())); break; case InstructionType.CALL_EI: - stack.Push(Core.Compilation.SpecialFunctions.Ei(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Ei(stack.Pop())); break; case InstructionType.CALL_LI: - stack.Push(Core.Compilation.SpecialFunctions.Li(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Li(stack.Pop())); break; case InstructionType.CALL_SI: - stack.Push(Core.Compilation.SpecialFunctions.Si(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Si(stack.Pop())); break; case InstructionType.CALL_CI: - stack.Push(Core.Compilation.SpecialFunctions.Ci(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Ci(stack.Pop())); break; case InstructionType.CALL_SHI: - stack.Push(Core.Compilation.SpecialFunctions.Shi(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Shi(stack.Pop())); break; case InstructionType.CALL_CHI: - stack.Push(Core.Compilation.SpecialFunctions.Chi(stack.Pop())); + stack.Push(AngouriMath.Numerics.SpecialFunctions.Chi(stack.Pop())); break; case InstructionType.CALL_FACTORIAL: // https://stackoverflow.com/a/15454784/5429648 diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs index e5d2b0594..46b76d12a 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs @@ -290,24 +290,24 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a) // The special functions, in double precision as the kernels define them. // https://github.com/asc-community/AngouriMath/issues/1607 - public static System.Numerics.Complex Erf(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Erf(a); - public static double Erf(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Erf(a); - public static System.Numerics.Complex Erfc(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfc(a); - public static double Erfc(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfc(a); - public static System.Numerics.Complex Erfi(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfi(a); - public static double Erfi(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Erfi(a); - public static System.Numerics.Complex Ei(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Ei(a); - public static double Ei(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Ei(a); - public static System.Numerics.Complex Li(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Li(a); - public static double Li(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Li(a); - public static System.Numerics.Complex Si(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Si(a); - public static double Si(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Si(a); - public static System.Numerics.Complex Ci(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Ci(a); - public static double Ci(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Ci(a); - public static System.Numerics.Complex Shi(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Shi(a); - public static double Shi(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Shi(a); - public static System.Numerics.Complex Chi(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.Chi(a); - public static double Chi(double a) => AngouriMath.Core.Compilation.SpecialFunctions.Chi(a); + 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); + public static double Erfc(double a) => AngouriMath.Numerics.SpecialFunctions.Erfc(a); + public static System.Numerics.Complex Erfi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erfi(a); + public static double Erfi(double a) => AngouriMath.Numerics.SpecialFunctions.Erfi(a); + public static System.Numerics.Complex Ei(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Ei(a); + public static double Ei(double a) => AngouriMath.Numerics.SpecialFunctions.Ei(a); + public static System.Numerics.Complex Li(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Li(a); + public static double Li(double a) => AngouriMath.Numerics.SpecialFunctions.Li(a); + public static System.Numerics.Complex Si(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Si(a); + public static double Si(double a) => AngouriMath.Numerics.SpecialFunctions.Si(a); + public static System.Numerics.Complex Ci(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Ci(a); + public static double Ci(double a) => AngouriMath.Numerics.SpecialFunctions.Ci(a); + public static System.Numerics.Complex Shi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Shi(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); /// /// True where System.Numerics.Complex.IsNaN is, spelled out because that overload is diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt index f539284af..9fda9de21 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt @@ -194,8 +194,8 @@ namespace AngouriMath.Core.Compilation.IntoLinq // The special functions, in double precision as the kernels define them. // https://github.com/asc-community/AngouriMath/issues/1607 <# foreach (var name in special) { #> - public static System.Numerics.Complex <#= name #>(System.Numerics.Complex a) => AngouriMath.Core.Compilation.SpecialFunctions.<#= name #>(a); - public static double <#= name #>(double a) => AngouriMath.Core.Compilation.SpecialFunctions.<#= name #>(a); + 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); <# } #> /// diff --git a/Sources/AngouriMath/Functions/Compilation/SpecialFunctions.cs b/Sources/AngouriMath/Numerics/SpecialFunctions.cs similarity index 96% rename from Sources/AngouriMath/Functions/Compilation/SpecialFunctions.cs rename to Sources/AngouriMath/Numerics/SpecialFunctions.cs index 26bc9023c..d19a454d1 100644 --- a/Sources/AngouriMath/Functions/Compilation/SpecialFunctions.cs +++ b/Sources/AngouriMath/Numerics/SpecialFunctions.cs @@ -8,11 +8,12 @@ using System; using NumericsComplex = System.Numerics.Complex; -namespace AngouriMath.Core.Compilation +namespace AngouriMath.Numerics { /// - /// The special functions in double precision, for both compilers: erf, erfc, - /// erfi, Ei, li, Si, Ci, Shi and Chi. + /// The special functions in double precision: erf, erfc, erfi, Ei, + /// li, Si, Ci, Shi and Chi. Both compilers call them, and they + /// sit beside the numerical evaluation so that its fast tier can take them up too. /// The interpreter's kernels, and the /// others, work in arbitrary precision. At 20 digits they take 0.7 to 5.5 ms a call. Newton's /// method makes thousands of calls. These take a microsecond or two. @@ -54,7 +55,7 @@ internal static class SpecialFunctions private static bool IsNaN(NumericsComplex z) => double.IsNaN(z.Real) || double.IsNaN(z.Imaginary); private static bool IsInfinite(NumericsComplex z) => double.IsInfinity(z.Real) || double.IsInfinity(z.Imaginary); - [ConstantField] private static readonly NumericsComplex NaN = new(double.NaN, double.NaN); + [AngouriMath.Core.ConstantField] private static readonly NumericsComplex NaN = new(double.NaN, double.NaN); public static NumericsComplex Erf(NumericsComplex z) { @@ -140,7 +141,7 @@ private static NumericsComplex Faddeeva(NumericsComplex z) } private const int WeidemanTerms = 40; - [ConstantField] private static readonly double WeidemanL = Math.Sqrt(WeidemanTerms / Math.Sqrt(2)); + [AngouriMath.Core.ConstantField] private static readonly double WeidemanL = Math.Sqrt(WeidemanTerms / Math.Sqrt(2)); /// /// The polynomial's coefficients, a_n = 1/(4N) sum_k f(t_k) cos(pi n k/(2N)) over @@ -148,7 +149,7 @@ private static NumericsComplex Faddeeva(NumericsComplex z) /// the discrete Fourier transform Weideman takes, written out, since f is even. /// Index 0 is unused. /// - [ConstantField] private static readonly double[] weideman = WeidemanCoefficients(); + [AngouriMath.Core.ConstantField] private static readonly double[] weideman = WeidemanCoefficients(); private static double[] WeidemanCoefficients() {