From fcb654bb849f007f071c0c373e4c7a3ed62d361d Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 22 Sep 2026 01:53:21 +0000 Subject: [PATCH] Implement ITrigonometricFunctions, plus Atan2 [minor] MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds the circular trigonometric functions and their inverses, none routing through double, completing #78's fourth interface surface. Sin/Cos/SinCos reduce modulo π/2 into [-π/4, π/4] with an octant index, so one halving-and-doubling kernel serves both; the reduction reads PiTo at the width the argument demands (per #79) rather than the capped Pi property, which is what keeps Sin of a large angle meaningful. Tan is sin/cos from one reduction. Atan reduces by the half-angle identity until the argument is small, then Taylor; Asin/Acos build on Atan and Sqrt, with the ±1 endpoints special-cased around the division by zero. The half-turn family reduces on the argument before multiplying by π, so SinPi(1e20) is well-defined where Sin(1e20·π) is not. Atan2 is a bespoke static rather than an interface member: it lives on IFloatingPointIeee754, which PreciseNumber does not implement because it has no NaN or infinity to give the interface's edge cases meaning. Tests pin published digits for the standard angles, the Pythagorean and double-angle identities across a sweep, the SinCos/Sin+Cos agreement, the inverse round trips, Atan2 in every quadrant and on the axes, and a 120-digit Sin(1000000) reference — the last of which cannot be produced from a short π, so it fails on an implementation that reduces against one. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01By7NnPN7STCZqeAftJ1BmH --- .../TrigonometryBenchmarks.cs | 120 +++ .../PreciseNumberTrigonometryTests.cs | 299 ++++++ PreciseNumber/PreciseNumber.Trigonometry.cs | 872 ++++++++++++++++++ README.md | 15 +- 4 files changed, 1304 insertions(+), 2 deletions(-) create mode 100644 PreciseNumber.Benchmarks/TrigonometryBenchmarks.cs create mode 100644 PreciseNumber.Test/PreciseNumberTrigonometryTests.cs create mode 100644 PreciseNumber/PreciseNumber.Trigonometry.cs diff --git a/PreciseNumber.Benchmarks/TrigonometryBenchmarks.cs b/PreciseNumber.Benchmarks/TrigonometryBenchmarks.cs new file mode 100644 index 0000000..22abc9d --- /dev/null +++ b/PreciseNumber.Benchmarks/TrigonometryBenchmarks.cs @@ -0,0 +1,120 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Benchmarks; + +using BenchmarkDotNet.Attributes; + +/// +/// Measures the circular trigonometric functions and their inverses. +/// +/// +/// These have no predecessor in the type — there was no trigonometry before — +/// so the and are a floor rather +/// than a thing replaced: they answer in about fifteen correct digits whatever the Digits +/// axis says. +/// +/// The Digits axis drives both the operand and the digits asked of the answer, since the +/// default precision follows the operand. The cost of a sine grows with the digit count on two +/// fronts: the series gains terms, and every term is a wider . +/// Read it against 's division, which sits inside every loop here. +/// +/// +/// against and +/// is the measurement that justifies the combined member: it does one argument reduction where the +/// two separate calls do two, so it should cost noticeably less than their sum. +/// +/// +/// is the case the wide π pays for. Its reduction reads π to the +/// width of the argument on top of the answer, so it should sit above by +/// the cost of that wider constant and its multiplication, and the gap should widen with the +/// Digits axis. +/// +/// +[MemoryDiagnoser] +public class TrigonometryBenchmarks +{ + private PreciseNumber angle = PreciseNumber.Zero; + private PreciseNumber largeAngle = PreciseNumber.Zero; + private PreciseNumber unitInterval = PreciseNumber.Zero; + private PreciseNumber tangent = PreciseNumber.Zero; + private PreciseNumber ordinate = PreciseNumber.Zero; + private PreciseNumber abscissa = PreciseNumber.Zero; + private double angleAsDouble; + + /// + /// Gets or sets the number of significant digits in the operand, and so in the answer. + /// + [Params(8, 30, 200)] + public int Digits { get; set; } + + /// + /// Prepares the operands. + /// + [GlobalSetup] + public void Setup() + { + // A little over one radian, so the reduction into an octant and the series both run. + angle = Operands.Number(Digits, -(Digits - 1)); + + // Several thousand radians, so the reduction reads π wider than the answer. + largeAngle = Operands.Number(Digits, -(Digits - 4), offset: 7); + + // In (0, 1), the domain of the inverse sine. + unitInterval = Operands.Number(Digits, -Digits, offset: 13); + + tangent = Operands.Number(Digits, -(Digits - 1), offset: 29); + ordinate = Operands.Number(Digits, -(Digits - 1), offset: 5); + abscissa = -Operands.Number(Digits, -(Digits - 1), offset: 17); + angleAsDouble = angle.To(); + } + + /// Takes the sine of an angle. + /// The result. + [Benchmark] + public PreciseNumber SinOfAnAngle() => PreciseNumber.Sin(angle); + + /// Takes the cosine of an angle. + /// The result. + [Benchmark] + public PreciseNumber CosOfAnAngle() => PreciseNumber.Cos(angle); + + /// Takes the sine and cosine of an angle from one reduction. + /// The result. + [Benchmark] + public (PreciseNumber Sin, PreciseNumber Cos) SinCosOfAnAngle() => PreciseNumber.SinCos(angle); + + /// Takes the tangent of an angle, which is one reduction and a division. + /// The result. + [Benchmark] + public PreciseNumber TanOfAnAngle() => PreciseNumber.Tan(angle); + + /// Takes the sine of a large angle, whose reduction reads a wide π. + /// The result. + [Benchmark] + public PreciseNumber SinOfALargeAngle() => PreciseNumber.Sin(largeAngle); + + /// Takes the arc tangent, whose half-angle reduction precedes the series. + /// The result. + [Benchmark] + public PreciseNumber AtanOfAValue() => PreciseNumber.Atan(tangent); + + /// Takes the arc sine, which is a square root and an arc tangent. + /// The result. + [Benchmark] + public PreciseNumber AsinOfAValue() => PreciseNumber.Asin(unitInterval); + + /// Takes the two-argument arc tangent, with quadrant dispatch. + /// The result. + [Benchmark] + public PreciseNumber Atan2OfAPoint() => PreciseNumber.Atan2(ordinate, abscissa); + + /// The sine, as a floor on the measurement. + /// The result. + [Benchmark(Baseline = true)] + public double DoubleSinBaseline() => Math.Sin(angleAsDouble); + + /// The two-argument arc tangent, as a floor on the measurement. + /// The result. + [Benchmark] + public double DoubleAtan2Baseline() => Math.Atan2(angleAsDouble, -angleAsDouble); +} diff --git a/PreciseNumber.Test/PreciseNumberTrigonometryTests.cs b/PreciseNumber.Test/PreciseNumberTrigonometryTests.cs new file mode 100644 index 0000000..ee4a9cc --- /dev/null +++ b/PreciseNumber.Test/PreciseNumberTrigonometryTests.cs @@ -0,0 +1,299 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber.Test; + +using System.Globalization; +using System.Numerics; + +/// +/// Covers on , together with +/// the bespoke . +/// +/// +/// The digit-for-digit assertions carry published values rather than values this library produced, +/// so a series that agrees with itself but not with mathematics still fails. +/// +/// Several assertions here exist to pin the type's precision claim against the constant it reduces +/// by. in particular carries a reference +/// far wider than a can hold: reducing an argument of magnitude 10^6 to +/// that many digits is only possible against a π of at least that width, so the same test that +/// confirms the answer also confirms the reduction reads π wide rather than capped. +/// +/// +[TestClass] +public class PreciseNumberTrigonometryTests +{ + /// The first fifty significant digits of the sine, cosine and tangent of one radian. + private const string Sin1Digits = "84147098480789650665250232163029899962256306079837"; + private const string Cos1Digits = "54030230586813971740093660744297660373231042061792"; + private const string Tan1Digits = "15574077246549022305069748074583601730872507723815"; + + /// The first fifty significant digits of π/6 and π/4. + private const string PiOverSixDigits = "52359877559829887307710723054658381403286156656252"; + private const string PiOverFourDigits = "78539816339744830961566084581987572104929234984378"; + + /// π/6 and π/3 as decimals, for value comparisons where a trailing zero would be normalised away. + private const string PiOverSix = "0.52359877559829887307710723054658381403286156656252"; + private const string PiOverThree = "1.0471975511965977461542144610931676280657231331250"; + + /// The first fifty significant digits of √3 / 2, the cosine of π/6. + private const string RootThreeOverTwoDigits = "86602540378443864676372317075293618347140262690519"; + + /// + /// The sine of one million radians, to a reference far wider than a holds. + /// + /// + /// A hand-rolled argument reduction against a short π cannot reach these digits: a π of only the + /// twenty-six digits the type once carried leaves about twenty-one usable digits of a + /// 10^6 argument, so a reference this wide is a direct test of the constant behind the + /// reduction as much as of the series in front of it. + /// + private const string SinOneMillionDigits = + "-0.3499935021712929521176524867807714690614066053287162738570590546446412263954505050656668976688940081127331690567910649695709417663"; + + private static PreciseNumber Parse(string text) => + PreciseNumber.Parse(text, CultureInfo.InvariantCulture); + + private static string Digits(PreciseNumber value) => + value.Significand.ToString(CultureInfo.InvariantCulture); + + /// + /// Asserts that two values agree to a number of significant digits, comparing relative to the + /// expected magnitude so the assertion means the same thing at every exponent. + /// + private static void AssertAgreesTo(PreciseNumber expected, PreciseNumber actual, int digits, string message) + { + PreciseNumber difference = PreciseNumber.Abs(actual - expected); + PreciseNumber tolerance = PreciseNumber.Abs(expected) * Parse($"1E-{digits.ToString(CultureInfo.InvariantCulture)}"); + + Assert.IsTrue( + difference <= tolerance, + $"{message}: expected {expected}, got {actual}, which differs by {difference}"); + } + + /// + /// A spread of angles across several revolutions and both signs, including magnitudes whose + /// reduction depends on a π far wider than the answer. + /// + private static PreciseNumber[] AngleSweep() => + [ + Parse("0.3"), + Parse("-0.3"), + Parse("1"), + Parse("2.7"), + Parse("-5.5"), + Parse("3.14159"), + Parse("100"), + Parse("-1000"), + Parse("1000000"), + Parse("1E-40"), + ]; + + [TestMethod] + public void TestSinCosTanMatchPublishedDigits() + { + Assert.AreEqual(Sin1Digits, Digits(PreciseNumber.Sin(PreciseNumber.One, 50)), "Sin(1) is wrong"); + Assert.AreEqual(Cos1Digits, Digits(PreciseNumber.Cos(PreciseNumber.One, 50)), "Cos(1) is wrong"); + Assert.AreEqual(Tan1Digits, Digits(PreciseNumber.Tan(PreciseNumber.One, 50)), "Tan(1) is wrong"); + } + + [TestMethod] + public void TestSineOfKnownAnglesFromTheirRadianConstants() + { + // sin(π/6) = 1/2 exactly, and cos(π/6) = √3/2, checked from the angle rather than asserted. + Assert.AreEqual("0.5", PreciseNumber.Sin(Parse("0." + PiOverSixDigits), 48).ToString()); + Assert.AreEqual(RootThreeOverTwoDigits, Digits(PreciseNumber.Cos(Parse("0." + PiOverSixDigits), 50))); + } + + [TestMethod] + public void TestSinOfZeroAndCosOfZero() + { + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Sin(PreciseNumber.Zero)); + Assert.AreEqual(PreciseNumber.One, PreciseNumber.Cos(PreciseNumber.Zero)); + } + + [TestMethod] + public void TestSinReducesALargeArgumentAgainstAWidePi() + { + // This is the assertion that ties the trigonometry to the width of Pi. A reduction against a + // short constant loses roughly one digit per digit of the argument, so a 10^6 argument held + // to 120 digits could not be produced from a 26-digit Pi at all. + PreciseNumber sine = PreciseNumber.Sin(Parse("1000000"), 130); + AssertAgreesTo(Parse(SinOneMillionDigits), sine, 120, "Sin(1000000) did not match its wide reference"); + } + + [TestMethod] + public void TestPythagoreanIdentityHoldsAcrossTheSweep() + { + foreach (PreciseNumber angle in AngleSweep()) + { + (PreciseNumber sin, PreciseNumber cos) = PreciseNumber.SinCos(angle, 60); + PreciseNumber identity = PreciseNumber.Add( + PreciseNumber.Multiply(sin, sin), + PreciseNumber.Multiply(cos, cos)); + AssertAgreesTo(PreciseNumber.One, identity, 58, $"sin² + cos² ({angle}) was not one"); + } + } + + [TestMethod] + public void TestSinCosSharesOneReductionWithSinAndCos() + { + // The member exists so a caller pays for the reduction once; both halves must equal the + // separate calls exactly, or SinCos would be answering a different question than Sin and Cos. + foreach (PreciseNumber angle in AngleSweep()) + { + (PreciseNumber sin, PreciseNumber cos) = PreciseNumber.SinCos(angle, 60); + Assert.AreEqual(PreciseNumber.Sin(angle, 60), sin, $"SinCos({angle}).Sin disagreed with Sin"); + Assert.AreEqual(PreciseNumber.Cos(angle, 60), cos, $"SinCos({angle}).Cos disagreed with Cos"); + } + } + + [TestMethod] + public void TestTanIsSinOverCos() + { + foreach (PreciseNumber angle in AngleSweep()) + { + (PreciseNumber sin, PreciseNumber cos) = PreciseNumber.SinCos(angle, 60); + PreciseNumber expected = PreciseNumber.Divide(sin, cos, 50); + AssertAgreesTo(expected, PreciseNumber.Tan(angle, 50), 49, $"Tan({angle}) was not sin/cos"); + } + } + + [TestMethod] + public void TestSinPiIsWellDefinedAtAHugeArgument() + { + // SinPi(1e20) is an even multiple of π, so exactly zero; Sin(1e20 * Pi) is not defined at all, + // which is the whole reason the half-turn family reduces on the argument before multiplying. + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.SinPi(Parse("1E20"), 50)); + Assert.AreEqual(PreciseNumber.One, PreciseNumber.SinPi(Parse("0.5"), 50)); + Assert.AreEqual(PreciseNumber.NegativeOne, PreciseNumber.CosPi(PreciseNumber.One, 50)); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.SinPi(Parse("2"), 50)); + } + + [TestMethod] + public void TestSinPiAgreesWithSinOfTheAngleInRadians() + { + foreach (PreciseNumber halfTurns in new[] { Parse("0.1"), Parse("0.37"), Parse("-0.8"), Parse("1.25") }) + { + PreciseNumber radians = PreciseNumber.Multiply(halfTurns, PreciseNumber.PiTo(70)); + AssertAgreesTo(PreciseNumber.Sin(radians, 55), PreciseNumber.SinPi(halfTurns, 50), 49, $"SinPi({halfTurns})"); + AssertAgreesTo(PreciseNumber.Cos(radians, 55), PreciseNumber.CosPi(halfTurns, 50), 49, $"CosPi({halfTurns})"); + } + } + + [TestMethod] + public void TestAtanMatchesPiOverFour() => + Assert.AreEqual(PiOverFourDigits, Digits(PreciseNumber.Atan(PreciseNumber.One, 50))); + + [TestMethod] + public void TestAsinAndAcosMatchTheirKnownAngles() + { + Assert.AreEqual(PiOverSixDigits, Digits(PreciseNumber.Asin(Parse("0.5"), 50)), "Asin(0.5) is not π/6"); + AssertAgreesTo(Parse(PiOverSix), PreciseNumber.Asin(Parse("0.5"), 50), 49, "Asin(0.5) is not π/6"); + AssertAgreesTo(Parse(PiOverThree), PreciseNumber.Acos(Parse("0.5"), 50), 49, "Acos(0.5) is not π/3"); + } + + [TestMethod] + public void TestAsinAndAcosAtTheEndpoints() + { + // The endpoints are where asin's √(1 - x²) is zero, so they are special-cased rather than + // divided by zero. + AssertAgreesTo(PreciseNumber.Divide(PreciseNumber.PiTo(60), Parse("2"), 55), PreciseNumber.Asin(PreciseNumber.One, 50), 49, "Asin(1)"); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Acos(PreciseNumber.One, 50)); + AssertAgreesTo(PreciseNumber.PiTo(60), PreciseNumber.Acos(PreciseNumber.NegativeOne, 50), 49, "Acos(-1)"); + } + + [TestMethod] + public void TestInverseFunctionsRoundTripSineAndCosine() + { + foreach (PreciseNumber value in new[] { Parse("0.1"), Parse("-0.4"), Parse("0.9"), Parse("0.999") }) + { + AssertAgreesTo(value, PreciseNumber.Sin(PreciseNumber.Asin(value, 60), 60), 49, $"Sin(Asin({value}))"); + AssertAgreesTo(value, PreciseNumber.Cos(PreciseNumber.Acos(value, 60), 60), 49, $"Cos(Acos({value}))"); + AssertAgreesTo(value, PreciseNumber.Tan(PreciseNumber.Atan(value, 60), 60), 49, $"Tan(Atan({value}))"); + } + } + + [TestMethod] + public void TestAtanReducesALargeArgument() + { + // atan of a large value approaches π/2, and the half-angle reduction is what gets it there + // without a series that never converges. + AssertAgreesTo( + PreciseNumber.Divide(PreciseNumber.PiTo(60), Parse("2"), 55), + PreciseNumber.Atan(Parse("1E30"), 50), + 30, + "Atan(1E30) did not approach π/2"); + } + + [TestMethod] + public void TestInverseFunctionsRejectValuesOutsideTheirDomain() + { + // There is no NaN to return, so the domain is enforced. + Assert.ThrowsExactly(() => PreciseNumber.Asin(Parse("1.5"))); + Assert.ThrowsExactly(() => PreciseNumber.Acos(Parse("-2"))); + Assert.ThrowsExactly(() => PreciseNumber.Sin(PreciseNumber.One, 0)); + } + + [TestMethod] + public void TestAtan2PlacesTheAngleInEveryQuadrant() + { + PreciseNumber pi = PreciseNumber.PiTo(55); + PreciseNumber quarterPi = PreciseNumber.Divide(pi, Parse("4"), 55); + PreciseNumber threeQuarterPi = PreciseNumber.Multiply(Parse("3"), quarterPi); + + AssertAgreesTo(quarterPi, PreciseNumber.Atan2(PreciseNumber.One, PreciseNumber.One, 50), 49, "atan2(+, +)"); + AssertAgreesTo(threeQuarterPi, PreciseNumber.Atan2(PreciseNumber.One, PreciseNumber.NegativeOne, 50), 49, "atan2(+, -)"); + AssertAgreesTo(-threeQuarterPi, PreciseNumber.Atan2(PreciseNumber.NegativeOne, PreciseNumber.NegativeOne, 50), 49, "atan2(-, -)"); + AssertAgreesTo(-quarterPi, PreciseNumber.Atan2(PreciseNumber.NegativeOne, PreciseNumber.One, 50), 49, "atan2(-, +)"); + } + + [TestMethod] + public void TestAtan2OnTheAxes() + { + PreciseNumber pi = PreciseNumber.PiTo(55); + PreciseNumber halfPi = PreciseNumber.Divide(pi, Parse("2"), 55); + + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Atan2(PreciseNumber.Zero, PreciseNumber.Zero, 50)); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.Atan2(PreciseNumber.Zero, PreciseNumber.One, 50)); + AssertAgreesTo(halfPi, PreciseNumber.Atan2(PreciseNumber.One, PreciseNumber.Zero, 50), 49, "atan2(+, 0)"); + AssertAgreesTo(-halfPi, PreciseNumber.Atan2(PreciseNumber.NegativeOne, PreciseNumber.Zero, 50), 49, "atan2(-, 0)"); + AssertAgreesTo(pi, PreciseNumber.Atan2(PreciseNumber.Zero, PreciseNumber.NegativeOne, 50), 49, "atan2(0, -)"); + } + + [TestMethod] + public void TestAgreesWithDoublePrecisionAsACheapRegressionNet() + { + foreach (double sample in new[] { 0.3, 2.7, -5.5, 1.0, 0.75 }) + { + PreciseNumber argument = Parse(sample.ToString("R", CultureInfo.InvariantCulture)); + AssertAgreesTo(Parse(Math.Sin(sample).ToString("R", CultureInfo.InvariantCulture)), PreciseNumber.Sin(argument, 20), 14, $"Sin({sample})"); + AssertAgreesTo(Parse(Math.Cos(sample).ToString("R", CultureInfo.InvariantCulture)), PreciseNumber.Cos(argument, 20), 14, $"Cos({sample})"); + AssertAgreesTo(Parse(Math.Atan(sample).ToString("R", CultureInfo.InvariantCulture)), PreciseNumber.Atan(argument, 20), 14, $"Atan({sample})"); + } + + AssertAgreesTo( + Parse(Math.Atan2(1.0, -1.0).ToString("R", CultureInfo.InvariantCulture)), + PreciseNumber.Atan2(PreciseNumber.One, PreciseNumber.NegativeOne, 20), + 14, + "Atan2(1, -1)"); + } + + [TestMethod] + public void TestDegreeAndRadianConversionUseACorrectlyRoundedPi() + { + AssertAgreesTo(PreciseNumber.PiTo(55), PreciseNumber.DegreesToRadians(Parse("180"), 50), 49, "180° in radians is π"); + Assert.AreEqual("180", PreciseNumber.RadiansToDegrees(PreciseNumber.PiTo(55), 50).ToString()); + AssertAgreesTo(Parse("45"), PreciseNumber.RadiansToDegrees(PreciseNumber.Divide(PreciseNumber.PiTo(55), Parse("4"), 55), 50), 48, "π/4 in degrees is 45"); + } + + [TestMethod] + public void TestInverseHalfTurnFamilyReturnsHalfTurns() + { + // asin/acos/atan divided by π: the endpoints are exact. + Assert.AreEqual(Parse("0.5"), PreciseNumber.AsinPi(PreciseNumber.One, 50)); + Assert.AreEqual(PreciseNumber.Zero, PreciseNumber.AcosPi(PreciseNumber.One, 50)); + Assert.AreEqual(PreciseNumber.One, PreciseNumber.AcosPi(PreciseNumber.NegativeOne, 50)); + AssertAgreesTo(Parse("0.25"), PreciseNumber.AtanPi(PreciseNumber.One, 50), 49, "AtanPi(1) is a quarter turn"); + } +} diff --git a/PreciseNumber/PreciseNumber.Trigonometry.cs b/PreciseNumber/PreciseNumber.Trigonometry.cs new file mode 100644 index 0000000..1757f7c --- /dev/null +++ b/PreciseNumber/PreciseNumber.Trigonometry.cs @@ -0,0 +1,872 @@ +// Copyright (c) 2023-2026 ktsu-dev contributors + +namespace ktsu.PreciseNumber; + +using System; +using System.Numerics; + +/// +/// Circular trigonometric functions and their inverses, plus , +/// none of which route through . +/// +/// +/// The accuracy of a trigonometric function is the accuracy of the constant it reduces by. Reducing +/// an argument of magnitude 10^d to n correct digits needs π to roughly d + n +/// digits, so every reduction here reads at the width the argument demands +/// rather than the property — which is exactly why +/// of a large angle is meaningful at all. +/// +/// and reduce modulo π/2 +/// into [-π/4, π/4] with an octant index, so one kernel serves both and the series stays +/// short. exists to do that reduction once for a caller that +/// needs both, and both individual members read from it. +/// +/// +/// The half-turn family ( and the rest) reduces modulo two on the +/// argument before multiplying by π, so the multiplication never magnifies the argument and +/// no wide π is needed: SinPi(1e20) is well-defined where Sin(1e20 · π) is not. +/// +/// +/// has no NaN, so and +/// throw outside [-1, 1] where a would +/// return NaN and carry on, the same way does. +/// +/// +public readonly partial record struct PreciseNumber + : ITrigonometricFunctions +{ + /// + /// Digits computed past the ones the caller asked for, so that the digit the final rounding + /// decision is made on is itself correct. + /// + /// + /// The same width the exponential series uses, and for the same reason: a series makes one + /// rounding decision per term, and the doubling that follows argument halving compounds whatever + /// relative error it is handed. + /// + private const int TrigonometricGuardDigits = 10; + + /// + /// Times the sine/cosine kernel may halve its argument before summing. + /// + /// + /// The argument reaching the kernel is at most π/4, so six halvings always bring it under + /// ; the rest is slack for an argument that landed a little + /// outside the octant because the reduction rounded. + /// + private const int MaximumTrigonometricHalvings = 12; + + /// + /// Times Atan may apply its half-angle reduction before summing. + /// + /// + /// The first reduction brings any magnitude down to about one, and each after that roughly halves + /// the argument, so reaching from one takes about six more. The + /// bound is generous against an argument that starts just above the limit. + /// + private const int MaximumAtanReductions = 64; + + /// The message carried by the exception thrown for an inverse sine or cosine outside [-1, 1]. + private const string InverseDomainMessage = "The inverse sine and cosine are only defined on the interval [-1, 1]."; + + /// One hundred and eighty, the degrees in a half turn. + private static PreciseNumber OneEighty { get; } = new(0, 180); + + /// + /// The magnitude at or below which Atan sums its series directly rather than reducing first. + /// + /// + /// One sixty-fourth, matching the exponential's . Below it the + /// arctangent series converges in a handful of terms. + /// + private static PreciseNumber AtanReductionLimit { get; } = SeriesArgumentLimit; + + /// + /// Returns the sine of an angle in radians. + /// + /// The angle, in radians. + /// The sine of . + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. + /// + public static PreciseNumber Sin(PreciseNumber x) => + Sin(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the sine of an angle in radians, to a chosen number of significant digits. + /// + /// The angle, in radians. + /// The number of significant digits to produce. + /// The sine of . + /// Thrown when is less than one. + public static PreciseNumber Sin(PreciseNumber x, int significantDigits) => + SinCos(x, significantDigits).Sin; + + /// + /// Returns the cosine of an angle in radians. + /// + /// The angle, in radians. + /// The cosine of . + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. + /// + public static PreciseNumber Cos(PreciseNumber x) => + Cos(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the cosine of an angle in radians, to a chosen number of significant digits. + /// + /// The angle, in radians. + /// The number of significant digits to produce. + /// The cosine of . + /// Thrown when is less than one. + public static PreciseNumber Cos(PreciseNumber x, int significantDigits) => + SinCos(x, significantDigits).Cos; + + /// + /// Returns the sine and cosine of an angle in radians. + /// + /// The angle, in radians. + /// A tuple of the sine and cosine of . + /// + /// Both come from one argument reduction and one kernel, which is the reason to prefer this over a + /// separate and when both are + /// wanted. + /// + public static (PreciseNumber Sin, PreciseNumber Cos) SinCos(PreciseNumber x) => + SinCos(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the sine and cosine of an angle in radians, to a chosen number of significant digits. + /// + /// The angle, in radians. + /// The number of significant digits to produce. + /// A tuple of the sine and cosine of . + /// Thrown when is less than one. + /// + /// x = q · π/2 + r with q the nearest integer and r in [-π/4, π/4]. + /// The kernel evaluates the sine and cosine of r, and q mod four selects which, and + /// with which sign, becomes the sine and cosine of x. The π/2 the reduction + /// subtracts is read wide enough that the integer part it cancels was present in the constant + /// rather than invented. + /// + public static (PreciseNumber Sin, PreciseNumber Cos) SinCos(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return (Zero, One); + } + + int working = significantDigits + TrigonometricGuardDigits; + + // The subtraction below cancels the integer part of x / (π/2), so the constant has to carry + // that many digits past the answer or the remainder is only as good as what was left over. + int argumentDigits = IntegerDigitCount(x); + int reductionDigits = working + argumentDigits + TrigonometricGuardDigits; + + PreciseNumber piOverTwo = Divide(PiTo(reductionDigits), Two, reductionDigits); + BigInteger quadrant = RoundToNearestInteger(Divide(x, piOverTwo, reductionDigits)); + PreciseNumber remainder = Subtract(x, Multiply(new(0, quadrant), piOverTwo)) + .ReduceSignificance(working); + + (PreciseNumber sinRemainder, PreciseNumber cosRemainder) = SmallAngleSinCos(remainder, working); + return SelectOctant(quadrant, sinRemainder, cosRemainder, significantDigits); + } + + /// + /// Returns the tangent of an angle in radians. + /// + /// The angle, in radians. + /// The tangent of . + /// + /// Produced to the significant digits of , and never fewer than + /// . Use to choose + /// that precision. + /// + public static PreciseNumber Tan(PreciseNumber x) => + Tan(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the tangent of an angle in radians, to a chosen number of significant digits. + /// + /// The angle, in radians. + /// The number of significant digits to produce. + /// The tangent of . + /// Thrown when is less than one. + /// Thrown when the cosine of is zero. + /// + /// sin / cos from one reduction and a single division, rather than two independent series. + /// + public static PreciseNumber Tan(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + int working = significantDigits + TrigonometricGuardDigits; + (PreciseNumber sin, PreciseNumber cos) = SinCos(x, working); + return Divide(sin, cos, significantDigits); + } + + /// + /// Returns the arc sine of a value, in radians. + /// + /// The value, which must lie in [-1, 1]. + /// The angle in radians whose sine is , in [-π/2, π/2]. + /// Thrown when lies outside [-1, 1]. + public static PreciseNumber Asin(PreciseNumber x) => + Asin(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the arc sine of a value, in radians, to a chosen number of significant digits. + /// + /// The value, which must lie in [-1, 1]. + /// The number of significant digits to produce. + /// The angle in radians whose sine is , in [-π/2, π/2]. + /// + /// Thrown when lies outside [-1, 1], or when + /// is less than one. + /// + /// + /// asin x = atan( x / √(1 - x²) ), with asin(±1) = ±π/2 special-cased because the + /// division is by zero exactly there. + /// + public static PreciseNumber Asin(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + PreciseNumber magnitude = Abs(x); + if (magnitude > One) + { + throw new ArgumentOutOfRangeException(nameof(x), x, InverseDomainMessage); + } + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + TrigonometricGuardDigits; + + if (magnitude == One) + { + PreciseNumber halfPi = Divide(PiTo(working), Two, significantDigits); + return x.Significand.Sign > 0 ? halfPi : -halfPi; + } + + PreciseNumber denominator = Sqrt(Subtract(One, Multiply(x, x)), working); + return Atan(Divide(x, denominator, working), significantDigits); + } + + /// + /// Returns the arc cosine of a value, in radians. + /// + /// The value, which must lie in [-1, 1]. + /// The angle in radians whose cosine is , in [0, π]. + /// Thrown when lies outside [-1, 1]. + public static PreciseNumber Acos(PreciseNumber x) => + Acos(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the arc cosine of a value, in radians, to a chosen number of significant digits. + /// + /// The value, which must lie in [-1, 1]. + /// The number of significant digits to produce. + /// The angle in radians whose cosine is , in [0, π]. + /// + /// Thrown when lies outside [-1, 1], or when + /// is less than one. + /// + /// + /// acos x = π/2 - asin x, with the endpoints returned exactly: acos(1) = 0 and + /// acos(-1) = π. + /// + public static PreciseNumber Acos(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + PreciseNumber magnitude = Abs(x); + if (magnitude > One) + { + throw new ArgumentOutOfRangeException(nameof(x), x, InverseDomainMessage); + } + + int working = significantDigits + TrigonometricGuardDigits; + + if (x == One) + { + return Zero; + } + + if (x == NegativeOne) + { + return PiTo(working).ReduceSignificance(significantDigits); + } + + PreciseNumber halfPi = Divide(PiTo(working), Two, working); + return Subtract(halfPi, Asin(x, working)).ReduceSignificance(significantDigits); + } + + /// + /// Returns the arc tangent of a value, in radians. + /// + /// The value. + /// The angle in radians whose tangent is , in (-π/2, π/2). + public static PreciseNumber Atan(PreciseNumber x) => + Atan(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the arc tangent of a value, in radians, to a chosen number of significant digits. + /// + /// The value. + /// The number of significant digits to produce. + /// The angle in radians whose tangent is , in (-π/2, π/2). + /// Thrown when is less than one. + /// + /// atan x = 2 · atan( x / (1 + √(1 + x²)) ), applied until the argument is small, then the + /// Taylor series. Each application roughly halves the argument, so any magnitude is brought into + /// range in a bounded number of steps and the result is scaled back by the matching power of two. + /// + public static PreciseNumber Atan(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + TrigonometricGuardDigits; + + int reductions = 0; + PreciseNumber argument = x; + while (Abs(argument) > AtanReductionLimit && reductions < MaximumAtanReductions) + { + PreciseNumber root = Sqrt(Add(One, Multiply(argument, argument)), working + reductions); + argument = Divide(argument, Add(One, root), working + reductions); + reductions++; + } + + PreciseNumber series = AtanSeries(argument, working + reductions); + PreciseNumber result = reductions == 0 + ? series + : Multiply(new(0, BigInteger.One << reductions), series); + + return result.ReduceSignificance(significantDigits); + } + + /// + /// Returns the angle in radians whose tangent is y / x, using the signs of both to place + /// the angle in the correct quadrant. + /// + /// The ordinate. + /// The abscissa. + /// The angle in radians, in (-π, π]. + /// + /// Not part of : that lives on + /// , which does not + /// implement because it has no NaN or infinity to give the interface's edge cases meaning. It is + /// offered here as a bespoke static because a conversion from Cartesian to polar coordinates needs + /// it. Produced to the greater of the two arguments' significant digits, and never fewer than + /// . + /// + public static PreciseNumber Atan2(PreciseNumber y, PreciseNumber x) => + Atan2(y, x, Math.Max(DefaultTrigonometricPrecision(y), DefaultTrigonometricPrecision(x))); + + /// + /// Returns the angle in radians whose tangent is y / x, to a chosen number of significant + /// digits, using the signs of both to place the angle in the correct quadrant. + /// + /// The ordinate. + /// The abscissa. + /// The number of significant digits to produce. + /// The angle in radians, in (-π, π]. + /// Thrown when is less than one. + /// + /// Quadrant dispatch on atan(y / x), with the axes handled explicitly: on the positive + /// abscissa the arc tangent stands alone; on the negative abscissa it is offset by ±π to + /// carry the angle into the correct half; and where the abscissa is zero the angle is ±π/2, + /// or zero when the ordinate is zero as well. + /// + public static PreciseNumber Atan2(PreciseNumber y, PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + int working = significantDigits + TrigonometricGuardDigits; + + if (x.Significand.IsZero) + { + if (y.Significand.IsZero) + { + return Zero; + } + + PreciseNumber halfPi = Divide(PiTo(working), Two, significantDigits); + return y.Significand.Sign > 0 ? halfPi : -halfPi; + } + + PreciseNumber baseAngle = Atan(Divide(y, x, working), working); + + if (x.Significand.Sign > 0) + { + return baseAngle.ReduceSignificance(significantDigits); + } + + PreciseNumber pi = PiTo(working); + PreciseNumber shifted = y.Significand.Sign < 0 + ? Subtract(baseAngle, pi) + : Add(baseAngle, pi); + return shifted.ReduceSignificance(significantDigits); + } + + /// + /// Returns the sine of a value given in half turns. + /// + /// The argument, in half turns, so that SinPi(x) = Sin(x · π). + /// The sine of x · π. + /// + /// The argument is reduced modulo two before the multiplication by π, so a large argument keeps + /// its meaning: SinPi(1e20) is well-defined where Sin(1e20 · π) is not. + /// + public static PreciseNumber SinPi(PreciseNumber x) => + SinPi(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the sine of a value given in half turns, to a chosen number of significant digits. + /// + /// The argument, in half turns, so that SinPi(x) = Sin(x · π). + /// The number of significant digits to produce. + /// The sine of x · π. + /// Thrown when is less than one. + public static PreciseNumber SinPi(PreciseNumber x, int significantDigits) => + SinCosPi(x, significantDigits).SinPi; + + /// + /// Returns the cosine of a value given in half turns. + /// + /// The argument, in half turns, so that CosPi(x) = Cos(x · π). + /// The cosine of x · π. + /// + /// The argument is reduced modulo two before the multiplication by π, so a large argument keeps + /// its meaning: CosPi(1e20) is well-defined where Cos(1e20 · π) is not. + /// + public static PreciseNumber CosPi(PreciseNumber x) => + CosPi(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the cosine of a value given in half turns, to a chosen number of significant digits. + /// + /// The argument, in half turns, so that CosPi(x) = Cos(x · π). + /// The number of significant digits to produce. + /// The cosine of x · π. + /// Thrown when is less than one. + public static PreciseNumber CosPi(PreciseNumber x, int significantDigits) => + SinCosPi(x, significantDigits).CosPi; + + /// + /// Returns the sine and cosine of a value given in half turns. + /// + /// The argument, in half turns. + /// A tuple of the sine and cosine of x · π. + public static (PreciseNumber SinPi, PreciseNumber CosPi) SinCosPi(PreciseNumber x) => + SinCosPi(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the sine and cosine of a value given in half turns, to a chosen number of significant + /// digits. + /// + /// The argument, in half turns. + /// The number of significant digits to produce. + /// A tuple of the sine and cosine of x · π. + /// Thrown when is less than one. + /// + /// x = q/2 + f with q the nearest integer number of quarter turns and f in + /// [-1/4, 1/4] half turns. Only f · π reaches a series, and it is small however + /// large x is, so π is needed only to the width of the answer — the argument's magnitude + /// went into the integer q, which the octant selection consumes exactly. + /// + public static (PreciseNumber SinPi, PreciseNumber CosPi) SinCosPi(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return (Zero, One); + } + + int working = significantDigits + TrigonometricGuardDigits; + + // q counts quarter turns; the subtraction of q/2 is exact, so the argument's magnitude never + // reaches the multiplication by π and no wide constant is needed. + BigInteger quarters = RoundToNearestInteger(Multiply(x, Two)); + PreciseNumber fraction = Subtract(x, Divide(new(0, quarters), Two, working)); + + if (fraction.Significand.IsZero) + { + return SelectOctant(quarters, Zero, One, significantDigits); + } + + PreciseNumber radians = Multiply(fraction, PiTo(working)).ReduceSignificance(working); + (PreciseNumber sinRadians, PreciseNumber cosRadians) = SmallAngleSinCos(radians, working); + return SelectOctant(quarters, sinRadians, cosRadians, significantDigits); + } + + /// + /// Returns the tangent of a value given in half turns. + /// + /// The argument, in half turns, so that TanPi(x) = Tan(x · π). + /// The tangent of x · π. + public static PreciseNumber TanPi(PreciseNumber x) => + TanPi(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the tangent of a value given in half turns, to a chosen number of significant digits. + /// + /// The argument, in half turns, so that TanPi(x) = Tan(x · π). + /// The number of significant digits to produce. + /// The tangent of x · π. + /// Thrown when is less than one. + /// Thrown when the cosine of x · π is zero. + public static PreciseNumber TanPi(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + int working = significantDigits + TrigonometricGuardDigits; + (PreciseNumber sin, PreciseNumber cos) = SinCosPi(x, working); + return Divide(sin, cos, significantDigits); + } + + /// + /// Returns the arc sine of a value, in half turns. + /// + /// The value, which must lie in [-1, 1]. + /// asin(x) / π, in [-1/2, 1/2]. + /// Thrown when lies outside [-1, 1]. + public static PreciseNumber AsinPi(PreciseNumber x) => + AsinPi(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the arc sine of a value, in half turns, to a chosen number of significant digits. + /// + /// The value, which must lie in [-1, 1]. + /// The number of significant digits to produce. + /// asin(x) / π, in [-1/2, 1/2]. + /// + /// Thrown when lies outside [-1, 1], or when + /// is less than one. + /// + public static PreciseNumber AsinPi(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x == One) + { + return Half; + } + + if (x == NegativeOne) + { + return -Half; + } + + int working = significantDigits + TrigonometricGuardDigits; + return Divide(Asin(x, working), PiTo(working), significantDigits); + } + + /// + /// Returns the arc cosine of a value, in half turns. + /// + /// The value, which must lie in [-1, 1]. + /// acos(x) / π, in [0, 1]. + /// Thrown when lies outside [-1, 1]. + public static PreciseNumber AcosPi(PreciseNumber x) => + AcosPi(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the arc cosine of a value, in half turns, to a chosen number of significant digits. + /// + /// The value, which must lie in [-1, 1]. + /// The number of significant digits to produce. + /// acos(x) / π, in [0, 1]. + /// + /// Thrown when lies outside [-1, 1], or when + /// is less than one. + /// + public static PreciseNumber AcosPi(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x == One) + { + return Zero; + } + + if (x == NegativeOne) + { + return One; + } + + int working = significantDigits + TrigonometricGuardDigits; + return Divide(Acos(x, working), PiTo(working), significantDigits); + } + + /// + /// Returns the arc tangent of a value, in half turns. + /// + /// The value. + /// atan(x) / π, in (-1/2, 1/2). + public static PreciseNumber AtanPi(PreciseNumber x) => + AtanPi(x, DefaultTrigonometricPrecision(x)); + + /// + /// Returns the arc tangent of a value, in half turns, to a chosen number of significant digits. + /// + /// The value. + /// The number of significant digits to produce. + /// atan(x) / π, in (-1/2, 1/2). + /// Thrown when is less than one. + public static PreciseNumber AtanPi(PreciseNumber x, int significantDigits) + { + RequireSignificantDigits(significantDigits); + + if (x.Significand.IsZero) + { + return Zero; + } + + int working = significantDigits + TrigonometricGuardDigits; + return Divide(Atan(x, working), PiTo(working), significantDigits); + } + + /// + /// Converts an angle in degrees to radians. + /// + /// The angle, in degrees. + /// The angle in radians. + /// + /// degrees · π / 180, with π read as the correctly-rounded literal + /// rather than derived from any other constant. Produced to the + /// significant digits of , and never fewer than + /// . + /// + public static PreciseNumber DegreesToRadians(PreciseNumber degrees) => + DegreesToRadians(degrees, DefaultTrigonometricPrecision(degrees)); + + /// + /// Converts an angle in degrees to radians, to a chosen number of significant digits. + /// + /// The angle, in degrees. + /// The number of significant digits to produce. + /// The angle in radians. + /// Thrown when is less than one. + public static PreciseNumber DegreesToRadians(PreciseNumber degrees, int significantDigits) + { + RequireSignificantDigits(significantDigits); + int working = significantDigits + TrigonometricGuardDigits; + return Divide(Multiply(degrees, PiTo(working)), OneEighty, significantDigits); + } + + /// + /// Converts an angle in radians to degrees. + /// + /// The angle, in radians. + /// The angle in degrees. + /// + /// radians · 180 / π, with π read as the correctly-rounded literal + /// rather than derived from any other constant. Produced to the + /// significant digits of , and never fewer than + /// . + /// + public static PreciseNumber RadiansToDegrees(PreciseNumber radians) => + RadiansToDegrees(radians, DefaultTrigonometricPrecision(radians)); + + /// + /// Converts an angle in radians to degrees, to a chosen number of significant digits. + /// + /// The angle, in radians. + /// The number of significant digits to produce. + /// The angle in degrees. + /// Thrown when is less than one. + public static PreciseNumber RadiansToDegrees(PreciseNumber radians, int significantDigits) + { + RequireSignificantDigits(significantDigits); + int working = significantDigits + TrigonometricGuardDigits; + return Divide(Multiply(radians, OneEighty), PiTo(working), significantDigits); + } + + /// + /// Gets the significant digits a trigonometric function produces when the caller does not choose. + /// + /// The value being operated on. + /// The significant digits of , or if that is more. + /// + /// The same rule , the roots and the + /// exponentials follow, so that a constant carrying digits does + /// not silently cap the expression at fifty. + /// + private static int DefaultTrigonometricPrecision(PreciseNumber value) => + Math.Max(value.SignificantDigits, MinimumDivisionPrecision); + + /// + /// Selects the sine and cosine of an angle from those of its reduced remainder and the quadrant + /// the reduction removed. + /// + /// The number of quarter turns the reduction subtracted. + /// The sine of the remainder, in [-π/4, π/4]. + /// The cosine of the remainder. + /// The number of significant digits to produce. + /// The sine and cosine of the original angle. + /// + /// Adding a quarter turn rotates (sin, cos) to (cos, -sin), so the quadrant modulo + /// four names one of four sign-and-swap patterns. + /// + private static (PreciseNumber Sin, PreciseNumber Cos) SelectOctant( + BigInteger quadrant, PreciseNumber sinRemainder, PreciseNumber cosRemainder, int significantDigits) + { + int octant = (int)(((quadrant % 4) + 4) % 4); + (PreciseNumber sin, PreciseNumber cos) = octant switch + { + 0 => (sinRemainder, cosRemainder), + 1 => (cosRemainder, -sinRemainder), + 2 => (-sinRemainder, -cosRemainder), + _ => (-cosRemainder, sinRemainder), + }; + + return (sin.ReduceSignificance(significantDigits), cos.ReduceSignificance(significantDigits)); + } + + /// + /// Computes the sine and cosine of an angle already reduced into [-π/4, π/4]. + /// + /// The reduced angle. + /// The significant digits to carry through the series. + /// The sine and cosine of . + /// + /// Halving the angle before the series and applying the double-angle identities afterwards trades + /// a handful of multiplications for most of the terms. Each doubling compounds whatever relative + /// error it is handed, so the series is carried one digit wider per halving. + /// + private static (PreciseNumber Sin, PreciseNumber Cos) SmallAngleSinCos(PreciseNumber angle, int workingDigits) + { + if (angle.Significand.IsZero) + { + return (Zero, One); + } + + int halvings = 0; + PreciseNumber reduced = angle; + while (halvings < MaximumTrigonometricHalvings && Abs(reduced) > SeriesArgumentLimit) + { + // Halving terminates, so this is exact whatever precision is asked of it. + reduced = Divide(reduced, Two, workingDigits + MaximumTrigonometricHalvings + TrigonometricGuardDigits); + halvings++; + } + + int series = workingDigits + halvings + TrigonometricGuardDigits; + (PreciseNumber sin, PreciseNumber cos) = SinCosSeries(reduced, series); + + for (int doubling = 0; doubling < halvings; doubling++) + { + PreciseNumber nextSin = Multiply(Two, Multiply(sin, cos)).ReduceSignificance(series); + PreciseNumber nextCos = Subtract(Multiply(cos, cos), Multiply(sin, sin)).ReduceSignificance(series); + sin = nextSin; + cos = nextCos; + } + + return (sin.ReduceSignificance(workingDigits), cos.ReduceSignificance(workingDigits)); + } + + /// + /// Sums the sine and cosine series for a small argument. + /// + /// The argument, whose magnitude is at or below . + /// The significant digits to carry through the sums. + /// The sine and cosine of . + /// Thrown when a series does not converge. + /// + /// sin x = Σ (-1)ⁿ x^(2n+1) / (2n+1)! and cos x = Σ (-1)ⁿ x^(2n) / (2n)!, summed + /// together so the shared power of x² is formed once per term. Both stop as soon as their + /// terms fall below the last digit being carried. + /// + private static (PreciseNumber Sin, PreciseNumber Cos) SinCosSeries(PreciseNumber x, int workingDigits) + { + PreciseNumber negativeXSquared = -Multiply(x, x).ReduceSignificance(workingDigits); + + PreciseNumber sinTerm = x; + PreciseNumber sinSum = x; + PreciseNumber cosTerm = One; + PreciseNumber cosSum = One; + + for (int k = 1; k <= SeriesIterationAllowance(workingDigits); k++) + { + // cos term k is the previous one times -x² / ((2k-1)(2k)); sin term k times -x² / ((2k)(2k+1)). + cosTerm = Divide(Multiply(cosTerm, negativeXSquared), new(0, (long)((2 * k) - 1) * (2 * k)), workingDigits) + .ReduceSignificance(workingDigits); + sinTerm = Divide(Multiply(sinTerm, negativeXSquared), new(0, (long)(2 * k) * ((2 * k) + 1)), workingDigits) + .ReduceSignificance(workingDigits); + + PreciseNumber nextCos = Add(cosSum, cosTerm).ReduceSignificance(workingDigits); + PreciseNumber nextSin = Add(sinSum, sinTerm).ReduceSignificance(workingDigits); + + bool settled = nextCos == cosSum && nextSin == sinSum; + bool exhausted = cosTerm.Significand.IsZero && sinTerm.Significand.IsZero; + + cosSum = nextCos; + sinSum = nextSin; + + if (settled || exhausted) + { + return (sinSum, cosSum); + } + } + + throw new ArithmeticException( + $"The trigonometric series did not converge to {workingDigits.ToString(InvariantCulture)} significant digits."); + } + + /// + /// Sums the arc tangent series for a small argument. + /// + /// The argument, whose magnitude is at or below . + /// The significant digits to carry through the sum. + /// atan z. + /// Thrown when the series does not converge. + /// + /// atan z = z - z³/3 + z⁵/5 - …. The terms alternate in sign, and the sum stops as soon as + /// one falls below the last digit being carried. + /// + private static PreciseNumber AtanSeries(PreciseNumber z, int workingDigits) + { + if (z.Significand.IsZero) + { + return Zero; + } + + PreciseNumber negativeZSquared = -Multiply(z, z).ReduceSignificance(workingDigits); + PreciseNumber term = z; + PreciseNumber sum = z; + + for (int k = 1; k <= SeriesIterationAllowance(workingDigits); k++) + { + term = Multiply(term, negativeZSquared).ReduceSignificance(workingDigits); + if (term.Significand.IsZero) + { + return sum; + } + + PreciseNumber next = Add(sum, Divide(term, new(0, (2 * k) + 1), workingDigits)) + .ReduceSignificance(workingDigits); + + if (next == sum) + { + return sum; + } + + sum = next; + } + + throw new ArithmeticException( + $"The arc tangent series did not converge to {workingDigits.ToString(InvariantCulture)} significant digits."); + } +} diff --git a/README.md b/README.md index 7601026..8d93a88 100644 --- a/README.md +++ b/README.md @@ -82,7 +82,7 @@ A high-precision numeric type for .NET that provides arbitrary precision arithme - **Value Type**: A `readonly record struct` whose `default` value is zero. Adding, subtracting, multiplying, and comparing allocate nothing when the operands and every intermediate and final significand fit in an `int`. Exponent alignment counts, so `1 + 0.0000000001` allocates because it scales 1 by 10^10, and `99999 * 99999` allocates because its product is 9,999,800,001. -- **Comprehensive Mathematical Support**: Includes advanced mathematical functions like exponential operations (Pow, Exp, Squared, Cubed), roots (Sqrt, Cbrt, RootN, Hypot) through `IRootFunctions`, constant values (Pi, E, Tau) with high precision, absolute value operations, and specialized numerical checks (isOdd, isEven, etc.)—all with arbitrary precision. +- **Comprehensive Mathematical Support**: Includes advanced mathematical functions like exponential operations (Pow, Exp, Squared, Cubed), roots (Sqrt, Cbrt, RootN, Hypot) through `IRootFunctions`, trigonometry (Sin, Cos, Tan, Asin, Acos, Atan, Atan2, and the half-turn family) through `ITrigonometricFunctions`, constant values (Pi, E, Tau) with high precision, absolute value operations, and specialized numerical checks (isOdd, isEven, etc.)—all with arbitrary precision. - **Balanced Performance**: The design prioritizes accuracy and precision while maintaining reasonable performance. For calculations where extreme precision matters more than raw speed, PreciseNumber delivers excellent results, though built-in numeric types remain faster for standard precision needs. @@ -458,9 +458,18 @@ The `…M1` and `…P1` variants — `ExpM1`, `LogP1` and their siblings — are than as `Exp(x) - 1` and `Log(1 + x)`, so they keep the digits of a small argument instead of cancelling them away: `ExpM1(1e-30)` is `1e-30`, not zero. +`Sin`, `Cos`, `Tan`, their inverses, and `Atan2` follow the same rule again, and the accuracy of +each is the accuracy of the `π` it reduces by: reducing an angle of magnitude `10^d` to `n` correct +digits reads `π` to roughly `d + n` digits, so a large angle stays meaningful — `Sin(1000000)` is +correct to far more digits than a `double` holds. `SinCos` does that reduction once for a caller +that needs both. The half-turn family — `SinPi`, `CosPi`, `TanPi` and the inverses — reduces on the +argument before multiplying by `π`, so `SinPi(1e20)` is well-defined where `Sin(1e20 · π)` is not. +`Atan2` is a bespoke static rather than an interface member, because `PreciseNumber` has no NaN or +infinity to give `IFloatingPointIeee754`'s edge cases meaning. + ## Limitations -- There is no NaN, so `Sqrt()` of a negative value, `RootN()` of a negative value at an even degree, `Log()` of a value that is not positive, and `Pow()` of a negative value with a fractional exponent, throw `ArgumentOutOfRangeException` where a `double` would return NaN and carry on +- There is no NaN, so `Sqrt()` of a negative value, `RootN()` of a negative value at an even degree, `Log()` of a value that is not positive, `Pow()` of a negative value with a fractional exponent, and `Asin()` or `Acos()` of a value outside `[-1, 1]`, throw `ArgumentOutOfRangeException` where a `double` would return NaN and carry on. `Tan()` and `TanPi()` throw `DivideByZeroException` where the cosine is exactly zero - There is no infinity, so an exponential whose result needs a decimal exponent outside the range of an `int` throws `OverflowException` rather than saturating @@ -484,6 +493,8 @@ cancelling them away: `ExpM1(1e-30)` is `1e-30`, not zero. - **Roots**: `Sqrt()`, `Cbrt()`, `RootN()`, `Hypot()`, each with an overload taking the significant digits to produce +- **Trigonometry**: `Sin()`, `Cos()`, `SinCos()`, `Tan()`, `Asin()`, `Acos()`, `Atan()`, `Atan2()`, the half-turn family (`SinPi()`, `CosPi()`, `SinCosPi()`, `TanPi()`, `AsinPi()`, `AcosPi()`, `AtanPi()`), and `DegreesToRadians()` / `RadiansToDegrees()`, each with an overload taking the significant digits to produce + - **Utility**: `ToString()`, `Parse()`, `TryParse()`, `To()` - **Generic Conversion**: `TryConvertFromChecked`, `TryConvertFromSaturating`, `TryConvertFromTruncating`, `TryConvertToChecked`, `TryConvertToSaturating`, and `TryConvertToTruncating`, reached through `CreateChecked`, `CreateSaturating`, and `CreateTruncating`