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`