Skip to content

Instantly share code, notes, and snippets.

@mellinoe

mellinoe/AcosTests.cs

Created March 7, 2017 01:38
Show Gist options
  • Save mellinoe/86dafbbcdf4ba6a5e8bd46d727adbcb4 to your computer and use it in GitHub Desktop.
Save mellinoe/86dafbbcdf4ba6a5e8bd46d727adbcb4 to your computer and use it in GitHub Desktop.
Download ZIP
Simple test program showing differences in actual and expected behavior in the .NET implementation of Complex.Acos.
Raw
AcosTests.cs
using System;
using System.Numerics;
using static System.Numerics.Complex;
namespace ConsoleApp
{
public class Program
{
// This is the largest x for which 2 x^2 will not overflow. It is used for branching inside Asin and Acos.
private static readonly double s_asinOverflowThreshold = Math.Sqrt(double.MaxValue) / 2.0;
static void Main(string[] args)
{
ShowDiff(0, 0, Math.PI / 2, 0);
ShowDiff(0, double.NaN, Math.PI / 2, double.NaN);
ShowDiff(10, double.PositiveInfinity, Math.PI / 2, double.NegativeInfinity);
ShowDiff(10, Double.NaN, double.NaN, double.NaN);
ShowDiff(double.NegativeInfinity, 10, Math.PI, double.NegativeInfinity);
ShowDiff(double.PositiveInfinity, 10, 0, double.NegativeInfinity);
ShowDiff(double.NegativeInfinity, double.PositiveInfinity, 3 * Math.PI / 4, double.NegativeInfinity);
ShowDiff(double.PositiveInfinity, double.PositiveInfinity, Math.PI / 4, double.NegativeInfinity);
ShowDiff(double.PositiveInfinity, double.NaN, double.NaN, double.PositiveInfinity);
ShowDiff(double.NaN, 10, double.NaN, double.NaN);
ShowDiff(double.NaN, double.PositiveInfinity, double.NaN, double.NegativeInfinity);
ShowDiff(double.NaN, double.NaN, double.NaN, double.NaN);
}
public static void ShowDiff(double r, double i, double rE, double iE)
{
ShowDiff(new Complex(r, i), new Complex(rE, iE));
}
public static void ShowDiff(Complex c, Complex expected)
{
Console.WriteLine($"Acos( {c.Real}, {c.Imaginary} )");
Console.WriteLine($" OLD: {Acos(c)}");
Console.WriteLine($" NEW: {Acos_New(c)}");
Console.WriteLine($" EXP: {expected}");
}
public static Complex Asin_Old(Complex value)
{
if ((value.Imaginary == 0 && value.Real < 0) || value.Imaginary > 0)
{
return -Asin_Old(-value);
}
return (-ImaginaryOne) * Log(ImaginaryOne * value + Sqrt(One - value * value));
}
public static Complex Asin_New(Complex value)
{
double b, bPrime, v;
Asin_Internal(Math.Abs(value.Real), Math.Abs(value.Imaginary), out b, out bPrime, out v);
double u;
if (bPrime < 0.0)
{
u = Math.Asin(b);
}
else
{
u = Math.Atan(bPrime);
}
if (value.Real < 0.0) u = -u;
if (value.Imaginary < 0.0) v = -v;
return new Complex(u, v);
}
public static Complex Acos_New(Complex value)
{
double b, bPrime, v;
Asin_Internal(Math.Abs(value.Real), Math.Abs(value.Imaginary), out b, out bPrime, out v);
double u;
if (bPrime < 0.0)
{
u = Math.Acos(b);
}
else
{
u = Math.Atan(1.0 / bPrime);
}
if (value.Real < 0.0) u = Math.PI - u;
if (value.Imaginary > 0.0) v = -v;
return new Complex(u, v);
}
private static void Asin_Internal(double x, double y, out double b, out double bPrime, out double v)
{
// This method for the inverse complex sine (and cosine) is described in
// Hull and Tang, "Implementing the Complex Arcsine and Arccosine Functions Using Exception Handling",
// ACM Transactions on Mathematical Software (1997)
// (https://www.researchgate.net/profile/Ping_Tang3/publication/220493330_Implementing_the_Complex_Arcsine_and_Arccosine_Functions_Using_Exception_Handling/links/55b244b208ae9289a085245d.pdf)
// First, the basics: start with sin(w) = (e^{iw} - e^{-iw})/(2i) = z. Here z is the input
// and w is the output. To solve for w, define t = e^{i w} and multiply through by t to
// get the quadratic equation t^2 - 2 i z t - 1 = 0. The solution is t = i z + sqrt(1 - z^2), so
// w = arcsin(z) = - i log( i z + sqrt(1 - z^2) )
// Decompose z = x + i y, multiply out i z + sqrt(1 - z^2), use log(s) = |s| + i arg(s), and do a
// bunch of algebra to get the components of w = arcsin(z) = u + i v
// u = arcsin(\beta) v = sign(y) log(\alpha + sqrt(\alpha^2 - 1))
// where
// \alpha = (\rho + \sigma)/2 \beta = (\rho - \sigma)/2
// \rho = sqrt((x + 1)^2 + y^2) \sigma = \sqrt((x - 1)^2 + y^2)
// This appears in DLMF section 4.23. (http://dlmf.nist.gov/4.23), along with the analogous
// arccos(w) = arccos(\beta) - i sign(y) log(\alpha + sqrt(\alpha^2 - 1))
// So \alpha and \beta together give us arcsin(w) and arccos(w).
// As written, \alpha is not susceptible to cancelation errors, but \beta is. To avoid cancelation, note
// \beta = (\rho^2 - \sigma^2)/(\rho + \sigma))/2 = (2 x)/(\rho + \sigma) = x / \alpha
// which is not subject to cancelation. Note \alpha >= 1 and |\beta| <= 1.
// For \alpha ~ 1, the argument of the log is near unity, so we compute (\alpha - 1) instead,
// and write the argument as 1 + (\alpha - 1) + \sqrt{(\alpha - 1)(\alpha + 1)}.
// For \beta ~ 1, arccos does not accurately resolve small angles, so we compute the tangent of the angle
// instead. Hull and Tang derive formulas for (\alpha - 1) and \beta' = \tan(u) that do not suffer
// cancelation for these cases.
// For simplicity, we assume all positive inputs and return all positive outputs. The caller should
// assign signs appropriate to the desired cut conventions. We return v directly since its magnitude
// is the same for both arcsin and arccos. Instead of u, we usually return \beta and sometimes \beta'.
// If \beta' is not computed, it is set to -1; if it is computed, it should be used instead of \beta
// to determine u. Compute u = \arcsin(\beta) or u = \arctan(\beta') for arcsin, u = \arccos(\beta)
// or \arctan(1/\beta') for arccos.
// For x or y large enough to overflow \alpha^2, we can simplify our formulas and avoid overflow.
if ((x > s_asinOverflowThreshold) || (y > s_asinOverflowThreshold))
{
b = -1.0;
bPrime = x / y;
double small, big;
if (x < y)
{
small = x;
big = y;
}
else
{
small = y;
big = x;
}
double ratio = small / big;
v = Math.Log(2.0) + Math.Log(big) + 0.5 * Log1P(ratio * ratio);
}
else
{
double r = Hypot((x + 1.0), y);
double s = Hypot((x - 1.0), y);
double a = (r + s) / 2.0;
b = x / a;
if (b > 0.75)
{
double amx;
if (x <= 1.0)
{
amx = (y * y / (r + (x + 1.0)) + (s + (1.0 - x))) / 2.0;
}
else
{
amx = y * y * (1.0 / (r + (x + 1.0)) + 1.0 / (s + (x - 1.0))) / 2.0;
}
bPrime = x / Math.Sqrt((a + x) * amx);
}
else
{
bPrime = -1.0;
}
if (a < 1.5)
{
double am1;
if (x < 1.0)
{
am1 = y * y / 2.0 * (1.0 / (r + (x + 1.0)) + 1.0 / (s + (1.0 - x)));
}
else
{
am1 = (y * y / (r + (x + 1.0)) + (s + (x - 1.0))) / 2.0;
}
v = Log1P(am1 + Math.Sqrt(am1 * (a + 1.0)));
}
else
{
// Because of the test above, we can be sure that a * a will not overflow.
v = Math.Log(a + Math.Sqrt((a - 1.0) * (a + 1.0)));
}
}
}
private static double Hypot(double a, double b)
{
// Using
// sqrt(a^2 + b^2) = |a| * sqrt(1 + (b/a)^2)
// we can factor out the larger component to dodge overflow even when a * a would overflow.
a = Math.Abs(a);
b = Math.Abs(b);
double small, large;
if (a < b)
{
small = a;
large = b;
}
else
{
small = b;
large = a;
}
if (small == 0.0)
{
return (large);
}
else if (double.IsPositiveInfinity(large) && !double.IsNaN(small))
{
// The NaN test is necessary so we don't return +inf when small=NaN and large=+inf.
// NaN in any other place returns NaN without any special handling.
return (double.PositiveInfinity);
}
else
{
double ratio = small / large;
return (large * Math.Sqrt(1.0 + ratio * ratio));
}
}
private static double Log1P(double x)
{
// Compute log(1 + x) without loss of accuracy when x is small.
// Our only use case so far is for positive values, so this isn't coded to handle negative values.
double xp1 = 1.0 + x;
if (xp1 == 1.0)
{
return x;
}
else if (x < 0.75)
{
// This is accurate to within 5 ulp with any floating-point system that uses a guard digit,
// as proven in Theorem 4 of "What Every Computer Scientist Should Know About Floating-Point
// Arithmetic" (https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html)
return x * Math.Log(xp1) / (xp1 - 1.0);
}
else
{
return Math.Log(xp1);
}
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment