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..9a89de0e1 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/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] 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..3d08fa302 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(AngouriMath.Numerics.SpecialFunctions.Erf(stack.Pop())); + break; + case InstructionType.CALL_ERFC: + stack.Push(AngouriMath.Numerics.SpecialFunctions.Erfc(stack.Pop())); + break; + case InstructionType.CALL_ERFI: + stack.Push(AngouriMath.Numerics.SpecialFunctions.Erfi(stack.Pop())); + break; + case InstructionType.CALL_EI: + stack.Push(AngouriMath.Numerics.SpecialFunctions.Ei(stack.Pop())); + break; + case InstructionType.CALL_LI: + stack.Push(AngouriMath.Numerics.SpecialFunctions.Li(stack.Pop())); + break; + case InstructionType.CALL_SI: + stack.Push(AngouriMath.Numerics.SpecialFunctions.Si(stack.Pop())); + break; + case InstructionType.CALL_CI: + stack.Push(AngouriMath.Numerics.SpecialFunctions.Ci(stack.Pop())); + break; + case InstructionType.CALL_SHI: + stack.Push(AngouriMath.Numerics.SpecialFunctions.Shi(stack.Pop())); + break; + case InstructionType.CALL_CHI: + stack.Push(AngouriMath.Numerics.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..46b76d12a 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.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 /// .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..9fda9de21 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.Numerics.SpecialFunctions.<#= name #>(a); + public static double <#= name #>(double a) => AngouriMath.Numerics.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/Numerics/SpecialFunctions.cs b/Sources/AngouriMath/Numerics/SpecialFunctions.cs new file mode 100644 index 000000000..d19a454d1 --- /dev/null +++ b/Sources/AngouriMath/Numerics/SpecialFunctions.cs @@ -0,0 +1,406 @@ +// +// 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.Numerics +{ + /// + /// 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. + /// 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); + [AngouriMath.Core.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; + [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 + /// |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. + /// + [AngouriMath.Core.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))); + } +}