diff --git a/src/Numerics/Trigonometry.cs b/src/Numerics/Trigonometry.cs index 70c9534e..7fe39db9 100644 --- a/src/Numerics/Trigonometry.cs +++ b/src/Numerics/Trigonometry.cs @@ -192,11 +192,39 @@ namespace MathNet.Numerics return new Complex(Tan(value.Real), 0.0); } - var cosr = Cos(value.Real); - var sinhi = Sinh(value.Imaginary); - var denom = (cosr * cosr) + (sinhi * sinhi); + // tan(x + j*y) = (tan(x) + j*cosh(y)*sinh(y)/cos^2(x))/(1 + sinh^2(y)/cos^2(x)) + // if |y| > asinh(sqrt(max)), tan(z) = 4*cos(x)*sin(x)*exp(-2*|y|) + j*sign(y) + // if exp(-|y|) = 0, tan(z) = j*sign(y) + // if tan(x) = +/- oo or 1/cos^2(x) = 1 + tan^2(x) = oo, tan(z) = j*cosh(y)/sinh(y) + // + // The algorithm is from: + // Kahan, W. Branch Cuts for Complex Elementary Functions, or Much Ado + // About Nothing's Sign Bit." In The State of the Art in Numerical Analysis: + // Proceedings of the Joint IMA / SIAM Conference on the State of the Art in + // Numerical Analysis Held at the UN (Ed.A.Iserles and M.J.D.Powell). + // New York: Clarendon Press, pp. 165 - 211, 1987. + + double upperLimit = Asinh(Math.Sqrt(double.MaxValue)); // 355.58450362725193 + + if (Math.Abs(value.Imaginary) > upperLimit) + { + double e = Math.Exp(-Math.Abs(value.Imaginary)); + return e == 0.0 + ? new Complex(0.0, Math.Sign(value.Imaginary)) + : new Complex(4.0 * Math.Cos(value.Real) * Math.Sin(value.Real) * e * e, Math.Sign(value.Imaginary)); + } + + double tanr = Math.Tan(value.Real); + double beta = 1 + tanr * tanr; // beta = 1/cos(x)^2 = 1 + t^2 + double sinhi = Math.Sinh(value.Imaginary); + double coshi = Math.Cosh(value.Imaginary); + + if (double.IsInfinity(tanr)) + return new Complex(0.0, coshi / sinhi); - return new Complex(Sin(value.Real) * cosr / denom, sinhi * Cosh(value.Imaginary) / denom); + double denom = 1.0 + beta * sinhi * sinhi; + + return new Complex(tanr / denom, beta * coshi * sinhi / denom); } /// @@ -221,11 +249,34 @@ namespace MathNet.Numerics return new Complex(Cot(value.Real), 0d); } - var sinr = Sin(value.Real); - var sinhi = Sinh(value.Imaginary); - var denom = (sinr * sinr) + (sinhi * sinhi); + // cot(x + j*y) = (cot(x) - j*cosh(y)*sinh(y)/sin^2(x))/(1 + sinh^2(y)/sin^2(x)) + // if |y| > asinh(sqrt(max)), cot(z) = 4*cos(x)*sin(x)*exp(-2*|y|) - j*sign(y) + // if exp(-|y|) = 0, cot(z) = -j*sign(y) + // if cot(x) = +/- oo or 1/sin^2(x) = 1 + cot^2(x) = oo, cot(z) = -j*cosh(y)/sinh(y) + // + // The algorithm is based on Kahan. + + double upperLimit = Asinh(Math.Sqrt(double.MaxValue)); // 355.58450362725193 + + if (Math.Abs(value.Imaginary) > upperLimit) + { + double e = Math.Exp(-Math.Abs(value.Imaginary)); + return e == 0.0 + ? new Complex(0.0, -Math.Sign(value.Imaginary)) + : new Complex(4.0 * Math.Cos(value.Real) * Math.Sin(value.Real) * e * e, -Math.Sign(value.Imaginary)); + } + + double cotr = Cot(value.Real); + double beta = 1 + cotr * cotr; // beta = 1/cos(x)^2 = 1 + t^2 + double sinhi = Sinh(value.Imaginary); + double coshi = Cosh(value.Imaginary); + + if (double.IsInfinity(cotr)) + return new Complex(0.0, - coshi / sinhi); - return new Complex(sinr * Cos(value.Real) / denom, -sinhi * Cosh(value.Imaginary) / denom); + double denom = 1.0 + beta * sinhi * sinhi; + + return new Complex(cotr / denom, - beta * coshi * sinhi / denom); } /// @@ -515,17 +566,33 @@ namespace MathNet.Numerics return new Complex(Tanh(value.Real), 0.0); } - var cosi = Cos(value.Imaginary); - var sinhr = Sinh(value.Real); + // tanh(x + j*y) = (cosh(x)*sinh(x)/cos^2(y) + j*tan(y))/(1 + sinh^2(x)/cos^2(y)) + // if |x| > asinh(sqrt(max)), tanh(z) = sign(x) + j*4*cos(y)*sin(y)*exp(-2*|x|) + // if exp(-|x|) = 0, tanh(z) = sign(x) + // if tan(y) = +/- oo or 1/cos^2(y) = 1 + tan^2(y) = oo, tanh(z) = cosh(x)/sinh(x) + // + // The algorithm is based on Kahan. + + double upperLimit = Asinh(Math.Sqrt(double.MaxValue)); // 355.58450362725193 - if (double.IsInfinity(sinhr)) + if (Math.Abs(value.Real) > upperLimit) { - return new Complex(double.IsPositiveInfinity(sinhr) ? 1 : -1, 0.0); + double e = Math.Exp(-Math.Abs(value.Real)); + return e == 0.0 + ? new Complex(Math.Sign(value.Real), 0.0) + : new Complex(Math.Sign(value.Real), 4.0 * Math.Cos(value.Imaginary) * Math.Sin(value.Imaginary) * e * e); } - var denom = (cosi * cosi) + (sinhr * sinhr); + double tani = Tan(value.Imaginary); + double beta = 1 + tani * tani; // beta = 1/cos^2(y) = 1 + t^2 + double sinhr = Sinh(value.Real); + double coshr = Cosh(value.Real); + + if (double.IsInfinity(tani)) + return new Complex(coshr / sinhr, 0.0); - return new Complex(Cosh(value.Real) * sinhr / denom, cosi * Sin(value.Imaginary) / denom); + double denom = 1.0 + beta * sinhr * sinhr; + return new Complex(beta * coshr * sinhr / denom, tani / denom); } /// @@ -562,17 +629,33 @@ namespace MathNet.Numerics return new Complex(Coth(value.Real), 0.0); } - var sini = Sin(value.Imaginary); - var sinhr = Sinh(value.Real); + // coth(x + j*y) = (cosh(x)*sinh(x)/sin^2(y) - j*cot(y))/(1 + sinh^2(x)/sin^2(y)) + // if |x| > asinh(sqrt(max)), coth(z) = sign(x) - j*4*cos(y)*sin(y)*exp(-2*|x|) + // if exp(-|x|) = 0, coth(z) = sign(x) + // if cot(y) = +/- oo or 1/sin^2(y) = 1 + cot^2(y) = oo, coth(z) = cosh(x)/sinh(x) + // + // The algorithm is based on Kahan. - if (double.IsInfinity(sinhr)) + double upperLimit = Asinh(Math.Sqrt(double.MaxValue)); // 355.58450362725193 + + if (Math.Abs(value.Real) > upperLimit) { - return new Complex(double.IsPositiveInfinity(sinhr) ? 1 : -1, 0.0); + double e = Math.Exp(-Math.Abs(value.Real)); + return e == 0.0 + ? new Complex(Math.Sign(value.Real), 0.0) + : new Complex(Math.Sign(value.Real), -4.0 * Math.Cos(value.Imaginary) * Math.Sin(value.Imaginary) * e * e); } - var denom = (sini * sini) + (sinhr * sinhr); + double coti = Cot(value.Imaginary); + double beta = 1 + coti * coti; + double sinhr = Math.Sinh(value.Real); + double coshr = Math.Cosh(value.Real); + + if (double.IsInfinity(coti)) + return new Complex(coshr / sinhr, 0.0); - return new Complex(sinhr * Cosh(value.Real) / denom, -sini * Cos(value.Imaginary) / denom); + double denom = 1.0 + beta * sinhr * sinhr; + return new Complex(beta * coshr * sinhr / denom, -coti / denom); } /// @@ -597,14 +680,36 @@ namespace MathNet.Numerics return new Complex(Sech(value.Real), 0.0); } - var exp = value.Exp(); + // sech(x + j*y) = (cosh(x)/cos(y) - j*sinh(x)*tan(y)/cos(y))/(1 + sinh^2(x)/cos^2(y)) + // if |x| > asinh(sqrt(max)), sech(z) = 4*cosh(x)*cos(y)*exp(-2*|x|) - j*4*sinh(x)*tan(y)*cos(y)*exp(-2*|x|) + // if exp(-|x|) = 0, sech(z) = 0 + // if tan(y) = +/- oo or 1/cos^2(y) = 1 + tan^2(y) = oo, sech(z) = -j*sign(tan(y))/sinh(x) + // + // The algorithm is based on Kahan. + + double upperLimit = Asinh(Math.Sqrt(double.MaxValue)); // 355.58450362725193 + + double tani = Tan(value.Imaginary); + double cosi = Cos(value.Imaginary); + double beta = 1.0 + tani * tani; + double sinhr = Math.Sinh(value.Real); + double coshr = Math.Cosh(value.Real); + + if (Math.Abs(value.Real) > upperLimit) + { + double e = Math.Exp(-Math.Abs(value.Real)); + return e == 0.0 + ? new Complex(0, 0) + : new Complex(4.0 * coshr * cosi * e * e, -4.0 * sinhr * tani * cosi * e * e); + } - if (exp.IsInfinity()) + if (double.IsInfinity(tani)) { - return Complex.Zero; + return new Complex(0.0, -Math.Sign(tani) / sinhr); } - return 2 * exp / (exp.Square() + 1); + double denom = 1.0 + beta * sinhr * sinhr; + return new Complex(coshr / cosi / denom, -sinhr * tani / cosi / denom); } /// @@ -629,14 +734,36 @@ namespace MathNet.Numerics return new Complex(Csch(value.Real), 0.0); } - var exp = value.Exp(); + // csch(x + j*y) = (sinh(x)*cot(y)/sin(y) - j*cosh(x)/sin(y))/(1 + sinh^2(x)/sin^2(y)) + // if |x| > asinh(sqrt(max)), csch(z) = 4*sinh(x)*cot(y)*sin(y)*exp(-2*|x|) - j*4*cosh(x)*sin(y)*exp(-2*|x|) + // if exp(-|x|) = 0, csch(z) = 0 + // if cot(y) = +/- oo or 1/sin^2(x) = 1 + cot^2(x) = oo, csch(z) = sign(cot(y))/sinh(x) + // + // The algorithm is based on Kahan. + + double upperLimit = Asinh(Math.Sqrt(double.MaxValue)); + + double coti = Cot(value.Imaginary); + double sini = Sin(value.Imaginary); + double beta = 1 + coti * coti; + double sinhr = Sinh(value.Real); + double coshr = Cosh(value.Real); + + if (Math.Abs(value.Real) > upperLimit) + { + double e = Math.Exp(-Math.Abs(value.Real)); + return e == 0.0 + ? new Complex(0, 0) + : new Complex(4.0 * sinhr * coti * sini * e * e, -4.0 * coshr * sini * e * e); + } - if (exp.IsInfinity()) + if (double.IsInfinity(coti)) { - return Complex.Zero; + return new Complex(Math.Sign(coti) / sinhr, 0.0); } - return 2 * exp / (exp.Square() - 1); + double denom = 1.0 + beta * sinhr * sinhr; + return new Complex(sinhr * coti / sini / denom, -coshr / sini / denom); } diff --git a/src/UnitTests/TrigonometryTest.cs b/src/UnitTests/TrigonometryTest.cs index a428cda6..2fb3194b 100644 --- a/src/UnitTests/TrigonometryTest.cs +++ b/src/UnitTests/TrigonometryTest.cs @@ -58,6 +58,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(-1.19209289550780998537e-7, 0.0, 0.99999999999999289, 0.0)] [TestCase(8.388608e6, 1.19209289550780998537e-7, -0.90175467375876572, -5.1528001100635277e-8)] [TestCase(-1.19209289550780998537e-7, -8.388608e6, double.PositiveInfinity, double.NegativeInfinity)] + [TestCase(512.0, 32.0, -39356457675156.6, -3139507838262.74)] + [TestCase(-512.0, -32.0, -39356457675156.6, -3139507838262.74)] + [TestCase(32.0, -512.0, +9.52855589474963e+221, +6.29843301304523e+221)] + [TestCase(-32.0, -512.0, +9.52855589474963e+221, -6.29843301304523e+221)] public void CanComputeComplexCosine(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Cos(); @@ -79,6 +83,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(-1.19209289550780998537e-7, 0.0, -1.19209289550780998537e-7, 0.0)] [TestCase(8.388608e6, 1.19209289550780998537e-7, 0.43224820225680083, -1.0749753400787824e-7)] [TestCase(-1.19209289550780998537e-7, -8.388608e6, double.NegativeInfinity, double.NegativeInfinity)] + [TestCase(512.0, 32.0, +3139507838262.74, -39356457675156.6)] + [TestCase(-512.0, -32.0, -3139507838262.74, +39356457675156.6)] + [TestCase(32.0, -512.0, +6.29843301304523e+221, -9.52855589474963e+221)] + [TestCase(-32.0, -512.0, -6.29843301304523e+221, -9.52855589474963e+221)] public void CanComputeComplexSine(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Sin(); @@ -100,6 +108,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(-1.19209289550780998537e-7, 0.0, -1.1920928955078157e-7, 0.0)] [TestCase(8.388608e6, 1.19209289550780998537e-7, -0.47934123862653449, 1.4659977233982276e-7)] [TestCase(-8.388608e6, -1.19209289550780998537e-7, 0.47934123862653449, -1.4659977233982276e-7)] + [TestCase(512.0, 32.0, -5.08515122860094E-29, 1)] + [TestCase(-512.0, -32.0, +5.08515122860094E-29, -1)] + [TestCase(32.0, 512.0, +3.52598650243739e-445, 1)] + [TestCase(32.0, -512.0, +3.52598650243739e-445, -1)] public void CanComputeComplexTangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Tan(); @@ -549,6 +561,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(-1.19209289550780998537e-7, 0.0, -8388607.999999978, 0.0)] [TestCase(8.388608e6, 1.19209289550780998537e-7, -2.0861964701080704, -6.3803383253713457e-7)] [TestCase(-8.388608e6, -1.19209289550780998537e-7, 2.0861964701080704, 6.3803383253713457e-7)] + [TestCase(512.0, 32.0, -5.08515122860094E-29, -1.0)] + [TestCase(-512.0, -32.0, 5.08515122860094E-29, +1.0)] + [TestCase(32.0, 512.0, +3.52598650243739e-445, -1.0)] + [TestCase(32.0, -512.0, +3.52598650243739e-445, +1.0)] public void CanComputeComplexCotangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Cot(); @@ -570,6 +586,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(-1.19209289550780998537e-7, 0.0, 1.0000000000000071, 0.0)] [TestCase(8.388608e6, 1.19209289550780998537e-7, -1.1089490624226177, 6.3367488045143761e-8)] [TestCase(-8.388608e6, -1.19209289550780998537e-7, -1.1089490624226177, 6.3367488045143761e-8)] + [TestCase(512.0, 32.0, -2.52481261731330570134e-14, +2.01407074478744036413e-15)] + [TestCase(-512.0, -32.0, -2.52481261731330570134e-14, +2.01407074478744036413e-15)] + [TestCase(32.0, 512.0, +7.30361056703505e-223, +4.82773070945482e-223)] + [TestCase(32.0, -512.0, +7.30361056703505e-223, -4.82773070945482e-223)] public void CanComputeComplexSecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Sec(); @@ -591,6 +611,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(-1.19209289550780998537e-7, 0.0, -8388608.0000000376, 0.0)] [TestCase(8.388608e6, 1.19209289550780998537e-7, 2.3134856195557596, 5.7534999050657057e-7)] [TestCase(-8.388608e6, -1.19209289550780998537e-7, -2.3134856195557596, -5.7534999050657057e-7)] + [TestCase(512.0, 32.0, 2.01407074478744e-15, 2.52481261731331e-14)] + [TestCase(-512.0, -32.0, -2.01407074478744e-15, -2.52481261731331e-14)] + [TestCase(32.0, 512.0, +4.82773070945482e-223, -7.30361056703505e-223)] + [TestCase(32.0, -512.0, +4.82773070945482e-223, +7.30361056703505e-223)] public void CanComputeComplexCosecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Csc(); @@ -616,6 +640,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(0.5, -0.5, 0.45730415318424922, -0.54061268571315335)] [TestCase(-0.5, 0.5, -0.45730415318424922, 0.54061268571315335)] [TestCase(-0.5, -0.5, -0.45730415318424922, -0.54061268571315335)] + [TestCase(512.0, 32.0, 9.52855589474962752845e221, 6.29843301304522728211e221)] + [TestCase(-512.0, -32.0, -9.52855589474962752845e221, -6.29843301304522728211e221)] + [TestCase(-32.0, 512.0, +3.93564576751566e13, +3.13950783826274e12)] + [TestCase(32.0, -512.0, -3.93564576751566e13, -3.13950783826274e12)] public void CanComputeComplexHyperbolicSine(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Sinh(); @@ -641,6 +669,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(0.5, -0.5, 0.9895848833999199, -0.24982639750046154)] [TestCase(-0.5, 0.5, 0.9895848833999199, -0.24982639750046154)] [TestCase(-0.5, -0.5, 0.9895848833999199, 0.24982639750046154)] + [TestCase(512.0, 32.0, 9.52855589474962752845e221, 6.29843301304522728211e221)] + [TestCase(-512.0, -32.0, 9.52855589474962752845e221, 6.29843301304522728211e221)] + [TestCase(-32.0, 512.0, -3.93564576751566e13, -3.13950783826274e12)] + [TestCase(32.0, -512.0, -3.93564576751566e13, -3.13950783826274e12)] public void CanComputeComplexHyperbolicCosine(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Cosh(); @@ -666,6 +698,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(0.5, -0.5, 0.56408314126749848, -0.40389645531602575)] [TestCase(-0.5, 0.5, -0.56408314126749848, 0.40389645531602575)] [TestCase(-0.5, -0.5, -0.56408314126749848, -0.40389645531602575)] + [TestCase(512.0, 32.0, 1.0, 0.0)] + [TestCase(-512.0, -32.0, -1.0, 0.0)] + [TestCase(-32.0, 512.0, -1.0, -5.08515122860093626173e-29)] + [TestCase(32.0, -512.0, +1.0, +5.08515122860093626173e-29)] public void CanComputeComplexHyperbolicTangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Tanh(); @@ -691,6 +727,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(0.5, -0.5, 1.1719451445243514, 0.8391395790248311)] [TestCase(-0.5, 0.5, -1.1719451445243514, -0.8391395790248311)] [TestCase(-0.5, -0.5, -1.1719451445243514, 0.8391395790248311)] + [TestCase(512.0, 32.0, 1.0, 0.0)] + [TestCase(-512.0, -32.0, -1.0, 0.0)] + [TestCase(-32.0, 512.0, -1.0, 5.08515122860093626173e-29)] + [TestCase(32.0, -512.0, +1.0, -5.08515122860093626173e-29)] public void CanComputeComplexHyperbolicCotangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Coth(); @@ -716,6 +756,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(0.5, -0.5, 0.94997886761549463, 0.23982763093808804)] [TestCase(-0.5, 0.5, 0.94997886761549463, 0.23982763093808804)] [TestCase(-0.5, -0.5, 0.94997886761549463, -0.23982763093808804)] + [TestCase(512.0, 32.0, 7.30361056703505051822e-223, -4.82773070945482080706e-223)] + [TestCase(-512.0, -32.0, 7.30361056703505051822e-223, -4.82773070945482080706e-223)] + [TestCase(-32.0, 512.0, -2.52481261731330570134e-14, +2.01407074478744036413e-15)] + [TestCase(-32.0, -512.0, -2.52481261731330570134e-14, -2.01407074478744036413e-15)] public void CanComputeComplexHyperbolicSecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Sech(); @@ -741,6 +785,10 @@ namespace MathNet.Numerics.UnitTests [TestCase(0.5, -0.5, 0.91207426403881078, 1.0782296946540223)] [TestCase(-0.5, 0.5, -0.91207426403881078, -1.0782296946540223)] [TestCase(-0.5, -0.5, -0.91207426403881078, 1.0782296946540223)] + [TestCase(512.0, 32.0, 7.30361056703505051822e-223, -4.82773070945482080706e-223)] + [TestCase(-512.0, -32.0, -7.30361056703505051822e-223, +4.82773070945482080706e-223)] + [TestCase(-32.0, 512.0, +2.52481261731330570134e-14, -2.01407074478744036413e-15)] + [TestCase(32.0, -512.0, -2.52481261731330570134e-14, +2.01407074478744036413e-15)] public void CanComputeComplexHyperbolicCosecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Csch();