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

Trig: updated Tan(z), Cot(z), Tanh(z), Coth(z), Csch(z), and Sech(z) based on the Kahan algorithm.

pull/540/head
diluculo 9 years ago
parent
11f3c54309
commit
740fce12e0
2 changed files with 203 additions and 28 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. 183
      src/Numerics/Trigonometry.cs
  2. 48
      src/UnitTests/TrigonometryTest.cs

183
src/Numerics/Trigonometry.cs
View File

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

48
src/UnitTests/TrigonometryTest.cs
View File

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

Write Preview
Loading…
Cancel
Save