Fix accuracy in angle computation (#224)

* Fix accuracy in angle computation

* Add unit test

* Enhance comparisons and expose QD 2*pi

CHANGES.rst

@@ -14,6 +14,13 @@ Release Notes
 - Fix problem reported where angle of 180 degrees results in an
   error with the C version of the code. [#223]
 
+- Use a different algorithm for square root in the ``qd`` library that
+  seems to be less prone to accuracy loss. This helps the bug in the
+  ``math_util.c`` due to which ``angle()`` returns a NaN for
+  coplanar vectors where the angle between the surface points should be
+  180 degrees and enhances the solution from #223. Also, use ``qd`` epsilon
+  instead of arbitrary value for rounding error check. [#224]
+
 
 1.2.22 (04-January-2022)
 ========================

libqd/include/qd/c_dd.h

@@ -90,6 +90,8 @@ void c_dd_comp(const double *a, const double *b, int *result);
 void c_dd_comp_dd_d(const double *a, double b, int *result);
 void c_dd_comp_d_dd(double a, const double *b, int *result);
 void c_dd_pi(double *a);
+void c_dd_2pi(double *a);
+double c_dd_epsilon(void);
 
 #ifdef __cplusplus
 }

libqd/include/qd/c_qd.h

@@ -7,7 +7,7 @@
  *
  * Copyright (c) 2000-2001
  *
- * Contains C wrapper function prototypes for quad-double precision 
+ * Contains C wrapper function prototypes for quad-double precision
  * arithmetic.  This can also be used from fortran code.
  */
 #ifndef _QD_C_QD_H
@@ -65,7 +65,7 @@ void c_qd_copy(const double *a, double *b);
 void c_qd_copy_dd(const double *a, double *b);
 void c_qd_copy_d(double a, double *b);
 
-void c_qd_sqrt(const double *a, double *b);
+int c_qd_sqrt(const double *a, double *b);
 void c_qd_sqr(const double *a, double *b);
 
 void c_qd_abs(const double *a, double *b);
@@ -111,6 +111,8 @@ void c_qd_comp(const double *a, const double *b, int *result);
 void c_qd_comp_qd_d(const double *a, double b, int *result);
 void c_qd_comp_d_qd(double a, const double *b, int *result);
 void c_qd_pi(double *a);
+void c_qd_2pi(double *a);
+double c_qd_epsilon(void);
 
 #ifdef __cplusplus
 }

libqd/include/qd/qd_real.h

@@ -9,8 +9,8 @@
  *
  * Quad-double precision (>= 212-bit significand) floating point arithmetic
  * package, written in ANSI C++, taking full advantage of operator overloading.
- * Uses similar techniques as that of David Bailey's double-double package 
- * and that of Jonathan Shewchuk's adaptive precision floating point 
+ * Uses similar techniques as that of David Bailey's double-double package
+ * and that of Jonathan Shewchuk's adaptive precision floating point
  * arithmetic package.  See
  *
  *   http://www.nersc.gov/~dhbailey/mpdist/mpdist.html
@@ -120,16 +120,16 @@ struct QD_API qd_real {
   static qd_real rand(void);
 
   void to_digits(char *s, int &expn, int precision = _ndigits) const;
-  void write(char *s, int len, int precision = _ndigits, 
+  void write(char *s, int len, int precision = _ndigits,
       bool showpos = false, bool uppercase = false) const;
-  std::string to_string(int precision = _ndigits, int width = 0, 
-      std::ios_base::fmtflags fmt = static_cast<std::ios_base::fmtflags>(0), 
+  std::string to_string(int precision = _ndigits, int width = 0,
+      std::ios_base::fmtflags fmt = static_cast<std::ios_base::fmtflags>(0),
       bool showpos = false, bool uppercase = false, char fill = ' ') const;
   static int read(const char *s, qd_real &a);
 
   /* Debugging methods */
   void dump(const std::string &name = "", std::ostream &os = std::cerr) const;
-  void dump_bits(const std::string &name = "", 
+  void dump_bits(const std::string &name = "",
                  std::ostream &os = std::cerr) const;
 
   static qd_real debug_rand();
@@ -150,7 +150,7 @@ namespace std {
 }
 
 QD_API qd_real polyeval(const qd_real *c, int n, const qd_real &x);
-QD_API qd_real polyroot(const qd_real *c, int n, 
+QD_API qd_real polyroot(const qd_real *c, int n,
     const qd_real &x0, int max_iter = 64, double thresh = 0.0);
 
 QD_API qd_real qdrand(void);
@@ -190,6 +190,7 @@ QD_API qd_real operator/(double a, const qd_real &b);
 
 QD_API qd_real sqr(const qd_real &a);
 QD_API qd_real sqrt(const qd_real &a);
+QD_API qd_real fsqrt(const qd_real &a, int &flag);
 QD_API qd_real pow(const qd_real &a, int n);
 QD_API qd_real pow(const qd_real &a, const qd_real &b);
 QD_API qd_real npwr(const qd_real &a, int n);

libqd/src/c_dd.cpp

@@ -283,7 +283,7 @@ void c_dd_comp(const double *a, const double *b, int *result) {
     *result = -1;
   else if (aa > bb)
     *result = 1;
-  else 
+  else
     *result = 0;
 }
 
@@ -293,7 +293,7 @@ void c_dd_comp_dd_d(const double *a, double b, int *result) {
     *result = -1;
   else if (aa > bb)
     *result = 1;
-  else 
+  else
     *result = 0;
 }
 
@@ -303,12 +303,20 @@ void c_dd_comp_d_dd(double a, const double *b, int *result) {
     *result = -1;
   else if (aa > bb)
     *result = 1;
-  else 
+  else
     *result = 0;
 }
 
 void c_dd_pi(double *a) {
   TO_DOUBLE_PTR(dd_real::_pi, a);
 }
 
+void c_dd_2pi(double *a) {
+  TO_DOUBLE_PTR(dd_real::_2pi, a);
+}
+
+double c_dd_epsilon(void) {
+    return (double) std::numeric_limits<dd_real>::epsilon();
+}
+
 }

libqd/src/c_qd.cpp

@@ -237,11 +237,14 @@ void c_qd_copy_d(double a, double *b) {
 }
 
 
-void c_qd_sqrt(const double *a, double *b) {
+int c_qd_sqrt(const double *a, double *b) {
+  int flag;
   qd_real bb;
-  bb = sqrt(qd_real(a));
+  bb = fsqrt(qd_real(a), flag);
   TO_DOUBLE_PTR(bb, b);
+  return flag;
 }
+
 void c_qd_sqr(const double *a, double *b) {
   qd_real bb;
   bb = sqr(qd_real(a));
@@ -419,7 +422,7 @@ void c_qd_comp(const double *a, const double *b, int *result) {
     *result = -1;
   else if (aa > bb)
     *result = 1;
-  else 
+  else
     *result = 0;
 }
 
@@ -429,7 +432,7 @@ void c_qd_comp_qd_d(const double *a, double b, int *result) {
     *result = -1;
   else if (aa > b)
     *result = 1;
-  else 
+  else
     *result = 0;
 }
 
@@ -439,12 +442,21 @@ void c_qd_comp_d_qd(double a, const double *b, int *result) {
     *result = -1;
   else if (a > bb)
     *result = 1;
-  else 
+  else
     *result = 0;
 }
 
 void c_qd_pi(double *a) {
   TO_DOUBLE_PTR(qd_real::_pi, a);
 }
 
+void c_qd_2pi(double *a) {
+  TO_DOUBLE_PTR(qd_real::_2pi, a);
+}
+
+double c_qd_epsilon(void) {
+    return (double) std::numeric_limits<qd_real>::epsilon();
+}
+
+
 }