diff --git a/src/Numerics/Trigonometry.cs b/src/Numerics/Trigonometry.cs index 7fe39db9..45b4b8fa 100644 --- a/src/Numerics/Trigonometry.cs +++ b/src/Numerics/Trigonometry.cs @@ -192,39 +192,10 @@ namespace MathNet.Numerics return new Complex(Tan(value.Real), 0.0); } - // 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); + // tan(z) = - j*tanh(j*z) - double denom = 1.0 + beta * sinhi * sinhi; - - return new Complex(tanr / denom, beta * coshi * sinhi / denom); + Complex z = Tanh(new Complex(-value.Imaginary, value.Real)); + return new Complex(z.Imaginary, -z.Real); } /// @@ -249,34 +220,10 @@ namespace MathNet.Numerics return new Complex(Cot(value.Real), 0d); } - // 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); - - double denom = 1.0 + beta * sinhi * sinhi; + // cot(z) = - j*coth(-j*z) - return new Complex(cotr / denom, - beta * coshi * sinhi / denom); + Complex z = Coth(new Complex(value.Imaginary, -value.Real)); + return new Complex(z.Imaginary, -z.Real); } /// @@ -567,13 +514,13 @@ namespace MathNet.Numerics } // 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 |x| > huge, 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 + double upperLimit = 177.6189650; // Asinh(double.MaxValue) / 4.0; if (Math.Abs(value.Real) > upperLimit) { @@ -629,33 +576,9 @@ namespace MathNet.Numerics return new Complex(Coth(value.Real), 0.0); } - // 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. - - double upperLimit = Asinh(Math.Sqrt(double.MaxValue)); // 355.58450362725193 - - if (Math.Abs(value.Real) > upperLimit) - { - 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); - } + // Coth(z) = 1/tanh(z) - 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); - - double denom = 1.0 + beta * sinhr * sinhr; - return new Complex(beta * coshr * sinhr / denom, -coti / denom); + return Complex.One / Tanh(value); } /// @@ -681,13 +604,13 @@ namespace MathNet.Numerics } // 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 |x| > huge, 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 upperLimit = 177.6189650; // Asinh(double.MaxValue) / 4.0; double tani = Tan(value.Imaginary); double cosi = Cos(value.Imaginary); @@ -735,13 +658,13 @@ namespace MathNet.Numerics } // 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 |x| > huge, 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 upperLimit = 177.6189650; // Asinh(double.MaxValue) / 4.0; double coti = Cot(value.Imaginary); double sini = Sin(value.Imaginary); diff --git a/src/UnitTests/TrigonometryTest.cs b/src/UnitTests/TrigonometryTest.cs index 2fb3194b..a7c4bc82 100644 --- a/src/UnitTests/TrigonometryTest.cs +++ b/src/UnitTests/TrigonometryTest.cs @@ -62,6 +62,8 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, -7.4732769989672537746e153, +2.25277102712558920722e149)] + [TestCase(355.58, 355.58, -1.11645148443766450218e154, +7.31511888973169437082e153)] public void CanComputeComplexCosine(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Cos(); @@ -87,6 +89,8 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, -2.25277102712558920722e149, -7.4732769989672537746e153)] + [TestCase(355.58, 355.58, -+7.31511888973169437082e153, -1.11645148443766450218e154)] public void CanComputeComplexSine(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Sin(); @@ -112,6 +116,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, +5.39739014654622322424e-313, +1.0)] + [TestCase(355.58, 355.58, +2.57308259096360322163e-309, +1.0)] + [TestCase(196.0, 16.0, +1.62991625102252204721e-14, +1.00000000000001938715)] public void CanComputeComplexTangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Tan(); @@ -565,6 +572,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, +5.39739014654622322424e-313, -1.0)] + [TestCase(355.58, 355.58, +2.57308259096360322163e-309, -1.0)] + [TestCase(196.0, 16.0, +1.62991625102252204721e-14, -1.00000000000001938715)] public void CanComputeComplexCotangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Cot(); @@ -590,6 +600,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, -1.33810107564527561878e-154, -4.03361916732891580715e-159)] + [TestCase(355.58, 355.58, -6.26665948010988663559e-155, -4.10598756663013051681e-155)] + [TestCase(196.0, 16.0, +7.70790453157574313253e-8, +2.11460357915116320189e-7)] public void CanComputeComplexSecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Sec(); @@ -615,6 +628,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, -4.03361916732891580715e-159, +1.33810107564527561878e-154)] + [TestCase(355.58, 355.58, -4.10598756663013051681e-155, +6.26665948010988663559e-155)] + [TestCase(196.0, 16.0, +2.11460357915116320189e-7, -7.70790453157574313253e-8)] public void CanComputeComplexCosecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Csc(); @@ -702,6 +718,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, +1.0, +5.39739014654622322424e-313)] + [TestCase(355.58, 355.58, +1.0, +2.57308259096360322163e-309)] + [TestCase(196.0, 16.0, +1.0, +6.29623407730470174659e-171)] public void CanComputeComplexHyperbolicTangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Tanh(); @@ -731,6 +750,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, +1.0, -5.39739014654622322424e-313)] + [TestCase(355.58, 355.58, +1.0, -2.57308259096360322163e-309)] + [TestCase(196.0, 16.0, +1.0, -6.29623407730470174659e-171)] public void CanComputeComplexHyperbolicCotangent(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Coth(); @@ -760,6 +782,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, -1.33810107564527561878e-154, +4.03361916732891580715e-159)] + [TestCase(355.58, 355.58, -6.26665948010988663559e-155, +4.10598756663013051681e-155)] + [TestCase(196.0, 16.0, -1.447180343166980306269e-85, +4.350690711792111509525e-86)] public void CanComputeComplexHyperbolicSecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Sech(); @@ -789,6 +814,9 @@ namespace MathNet.Numerics.UnitTests [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)] + [TestCase(355.0, 355.0, -1.33810107564527561878e-154, +4.03361916732891580715e-159)] + [TestCase(355.58, 355.58, -6.26665948010988663559e-155, +4.10598756663013051681e-155)] + [TestCase(196.0, 16.0, -1.447180343166980306269e-85, +4.350690711792111509525e-86)] public void CanComputeComplexHyperbolicCosecant(double real, double imag, double expectedReal, double expectedImag) { var actual = new Complex(real, imag).Csch();