ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

To avoid overflow, widen the range of x to apply the approximation of sinh^2(x). 

Simplify Coth(z), Tan(z), and Cot(z)
spatial
diluculo 9 years ago
parent
740fce12e0
commit
4aca8448bc
2 changed files with 42 additions and 91 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 105
      src/Numerics/Trigonometry.cs
  2. 28
      src/UnitTests/TrigonometryTest.cs

105
src/Numerics/Trigonometry.cs
View File

@ -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);
}
/// <summary>
@ -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);
}
/// <summary>
@ -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);
}
/// <summary>
@ -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);

28
src/UnitTests/TrigonometryTest.cs
View File

@ -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();

Write Preview
Loading…
Cancel
Save