Skip to content

The sine, the cosine and the arctangent are series in the same fixed point as the logarithm and the exponential #1347

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Sep 15, 2026
Merged

The sine, the cosine and the arctangent are series in the same fixed point as the logarithm and the exponential #1347

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 26 additions & 0 deletions Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,6 +256,32 @@ cent — its measured run-to-run spread on mean time is up to 51.8%. A move of 8
outside that band by a wide margin and agrees in sign and rough size with the allocation column
beside it. The small rows from the same run are still not worth reading, and are not quoted.

## The 2055th: the sine, cosine and arctangent in the same fixed point

The 2050th made the trigonometric functions argument reduction and short series, in
`EDecimal`; the 2054th found that a series in `EDecimal` pays two to three microseconds a term
for the alignment in each addition and moved the logarithm and the exponential to fixed point.
The sine, the cosine and the arctangent are there too now: the reduction to `[-pi, pi]` is an
integer division by `2 pi` in fixed point, the folds and halvings are subtractions and shifts,
the square root in the arctangent's halving `x/(1 + sqrt(1 + x^2))` is the integer square root
of the value shifted up, exactly floored, and the halvings come back as a shift. The arcsine
and arccosine are the arctangent's, and the tangent the quotient. The hundred-digit values
are the ones pinned before, to the last digit.

| benchmark | 2054th | 2055th | allocation | time |
|---|--:|--:|--:|--:|
| `EvalTrig` | 1,192,089 | **260,728** | **−78.1%** | 552 → **153 µs** |
| `EvalTrigPrecise` | 7,935,226 | **2,525,707** | **−68.2%** | 13.8 → **2.16 ms** |
| `EvalTranscendentalFresh` | 423,032 | **128,664** | **−69.6%** | 245 → **83 µs** |
| `SolveEasy` | 4,425,194 | **1,675,487** | **−62.1%** | 2.88 → **1.19 ms** |
| every other entry | | | within 1.5% | |

Bytes allocated per call, same machine, both columns by the gate in one session; the gate's
baseline was taken from this run. Alone, at a hundred digits: `sin(x)` is 45 µs where it was
176, `arctan(x)` 46 where it was 202, `arcsin(x)` 59 where it was 238; at five hundred digits
`arctan(x)` is 1.2 ms where it was 9.3. Against the 2049th, before any of this, `EvalTrig` is
709 → 153 µs and `EvalTrigPrecise` 22.5 → 2.2 ms.

## The 2054th: the logarithm and the exponential as series in fixed point

`EvalTranscendentalFresh` -- `ln(1 + x^2) arctan(x)` at a fresh point -- was 1.3 ms after the
Expand Down
129 changes: 64 additions & 65 deletions Sources/AngouriMath/Functions/InternalAMExtensions.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,9 @@ public static ConstantCache Lookup(EContext context)
FixedOne = EInteger.One.ShiftLeft(FixedBits);
FiveToFixedBits = EInteger.FromInt32(5).Pow(FixedBits);
FixedSqrt2 = ToFixed(EDecimal.FromString("1.41421356237309504880168872420969807856967187537694"), FixedBits);
FixedPi = ToFixed(Pi, FixedBits);
FixedTwoPi = FixedPi.ShiftLeft(1);
FixedHalfPi = FixedPi.ShiftRight(1);
// ln 2 = 2 artanh(1/3), and ln 10 = ln 8 + ln 5/4 = 3 ln 2 + 2 artanh(1/9).
FixedLn2 = TwiceArtanh(FixedOne.Divide(3), FixedBits);
FixedLn10 = FixedLn2.Multiply(3).Add(TwiceArtanh(FixedOne.Divide(9), FixedBits));
Expand All @@ -219,6 +222,12 @@ public static ConstantCache Lookup(EContext context)
public EInteger FiveToFixedBits { get; }
/// <summary>The square root of 2, in fixed point, to fifty digits</summary>
public EInteger FixedSqrt2 { get; }
/// <summary>Pi in fixed point, to the context's digits</summary>
public EInteger FixedPi { get; }
/// <summary>2 pi in fixed point</summary>
public EInteger FixedTwoPi { get; }
/// <summary>Pi / 2 in fixed point</summary>
public EInteger FixedHalfPi { get; }
/// <summary>The natural logarithm of 2, in fixed point</summary>
public EInteger FixedLn2 { get; }
/// <summary>The natural logarithm of 10, in fixed point</summary>
Expand Down Expand Up @@ -406,69 +415,68 @@ internal static (EDecimal Sin, EDecimal Cos) SinAndCos(EDecimal x, EContext cont
{
var working = WithGuardDigits(context, 8);
var consts = ConstantCache.Lookup(working);
var bits = consts.FixedBits;
var one = consts.FixedOne;

//truncating to [-2*PI;2*PI]
TruncateToPeriodicInterval(ref x, consts, working);
// now x in (-2pi, 2pi); into [-pi, pi], with the signs the shift costs
// In fixed point from here: the reduction is integer arithmetic on pi's digits,
// and a halving is a shift.
var fixedX = ToFixed(x, bits);
// Into [-pi, pi]: the nearest multiple of 2 pi off.
var shifted = fixedX.Add(consts.FixedPi);
var turns = shifted.Divide(consts.FixedTwoPi);
if (shifted.Sign < 0 && !shifted.Remainder(consts.FixedTwoPi).IsZero)
turns = turns.Subtract(1); // the division truncates, and this floors
if (!turns.IsZero)
fixedX = fixedX.Subtract(consts.FixedTwoPi.Multiply(turns));
// Into [-pi/2, pi/2], with the sign the fold costs the cosine.
var negateCos = false;
if (x.GreaterThan(consts.Pi))
{
x = x.Subtract(consts.TwoPi, working);
}
else if (x.LessThan(-consts.Pi))
{
x = x.Add(consts.TwoPi, working);
}
if (x.GreaterThan(consts.HalfPi))
if (fixedX.CompareTo(consts.FixedHalfPi) > 0)
{
// sin(x) = sin(pi - x), cos(x) = -cos(pi - x)
x = consts.Pi.Subtract(x, working);
fixedX = consts.FixedPi.Subtract(fixedX);
negateCos = true;
}
else if (x.LessThan(-consts.HalfPi))
else if (fixedX.CompareTo(consts.FixedHalfPi.Negate()) < 0)
{
// sin(x) = sin(-pi - x), cos(x) = -cos(-pi - x)
x = consts.Pi.Negate().Subtract(x, working);
fixedX = consts.FixedPi.Negate().Subtract(fixedX);
negateCos = true;
}

// Halve until |x| < 1/20.
var halvings = 0;
var twentieth = EDecimal.Create(5, -2);
while (x.Abs().GreaterThan(twentieth))
var twentieth = one.Divide(20);
while (fixedX.Abs().CompareTo(twentieth) > 0)
{
x = x.Multiply(consts.Half, working);
fixedX = fixedX.Sign < 0 ? fixedX.Negate().ShiftRight(1).Negate() : fixedX.ShiftRight(1);
halvings++;
}

var xx = x.Multiply(x, working);
var square = MultiplyFixed(fixedX, fixedX, bits);
// sin: x - x^3/3! + ...; cos: 1 - x^2/2! + ...
var sinTerm = x;
var sin = x;
var cosTerm = EDecimal.One;
var cos = EDecimal.One;
for (var i = 1; i < 200; i++)
var sinTerm = fixedX;
var sin = fixedX;
var cosTerm = one;
var cos = one;
for (var i = 1; i < 100000; i++)
{
cosTerm = -cosTerm.Multiply(xx, working).Divide((2 * i - 1) * (2 * i), working);
sinTerm = -sinTerm.Multiply(xx, working).Divide((2 * i) * (2 * i + 1), working);
var nextCos = cos.Add(cosTerm, working);
var nextSin = sin.Add(sinTerm, working);
var settled = nextCos.Equals(cos) && nextSin.Equals(sin);
cos = nextCos;
sin = nextSin;
if (settled)
cosTerm = MultiplyFixed(cosTerm, square, bits).Negate().Divide((2 * i - 1) * (2 * i));
sinTerm = MultiplyFixed(sinTerm, square, bits).Negate().Divide((2 * i) * (2 * i + 1));
if (cosTerm.IsZero && sinTerm.IsZero)
break;
cos = cos.Add(cosTerm);
sin = sin.Add(sinTerm);
}
for (var i = 0; i < halvings; i++)
{
// sin(2y) = 2 sin(y) cos(y), cos(2y) = 2 cos(y)^2 - 1
var doubledSin = sin.Multiply(cos, working).Multiply(2, working);
cos = cos.Multiply(cos, working).Multiply(2, working).Decrement();
var doubledSin = MultiplyFixed(sin, cos, bits).ShiftLeft(1);
cos = MultiplyFixed(cos, cos, bits).ShiftLeft(1).Subtract(one);
sin = doubledSin;
}
if (negateCos)
cos = -cos;
return (sin.RoundToPrecision(context), cos.RoundToPrecision(context));
cos = cos.Negate();
return (FromFixed(sin, consts, working).RoundToPrecision(context), FromFixed(cos, consts, working).RoundToPrecision(context));
}

/// <summary>Analogy of <see cref="Math.Tan(double)"/></summary>
Expand Down Expand Up @@ -533,21 +541,6 @@ public static EDecimal Arccosecant(this EDecimal x, EContext context)
}

/// <summary>Truncates <paramref name="x"/> to [-2*<see cref="Math.PI"/>, 2*<see cref="Math.PI"/>] </summary>
private static void TruncateToPeriodicInterval(ref EDecimal x, ConstantCache consts, EContext context)
{
while (x.GreaterThanOrEquals(consts.TwoPi))
{
var divide = x.Divide(consts.TwoPi, context).Floor().Abs();
x = divide.MultiplyAndAdd(-consts.TwoPi, x, context);
}

while (x.LessThanOrEquals(-consts.TwoPi))
{
var divide = x.Divide(consts.TwoPi, context).Floor().Abs();
x = divide.MultiplyAndAdd(consts.TwoPi, x, context);
}
}

/// <summary>Analogy of <see cref="Math.Asin(double)"/></summary>
/// <remarks>
/// <c>arcsin(x) = arctan(x/sqrt(1 - x^2))</c>, with <c>1 - x^2</c> as <c>(1 - x)(1 + x)</c>
Expand Down Expand Up @@ -598,27 +591,33 @@ public static EDecimal Arctan(this EDecimal x, EContext context)
if (x.GreaterThan(EDecimal.One))
return ConstantCache.Lookup(working).HalfPi.Subtract(Arctan(EDecimal.One.Divide(x, working), working), working).RoundToPrecision(context);

// In fixed point from here: the square root of a fixed-point number is the
// integer square root of it shifted up, exactly floored, and the halvings come
// back as a shift.
var workingConsts = ConstantCache.Lookup(working);
var bits = workingConsts.FixedBits;
var one = workingConsts.FixedOne;
var fixedX = ToFixed(x, bits);
var halvings = 0;
var twentieth = EDecimal.Create(5, -2);
while (x.GreaterThan(twentieth))
var twentieth = one.Divide(20);
while (fixedX.CompareTo(twentieth) > 0)
{
x = x.Divide(EDecimal.One.Add(x.MultiplyAndAdd(x, EDecimal.One, working).Sqrt(working), working), working);
var root = one.Add(MultiplyFixed(fixedX, fixedX, bits)).ShiftLeft(bits).Sqrt();
fixedX = fixedX.ShiftLeft(bits).Divide(one.Add(root));
halvings++;
}
var xx = x.Multiply(x, working);
var power = x;
var sum = x;
for (var i = 1; i < 400; i++)
var square = MultiplyFixed(fixedX, fixedX, bits);
var power = fixedX;
var sum = fixedX;
for (var i = 1; i < 100000; i++)
{
power = -power.Multiply(xx, working);
var next = sum.Add(power.Divide(2 * i + 1, working), working);
if (next.Equals(sum))
power = MultiplyFixed(power, square, bits).Negate();
var term = power.Divide(2 * i + 1);
if (term.IsZero)
break;
sum = next;
sum = sum.Add(term);
}
if (halvings > 0)
sum = sum.Multiply(EDecimal.FromInt32(1 << halvings), working);
return sum.RoundToPrecision(context);
return FromFixed(sum.ShiftLeft(halvings), workingConsts, working).RoundToPrecision(context);
}
/// <summary>Analogy of <see cref="Math.Acos(double)"/></summary>
public static EDecimal Acos(this EDecimal x, EContext context)
Expand Down
64 changes: 32 additions & 32 deletions Sources/Tests/DotnetBenchmark/performance-baseline.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,94 +1,94 @@
{
"comment": "Allocated bytes and mean nanoseconds per operation for the popular use cases of https://github.com/asc-community/AngouriMath/issues/746. Checked by PerformanceGate.cs; see Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md for when updating it is legitimate.",
"benchmark": "CommonFunctionsInterVersion",
"commit": "bfb60e9cdcdc855c13b5ac58001db0ec98e3f3b2",
"commit": "95082e2f814874b99f152fbb1072fc977735337f",
"measuredOn": "2026-09-15",
"runtime": ".NET 10.0.10",
"machine": "Ubuntu 26.04 LTS, X64, 8 logical cores",
"cases": {
"CompileEasy": {
"allocatedBytes": 11023,
"meanNanoseconds": 184521.07338867188
"allocatedBytes": 11031,
"meanNanoseconds": 184131.37203543526
},
"CompileHard": {
"allocatedBytes": 20737,
"meanNanoseconds": 309526.83076171874
"allocatedBytes": 20434,
"meanNanoseconds": 310381.04733072914
},
"Derivate": {
"allocatedBytes": 53099,
"meanNanoseconds": 14289.474793570382
"allocatedBytes": 53111,
"meanNanoseconds": 14270.049594116212
},
"EvalEasy": {
"allocatedBytes": 0,
"meanNanoseconds": 1.8900332285889558
"meanNanoseconds": 1.883551713079214
},
"EvalPolynomialFresh": {
"allocatedBytes": 12328,
"meanNanoseconds": 4151.610092798869
"meanNanoseconds": 4164.772733835073
},
"EvalPolynomialFresh15Digits": {
"allocatedBytes": 11656,
"meanNanoseconds": 4141.208742959158
"meanNanoseconds": 4157.980944497244
},
"EvalTranscendentalFresh": {
"allocatedBytes": 423032,
"meanNanoseconds": 245377.50509207588
"allocatedBytes": 128664,
"meanNanoseconds": 83463.06102643695
},
"EvalTrig": {
"allocatedBytes": 1192089,
"meanNanoseconds": 551794.5797293527
"allocatedBytes": 260728,
"meanNanoseconds": 153426.1548339844
},
"EvalTrigPrecise": {
"allocatedBytes": 7935226,
"meanNanoseconds": 13781485.40625
"allocatedBytes": 2525707,
"meanNanoseconds": 2164306.404296875
},
"ParseEasy": {
"allocatedBytes": 18125,
"meanNanoseconds": 5981.867161090558
"allocatedBytes": 18104,
"meanNanoseconds": 5932.770784231333
},
"ParseHard": {
"allocatedBytes": 3555057,
"meanNanoseconds": 1367037.3726283482
"meanNanoseconds": 1392708.7341145833
},
"RunEasy": {
"allocatedBytes": 0,
"meanNanoseconds": 20.056596242464504
"meanNanoseconds": 20.095689982175827
},
"RunHard": {
"allocatedBytes": 0,
"meanNanoseconds": 295.51566811970304
"meanNanoseconds": 295.40826177597046
},
"RunMedium": {
"allocatedBytes": 0,
"meanNanoseconds": 163.17946696281433
"meanNanoseconds": 162.8344373519604
},
"SimplifyEasy": {
"allocatedBytes": 78578,
"meanNanoseconds": 45092.335445149736
"meanNanoseconds": 44253.979728190105
},
"SimplifyHard": {
"allocatedBytes": 189726120,
"meanNanoseconds": 109370614.94736843
"allocatedBytes": 189604392,
"meanNanoseconds": 111491250.2413793
},
"SolveEasy": {
"allocatedBytes": 4425194,
"meanNanoseconds": 2884867.8782552085
"allocatedBytes": 1675487,
"meanNanoseconds": 1185599.8670479911
},
"SolveEasyMedium": {
"allocatedBytes": 65184,
"meanNanoseconds": 18597.25531768799
"meanNanoseconds": 18840.902629307337
},
"SolveHard": {
"allocatedBytes": 11706656,
"meanNanoseconds": 64677752.72
"meanNanoseconds": 63240162.666666664
},
"SolveMedium": {
"allocatedBytes": 446445,
"meanNanoseconds": 266309.50540597097
"allocatedBytes": 446414,
"meanNanoseconds": 259625.15890066963
},
"SolveMediumHard": {
"allocatedBytes": 1367798,
"meanNanoseconds": 8561109.407291668
"meanNanoseconds": 8776522.877083333
}
},
"ungated": {
Expand Down
Loading