angle-comparison.cpp

#include <cmath>
#include <cstdio>

float angle_via_law_of_cosines(float a, float b, float c)
{
 return std::acos( (b+a-c) / (float(2)*std::sqrt(b*a)) );
}

float angle_via_kahan(float a, float b, float c)
{
  float A = std::sqrt(a);
  float B = std::sqrt(b);
  float C = std::sqrt(c);
  // If necessary, swap a and b so that a ≥ b . 
  if(A<B){ 
    float tmp = B;
    B = A;
    A = tmp;
  }
  float mu;
  if(B>=C && C>=0)
  {
    mu = C-(A-B);
  }else if(C>B && B>=0)
  {
    mu = B-(A-C);
  }else
  {
    return 0.0/0.0;
  }
  return 2.0*std::atan(std::sqrt(((A-B)+C)*mu/((A+(B+C))*((A-C)+B)) ));
}

int main()
{
  // 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;


  printf("Law of Cosines: %0.17g %0.17g %0.17g\n",
    angle_via_law_of_cosines(b,c,a),
    angle_via_law_of_cosines(c,a,b),
    angle_via_law_of_cosines(a,b,c));

  printf("Kahan's method: %0.17g %0.17g %0.17g\n",
    angle_via_kahan(b,c,a),
    angle_via_kahan(c,a,b),
    angle_via_kahan(a,b,c));
}