How to compute the angles of a triangle accurately given squared edge lengths?

Question

Given the squared edge-lengths of a triangle as a, b, c > 0, how can we accurately compute the interior angles?

Kahan's oft cited method which carefully rearranges parentheses assumes (unsquared) edge lengths.

Another popular method (pp 15) assumes edge vectors.

My best bet so far is to use the Law of Cosines:

float angle = acos( (b+a-c) / (2.*sqrt(b*a)) );

but for some non-degenerate inputs this will return 0.0.

Admittedly, the way I'm generating valid input is by computing squared lengths of edge-vectors, but for the purpose of this question I'd like to ignore access to the original edge vectors. Here's a tricky test case:

// Triangle (0,0) (1,1) (1+e,1)
float e = 1e-7;
float a = ((1+e)-1)*((1+e)-1);
float b = (1+e)*(1+e) + 1*1;
float c = 1*1 + 1*1;

The ground-truth angles are approximately:

4.9999997529193436e-08 2.3561944901923448 0.7853981133974508

If I immediately take sqrt of a,b,c and apply Kahan's method I get:

0            3.14159  0

If I apply Law of Cosines, I get:

0            2.35619  0.785398

which is better, but I don't like the exact ==0 where the correct value is >0.

Further, if I scale a,b,c by a non-power of two, then Law of Cosines gives a very wrong result:

0            2.02856  1.11303

Is there a numerically more accurate way to compute angles directly from squared edge lengths?

Or...

Am I simply asking too much of floats here?

Are my squared edge-lengths in float violating some kind of triangle inequality?

My working implementations, numbers above using clang++ main.cpp -g && ./a.out on mac os.

Answer 1

After experimenting for many hours re-arranging floating-point expressions, exploiting fused multiply-add (FMA) and compensated sums, I concluded that computing via the Law of Cosines is not going to work robustly, not even when switching to double-float arithmetic (also known as pair arithmetic).

As already noted in the question, Kahan's numerically superior arrangement by itself is insufficient since the need to take square roots already injects numerical error before the angle computation starts. However, I observed that performing intermediate computation in double arithmetic lead to a robust implementation. Since the asker's requirements preclude the use of double computation, this leaves us with double-float computation as a back-up, which of course has a significant performance impact even on platforms with FMA support. A "smoke" test indicates that this leads to an implementation capable of providing all angles of a triangle with relative error of less than 2^-23 by straightforward translation of Kahan's algorithm specification.

For the C++11 code below I assumed that the target platform has support for FMA which is thus available to accelerate the double-float functions. My test framework is based on a very old arbitrary-precision library that I have been using for the past 30 years: R. P. Brent's MP from 1978. I configure it for 250-bit precision. The reference function using the MP library computes the angles with the Law of Cosines to provide a robust unit test. This portion of the code would have to be replaced to utilize commonly-used modern libraries of this kind. I built with the Intel C/C++ compiler with full optimization and maximum conformance to IEEE-754 floating-point requirements (/fp:strict).

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

// data structures and functions for double-float computation

typedef struct float2 {
    float x;
    float y;
} float2;

float2 make_float2 (float head, float tail);
float2 add_dblflt (float2 a, float2 b);
float2 sub_dblflt (float2 a, float2 b);
float2 mul_dblflt (float2 a, float2 b);
float2 div_dblflt (float2 a, float2 b);
float2 sqrt_dblflt (float2 a);

// Compute angle C of triangle with squared edges asq, bsq, csq
float angle_kahanf (float asq, float bsq, float csq)
{
    if (asq < bsq) { float t = bsq; bsq = asq; asq = t; } // ensure asq >= bsq
    float2 a = sqrt_dblflt (make_float2 (asq, 0));
    float2 b = sqrt_dblflt (make_float2 (bsq, 0));
    float2 c = sqrt_dblflt (make_float2 (csq, 0));
    float2 mu = {INFINITY / INFINITY, INFINITY / INFINITY};
    if      ((bsq >= csq) && (csq >= 0)) mu = sub_dblflt (c, sub_dblflt (a, b));
    else if ((csq >    0) && (bsq >= 0)) mu = sub_dblflt (b, sub_dblflt (a, c));
    else    fprintf (stderr, "angle_kahanf: not a real triangle\n");
    float2 fact_0 = add_dblflt (sub_dblflt (a, b), c);
    float2 num = mul_dblflt (fact_0, mu);
    float2 fact_1 = add_dblflt (a, add_dblflt (b, c));
    float2 fact_2 = add_dblflt (b, sub_dblflt (a, c));
    float2 den = mul_dblflt (fact_1, fact_2);
    float2 ratio = div_dblflt (num, den);
    float2 root = sqrt_dblflt (ratio);
    float atan_val = atanf (root.y);
    float angle = 2.0f * atan_val;
    return angle;
}

float2 make_float2 (float head, float tail)
{
    float2 r;
    r.x = tail;  // least signficant
    r.y = head;  // most signficant
    return r;
}

float2 add_dblflt (float2 a, float2 b)
{
    float2 z;
    float t1, t2, t3, t4, t5, e;
    t1 = a.y + b.y;
    t2 = t1 - a.y;
    t3 = (a.y + (t2 - t1)) + (b.y - t2);
    t4 = a.x + b.x;
    t2 = t4 - a.x;
    t5 = (a.x + (t2 - t4)) + (b.x - t2);
    t3 = t3 + t4;
    t4 = t1 + t3;
    t3 = (t1 - t4) + t3;
    t3 = t3 + t5;
    z.y = e = t4 + t3;
    z.x = (t4 - e) + t3;
    return z;
}

float2 sub_dblflt (float2 a, float2 b)
{
    float2 z;
    float t1, t2, t3, t4, t5, e;
    t1 = a.y - b.y;
    t2 = t1 - a.y;
    t3 = (a.y + (t2 - t1)) - (b.y + t2);
    t4 = a.x - b.x;
    t2 = t4 - a.x;
    t5 = (a.x + (t2 - t4)) - (b.x + t2);
    t3 = t3 + t4;
    t4 = t1 + t3;
    t3 = (t1 - t4) + t3;
    t3 = t3 + t5;
    z.y = e = t4 + t3;
    z.x = (t4 - e) + t3;
    return z;
}
 
float2 mul_dblflt (float2 a, float2 b)

{
    float2 t, z;
    float e;
    t.y = a.y * b.y;
    t.x = fmaf (a.y, b.y, -t.y);
    t.x = fmaf (a.x, b.x, t.x);
    t.x = fmaf (a.y, b.x, t.x);
    t.x = fmaf (a.x, b.y, t.x);
    z.y = e = t.y + t.x;
    z.x = (t.y - e) + t.x;
    return z;
}

float2 div_dblflt (float2 a, float2 b)
{
    float2 t, z;
    float e, r;
    r = 1.0f / b.y;
    t.y = a.y * r;
    e = fmaf (b.y, -t.y, a.y);
    t.y = fmaf (r, e, t.y);
    t.x = fmaf (b.y, -t.y, a.y);
    t.x = a.x + t.x;
    t.x = fmaf (b.x, -t.y, t.x);
    e = r * t.x;
    t.x = fmaf (b.y, -e, t.x);
    t.x = fmaf (r, t.x, e);
    z.y = e = t.y + t.x;
    z.x = (t.y - e) + t.x;
    return z;
}

float2 sqrt_dblflt (float2 a)
{
    float2 t, z;
    float e, y, s, r;
    r = 1.0f / sqrtf (a.y);
    if (a.y == 0.0f) r = 0.0f;
    y = a.y * r;
    s = fmaf (y, -y, a.y);
    r = 0.5f * r;
    z.y = e = s + a.x;
    z.x = (s - e) + a.x;
    t.y = r * z.y;
    t.x = fmaf (r, z.y, -t.y);
    t.x = fmaf (r, z.x, t.x);
    r = y + t.y;
    s = (y - r) + t.y;
    s = s + t.x;
    z.y = e = r + s;
    z.x = (r - e) + s;
    return z;
}

#include "mpglue.h"  // for MP library

// Compute angle C of triangle with squared edges asq, bsq, csq
double lawOfCosines_ref (double asq, double bsq, double csq)
{
    struct mpNum mpAsq, mpBsq, mpCsq, mpTmp0, mpTmp1;
    double angle;

    mpDoubleToMp (asq, &mpAsq);        // asq
    mpDoubleToMp (bsq, &mpBsq);        // bsq
    mpDoubleToMp (csq, &mpCsq);        // csq
    mpAdd (&mpAsq, &mpBsq, &mpTmp0);   // asq+bsq
    mpSub (&mpTmp0, &mpCsq, &mpTmp0);  // asq+bsq-csq
    mpMul (&mpAsq, &mpBsq, &mpTmp1);   // asq*bsq
    mpSqrt (&mpTmp1, &mpTmp1);         // sqrt(asq*bsq)
    mpMul (mpTwo(), &mpTmp1, &mpTmp1); // 2*sqrt(asq*bsq)
    mpDiv (&mpTmp0, &mpTmp1, &mpTmp0); // (asq+bsq-csq)/(2*sqrt(asq*bsq))
    mpAcos (&mpTmp0, &mpTmp0);         // acos((asq+bsq-csq)/(2*sqrt(asq*bsq)))
    mpMpToDouble (&mpTmp0, &angle);    //
    return angle;
}

// George Marsaglia's KISS PRNG, period 2**123. Newsgroup sci.math, 21 Jan 1999
// Bug fix: Greg Rose, "KISS: A Bit Too Simple" http://eprint.iacr.org/2011/007
static uint32_t kiss_z=362436069, kiss_w=521288629;
static uint32_t kiss_jsr=123456789, kiss_jcong=380116160;
#define znew (kiss_z=36969*(kiss_z&65535)+(kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&65535)+(kiss_w>>16))
#define MWC  ((znew<<16)+wnew )
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17), \
              kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)

int randint (int a)
{
    return KISS % a;
}

#define MIN(x,y)     (fmin(x,y))
#define MAX(x,y)     (fmax(x,y))
#define MIN3(x,y,z)  MIN(x,MIN(y,z))
#define MAX3(x,y,z)  MAX(x,MAX(y,z))
#define MED3(x,y,z)  MIN(MAX(MIN(y,z),x),MAX(y,z))

#define ERR_LIMIT (0x1.0p-23)
#define SCALE     (0.00001357)

int main (void)
{
    mpInit();  // initialize MP library

    unsigned long long int count = 0;
    double A_ref = 0, B_ref = 0, C_ref = 0;
    double relerrA, relerrB, relerrC;
    float A = 0, B = 0, C = 0;

    do {
        double a, b, c, aa, bb, cc;
        float asq, bsq, csq;

        if ((count & 0xfff) == 0) printf ("\rcount=%llu", count);
        do {
            a = (randint (1 << 23) + 1) * SCALE;
            b = (randint (1 << 23) + 1) * SCALE;
            c = (randint (1 << 23) + 1) * SCALE;
            // sort edges by length, ascending
            aa = MIN3 (a, b, c);
            bb = MED3 (a, b, c);
            cc = MAX3 (a, b, c);
        } while ((aa + bb) <= (1.000001 * cc)); // ensure valid triangle

        asq = (float)(a * a);
        bsq = (float)(b * b);
        csq = (float)(c * c);

        // function under test
        A = angle_kahanf (bsq, csq, asq);
        B = angle_kahanf (csq, asq, bsq);
        C = angle_kahanf (asq, bsq, csq);

        // reference
        A_ref = lawOfCosines_ref ((double)bsq, (double)csq, (double)asq);
        B_ref = lawOfCosines_ref ((double)csq, (double)asq, (double)bsq);
        C_ref = lawOfCosines_ref ((double)asq, (double)bsq, (double)csq);

        // compute relative error compaored to reference
        relerrA = fabs ((A - A_ref) / A_ref);
        relerrB = fabs ((B - B_ref) / B_ref);
        relerrC = fabs ((C - C_ref) / C_ref);

        if (relerrA > ERR_LIMIT) {
            printf ("!!!! asq=%15.8e bsq=%15.8e csq=%15.8e A=%15.8e A_ref=%15.8e relerr=%15.8e\n", 
                    asq, bsq, csq, A, A_ref, relerrA);
        }
        if (relerrB > ERR_LIMIT) {
            printf ("!!!! asq=%15.8e bsq=%15.8e csq=%15.8e B=%15.8e B_ref=%15.8e relerr=%15.8e\n", 
                    asq, bsq, csq, B, B_ref, relerrB);
        }
        if (relerrC > ERR_LIMIT) {
            printf ("!!!! asq=%15.8e bsq=%15.8e csq=%15.8e C=%15.8e C_ref=%15.8e relerr=%15.8e\n", 
                    asq, bsq, csq, C, C_ref, relerrC);
        }
        count++;
    } while (count < 1000000);
    
    return EXIT_SUCCESS;
}