给定三角形的squared边长为
a, b, c > 0Kahan 的经常引用的方法仔细地重新排列括号，假设（非平方）边长。
另一种流行方法（第15页）假设边缘向量。
到目前为止，我最好的选择是使用余弦定理：
float angle = acos( (b+a-c) / (2.*sqrt(b*a)) );
但对于某些非退化输入，这将返回
0.0诚然，我生成有效输入的方式是通过计算边缘向量的平方长度，但为了这个问题的目的，我想忽略对原始边缘向量的访问。这是一个棘手的测试用例：
// 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;
真实角度约为：
4.9999997529193436e-08 2.3561944901923448 0.7853981133974508
如果我立即采取
sqrta,b,c0 3.14159 0
如果我应用余弦定理，我会得到：
0 2.35619 0.785398
哪个更好，但我不喜欢确切的
==0>0此外，如果我按非 2 的幂缩放
a,b,c0 2.02856 1.11303
有没有一种在数值上更准确的方法来直接从平方边长计算角度？
或者...
我是不是对浮动要求太多了？
我的浮动边长平方是否违反了某种三角形不等式？
我的工作实现，上面的数字在 mac os 上使用
clang++ main.cpp -g && ./a.out在重新排列浮点表达式、利用融合乘加 (FMA) 和补偿和进行了多个小时的实验之后，我得出的结论是，通过余弦定律进行的计算不会稳健地工作，即使切换到双
float 算术（也称为对算术）。
正如问题中已经指出的那样，卡汉在数值上优越的排列本身是不够的，因为在角度计算开始之前取平方根的需要已经注入了数值误差。然而，我观察到，在
doubledoublefloat对于下面的 C++11 代码，我假设目标平台支持 FMA，因此可用于加速双
float/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;
}