Skip to content

Commit fb768a9

Browse files
Rafael-SOWNetclaude
Rafael-SOWNet
and
claude
authored
The special functions compile, in double precision (#1613)
erf, erfc, erfi, Ei, li, Si, Ci, Shi and Chi had no compiled form. So Compile and Compile<TIn, TOut> 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<TIn, TOut>, 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. SpecialFunctions lives in AngouriMath.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. Part of #1607. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
1 parent fc9c0c2 commit fb768a9

10 files changed

Lines changed: 773 additions & 10 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2339,9 +2339,9 @@ resort, Newton's method, compiles the expression and its derivative. They compil
23392339
not ordered;
23402340
- `max` and `min` differentiate where their arguments are real, as `|f|` does.
23412341

2342-
The search still declines what has no compiled form, so `erf(x) = x/2 + 1/7` is left unsolved. It
2343-
now also declines a set of roots that are not isolated: where the equation holds midway between two
2344-
of them, the points are samples of a continuum and not the set
2342+
The search still declines what has no compiled form, so `binomial(x, 3/2) = x/2` is left unsolved.
2343+
It now also declines a set of roots that are not isolated: where the equation holds midway between
2344+
two of them, the points are samples of a continuum and not the set
23452345
([#1603](https://github.com/asc-community/AngouriMath/issues/1603),
23462346
[#1420](https://github.com/asc-community/AngouriMath/issues/1420)).
23472347

@@ -2410,7 +2410,10 @@ e^(-t^2) dt`, its complement `erfc(z) = 1 - erf(z)` and the imaginary error func
24102410
`erfi(z) = -i erf(i z)` are nodes (`MathS.Erf`, `MathS.Erfc`, `MathS.Erfi`), the first of #1501's
24112411
special functions. Each is differentiated, printed in LaTeX as `\operatorname{erf}` and so on, and
24122412
evaluated to the working precision anywhere in the complex plane (`Entity.Number.Erf`, `Erfc`,
2413-
`Erfi`). An equation in one of them is left unsolved rather than answered with no roots. `erf`
2413+
`Erfi`). They compile, through `Compile` and `Compile<TIn, TOut>`, in double precision
2414+
([#1607](https://github.com/asc-community/AngouriMath/issues/1607)). So an equation the solver's
2415+
numerical search reaches, `erf(x) = x/2 + 1/7`, is answered with the roots it finds, while one that
2416+
isolates `erf`, `erf(x) = 1/2`, is left unsolved rather than answered with no roots. `erf`
24142417
was refused by name, but `erfc` and `erfi` were not, so each was read as a product with an
24152418
undeclared variable.
24162419

@@ -2421,6 +2424,8 @@ undeclared variable.
24212424
| `"erfi(x)"` | `erfi * x` | `erfi(x)` |
24222425
| `"erfc(1)".Simplify()` | `erfc` | `erfc(1)`, and `0.157299207050285…` evaluated |
24232426
| `"erf(x)".Differentiate("x")` | `UnrecognizedFunctionParseException` | `2 / sqrt(pi) * e ^ (-x ^ 2)` |
2427+
| `"erf(x)".Compile("x")`, and `erfc`, `erfi` | `UnrecognizedFunctionParseException`; for `erfc` and `erfi`, `UncompilableNodeException`: not among the variables | compiled |
2428+
| `"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 |
24242429

24252430
### `Ei` and `li` are functions
24262431

@@ -2440,6 +2445,7 @@ with an undeclared variable
24402445
| `"li(2)".Simplify()` | `2 * li` | `li(2)`, and `1.045163780117492…` evaluated |
24412446
| `"Ei(1)".Simplify()` | `Ei` | `Ei(1)`, and `1.895117816355936…` evaluated |
24422447
| `"Ei(x)".Differentiate("x")` | `Ei` | `e ^ x / x` |
2448+
| `"Ei(x)".Compile("x")`, and `li` | `UncompilableNodeException`: `Ei` is not among the variables | compiled ([#1607](https://github.com/asc-community/AngouriMath/issues/1607)) |
24432449

24442450
### The exponential and logarithmic integrals are integrated to
24452451

@@ -2478,6 +2484,7 @@ was refused, so each was read as a product with an undeclared variable
24782484
| `"Ci(2)".Simplify()` | `2 * Ci` | `Ci(2)`, and `0.4229808287748649…` evaluated |
24792485
| `"Chi(1)".Simplify()` | `Chi` | `Chi(1)`, and `0.8378669409802082…` evaluated |
24802486
| `"Si(x)".Differentiate("x")` | `Si` | `sin(x) / x` |
2487+
| `"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)) |
24812488

24822489
### A quadratic below a sine or a cosine is split over its two roots
24832490

‎Sources/.editorconfig‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -342,6 +342,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed
342342
[AngouriMath/Core/Entity/Continuous/Number/ErrorFunctions.cs]
343343
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
344344

345+
[AngouriMath/Numerics/SpecialFunctions.cs]
346+
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
347+
345348
[AngouriMath/Core/Entity/Continuous/Entity.Continuous.ErrorFunctions.Classes.cs]
346349
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
347350

‎Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -442,6 +442,107 @@ private protected override void CompileNode(Compiler compiler)
442442
compiler.Instructions.Add(new(InstructionType.CALL_PHI));
443443
}
444444
}
445+
446+
// The special functions, in double precision as the kernels define them.
447+
// https://github.com/asc-community/AngouriMath/issues/1607
448+
public partial record Erff
449+
{
450+
private protected override bool CompilesItself => true;
451+
452+
private protected override void CompileNode(Compiler compiler)
453+
{
454+
Argument.InnerCompile(compiler);
455+
compiler.Instructions.Add(new(InstructionType.CALL_ERF));
456+
}
457+
}
458+
459+
public partial record Erfcf
460+
{
461+
private protected override bool CompilesItself => true;
462+
463+
private protected override void CompileNode(Compiler compiler)
464+
{
465+
Argument.InnerCompile(compiler);
466+
compiler.Instructions.Add(new(InstructionType.CALL_ERFC));
467+
}
468+
}
469+
470+
public partial record Erfif
471+
{
472+
private protected override bool CompilesItself => true;
473+
474+
private protected override void CompileNode(Compiler compiler)
475+
{
476+
Argument.InnerCompile(compiler);
477+
compiler.Instructions.Add(new(InstructionType.CALL_ERFI));
478+
}
479+
}
480+
481+
public partial record Eif
482+
{
483+
private protected override bool CompilesItself => true;
484+
485+
private protected override void CompileNode(Compiler compiler)
486+
{
487+
Argument.InnerCompile(compiler);
488+
compiler.Instructions.Add(new(InstructionType.CALL_EI));
489+
}
490+
}
491+
492+
public partial record Lif
493+
{
494+
private protected override bool CompilesItself => true;
495+
496+
private protected override void CompileNode(Compiler compiler)
497+
{
498+
Argument.InnerCompile(compiler);
499+
compiler.Instructions.Add(new(InstructionType.CALL_LI));
500+
}
501+
}
502+
503+
public partial record Sif
504+
{
505+
private protected override bool CompilesItself => true;
506+
507+
private protected override void CompileNode(Compiler compiler)
508+
{
509+
Argument.InnerCompile(compiler);
510+
compiler.Instructions.Add(new(InstructionType.CALL_SI));
511+
}
512+
}
513+
514+
public partial record Cif
515+
{
516+
private protected override bool CompilesItself => true;
517+
518+
private protected override void CompileNode(Compiler compiler)
519+
{
520+
Argument.InnerCompile(compiler);
521+
compiler.Instructions.Add(new(InstructionType.CALL_CI));
522+
}
523+
}
524+
525+
public partial record Shif
526+
{
527+
private protected override bool CompilesItself => true;
528+
529+
private protected override void CompileNode(Compiler compiler)
530+
{
531+
Argument.InnerCompile(compiler);
532+
compiler.Instructions.Add(new(InstructionType.CALL_SHI));
533+
}
534+
}
535+
536+
public partial record Chif
537+
{
538+
private protected override bool CompilesItself => true;
539+
540+
private protected override void CompileNode(Compiler compiler)
541+
{
542+
Argument.InnerCompile(compiler);
543+
compiler.Instructions.Add(new(InstructionType.CALL_CHI));
544+
}
545+
}
445546
}
446547
}
447548

‎Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,15 @@ internal enum InstructionType
4040
CALL_FLOOR,
4141
CALL_CEIL,
4242
CALL_ROUND,
43+
CALL_ERF,
44+
CALL_ERFC,
45+
CALL_ERFI,
46+
CALL_EI,
47+
CALL_LI,
48+
CALL_SI,
49+
CALL_CI,
50+
CALL_SHI,
51+
CALL_CHI,
4352

4453
// 2-arg functions
4554
CALL_SUM = 100,
@@ -219,6 +228,35 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value
219228
case InstructionType.CALL_ARCCOSECANT:
220229
stack.Push(Core.Compilation.ComplexBranches.Arccosecant(stack.Pop()));
221230
break;
231+
// In double precision, as the kernels define them.
232+
// https://github.com/asc-community/AngouriMath/issues/1607
233+
case InstructionType.CALL_ERF:
234+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Erf(stack.Pop()));
235+
break;
236+
case InstructionType.CALL_ERFC:
237+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Erfc(stack.Pop()));
238+
break;
239+
case InstructionType.CALL_ERFI:
240+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Erfi(stack.Pop()));
241+
break;
242+
case InstructionType.CALL_EI:
243+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Ei(stack.Pop()));
244+
break;
245+
case InstructionType.CALL_LI:
246+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Li(stack.Pop()));
247+
break;
248+
case InstructionType.CALL_SI:
249+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Si(stack.Pop()));
250+
break;
251+
case InstructionType.CALL_CI:
252+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Ci(stack.Pop()));
253+
break;
254+
case InstructionType.CALL_SHI:
255+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Shi(stack.Pop()));
256+
break;
257+
case InstructionType.CALL_CHI:
258+
stack.Push(AngouriMath.Numerics.SpecialFunctions.Chi(stack.Pop()));
259+
break;
222260
case InstructionType.CALL_FACTORIAL:
223261
// https://stackoverflow.com/a/15454784/5429648
224262
const int g = 7;

‎Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -137,6 +137,18 @@ Arccosecantf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef
137137
Floorf => Expression.Call(GetDef("Floor", 1, e.Type), e),
138138
Ceilf => Expression.Call(GetDef("Ceil", 1, e.Type), e),
139139
Roundf => Expression.Call(GetDef("Round", 1, e.Type), e),
140+
141+
// The special functions, widened to double first as the trigonometry is.
142+
// https://github.com/asc-community/AngouriMath/issues/1607
143+
Erff when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Erf", 1, newE.Type), newE),
144+
Erfcf when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Erfc", 1, newE.Type), newE),
145+
Erfif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Erfi", 1, newE.Type), newE),
146+
Eif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Ei", 1, newE.Type), newE),
147+
Lif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Li", 1, newE.Type), newE),
148+
Sif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Si", 1, newE.Type), newE),
149+
Cif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Ci", 1, newE.Type), newE),
150+
Shif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Shi", 1, newE.Type), newE),
151+
Chif when ShouldBeAtLeastDouble(e) is var newE => Expression.Call(GetDef("Chi", 1, newE.Type), newE),
140152

141153
Notf => Expression.Not(e),
142154

‎Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -288,6 +288,27 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a)
288288
public static int Min(int a, int b) => Math.Min(a, b);
289289
public static BigInteger Min(BigInteger a, BigInteger b) => BigInteger.Min(a, b);
290290

291+
// The special functions, in double precision as the kernels define them.
292+
// https://github.com/asc-community/AngouriMath/issues/1607
293+
public static System.Numerics.Complex Erf(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erf(a);
294+
public static double Erf(double a) => AngouriMath.Numerics.SpecialFunctions.Erf(a);
295+
public static System.Numerics.Complex Erfc(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erfc(a);
296+
public static double Erfc(double a) => AngouriMath.Numerics.SpecialFunctions.Erfc(a);
297+
public static System.Numerics.Complex Erfi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Erfi(a);
298+
public static double Erfi(double a) => AngouriMath.Numerics.SpecialFunctions.Erfi(a);
299+
public static System.Numerics.Complex Ei(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Ei(a);
300+
public static double Ei(double a) => AngouriMath.Numerics.SpecialFunctions.Ei(a);
301+
public static System.Numerics.Complex Li(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Li(a);
302+
public static double Li(double a) => AngouriMath.Numerics.SpecialFunctions.Li(a);
303+
public static System.Numerics.Complex Si(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Si(a);
304+
public static double Si(double a) => AngouriMath.Numerics.SpecialFunctions.Si(a);
305+
public static System.Numerics.Complex Ci(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Ci(a);
306+
public static double Ci(double a) => AngouriMath.Numerics.SpecialFunctions.Ci(a);
307+
public static System.Numerics.Complex Shi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Shi(a);
308+
public static double Shi(double a) => AngouriMath.Numerics.SpecialFunctions.Shi(a);
309+
public static System.Numerics.Complex Chi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Chi(a);
310+
public static double Chi(double a) => AngouriMath.Numerics.SpecialFunctions.Chi(a);
311+
291312
/// <summary>
292313
/// True where <c>System.Numerics.Complex.IsNaN</c> is, spelled out because that overload is
293314
/// .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)
438459
{ ("Round", 1, typeof(long)), Def<long>(a => Round(a)) },
439460
{ ("Round", 1, typeof(int)), Def<int>(a => Round(a)) },
440461
{ ("Round", 1, typeof(BigInteger)), Def<BigInteger>(a => Round(a)) },
462+
{ ("Erf", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Erf(a)) },
463+
{ ("Erf", 1, typeof(double)), Def<double>(a => Erf(a)) },
464+
{ ("Erfc", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Erfc(a)) },
465+
{ ("Erfc", 1, typeof(double)), Def<double>(a => Erfc(a)) },
466+
{ ("Erfi", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Erfi(a)) },
467+
{ ("Erfi", 1, typeof(double)), Def<double>(a => Erfi(a)) },
468+
{ ("Ei", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Ei(a)) },
469+
{ ("Ei", 1, typeof(double)), Def<double>(a => Ei(a)) },
470+
{ ("Li", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Li(a)) },
471+
{ ("Li", 1, typeof(double)), Def<double>(a => Li(a)) },
472+
{ ("Si", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Si(a)) },
473+
{ ("Si", 1, typeof(double)), Def<double>(a => Si(a)) },
474+
{ ("Ci", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Ci(a)) },
475+
{ ("Ci", 1, typeof(double)), Def<double>(a => Ci(a)) },
476+
{ ("Shi", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Shi(a)) },
477+
{ ("Shi", 1, typeof(double)), Def<double>(a => Shi(a)) },
478+
{ ("Chi", 1, typeof(System.Numerics.Complex)), Def<System.Numerics.Complex>(a => Chi(a)) },
479+
{ ("Chi", 1, typeof(double)), Def<double>(a => Chi(a)) },
441480
};
442481

443482
// An expression tree is how C# spells "the MethodInfo of this method" without a string:

‎Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,9 @@
99
"Cot", "Sec", "Csc", "Acot", "Asec", "Acsc",
1010
"Abs", "Sgn", "Floor", "Ceil", "Round" };
1111
var binary = new[]{ "Log", "Pow", "Max", "Min" };
12+
// Defined for Complex and double alone: an integer or a float is widened to double first, as
13+
// for the trigonometry.
14+
var special = new[]{ "Erf", "Erfc", "Erfi", "Ei", "Li", "Si", "Ci", "Shi", "Chi" };
1215
var types = new[]{ "System.Numerics.Complex", "double", "float", "long", "int", "BigInteger" };
1316
#>
1417
//
@@ -188,6 +191,13 @@ namespace AngouriMath.Core.Compilation.IntoLinq
188191
public static int Min(int a, int b) => Math.Min(a, b);
189192
public static BigInteger Min(BigInteger a, BigInteger b) => BigInteger.Min(a, b);
190193

194+
// The special functions, in double precision as the kernels define them.
195+
// https://github.com/asc-community/AngouriMath/issues/1607
196+
<# foreach (var name in special) { #>
197+
public static System.Numerics.Complex <#= name #>(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.<#= name #>(a);
198+
public static double <#= name #>(double a) => AngouriMath.Numerics.SpecialFunctions.<#= name #>(a);
199+
<# } #>
200+
191201
/// <summary>
192202
/// True where <c>System.Numerics.Complex.IsNaN</c> is, spelled out because that overload is
193203
/// .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
221231
<# foreach (var type in types) { #>
222232
{ ("<#= name #>", 1, typeof(<#= type #>)), Def<<#= type #>>(a => <#= name #>(a)) },
223233
<# } #>
234+
<# } #>
235+
<# foreach (var name in special) { #>
236+
<# foreach (var type in new[]{ "System.Numerics.Complex", "double" }) { #>
237+
{ ("<#= name #>", 1, typeof(<#= type #>)), Def<<#= type #>>(a => <#= name #>(a)) },
238+
<# } #>
224239
<# } #>
225240
};
226241

0 commit comments

Comments
 (0)