diff mbox series

[v4,05/12] math: Use an improved algorithm for hypotl (ldbl-128)

Message ID 20211203000103.737833-6-adhemerval.zanella@linaro.org
State Superseded
Headers show
Series Improve hypot | expand

Checks

Context Check Description
dj/TryBot-apply_patch success Patch applied to master at the time it was sent

Commit Message

Adhemerval Zanella Dec. 3, 2021, midnight UTC
This implementation is based on 'An Improved Algorithm for hypot(a,b)'
by Carlos F. Borges [1] using the MyHypot3 with the following changes:

  - Handle qNaN and sNaN.
  - Tune the 'widely varying operands' to avoid spurious underflow
    due the multiplication and fix the return value for upwards
    rounding mode.
  - Handle required underflow exception for subnormal results.

The main advantage of the new algorithm is its precision.  With a
random 1e9 input pairs in the range of [LDBL_MIN, LDBL_MAX], glibc
current implementation shows around 0.05% results with an error of
1 ulp (453266 results) while the new implementation only shows
0.0001% of total (1280).

Checked on aarch64-linux-gnu and x86_64-linux-gnu.

[1] https://arxiv.org/pdf/1904.09481.pdf
---
 sysdeps/ieee754/ldbl-128/e_hypotl.c | 226 ++++++++++++----------------
 1 file changed, 96 insertions(+), 130 deletions(-)

Comments

Wilco Dijkstra Dec. 6, 2021, 11:58 a.m. UTC | #1
Hi Adhemerval,


+  /* If ay is tiny, scale both inputs up.  */
+  if (__glibc_unlikely (ay < TINY_VAL))
+    {
+      if (__glibc_unlikely (ax >= ay / EPS))
+       return ax;

Why not return ax + ay like in the other special cases? It's more accurate for rounding
upwards, and it's best to keep the code as similar as possible as the double version.

Cheers,
Wilco
Adhemerval Zanella Dec. 6, 2021, 12:22 p.m. UTC | #2
On 06/12/2021 08:58, Wilco Dijkstra wrote:
> Hi Adhemerval,
> 
> 
> +  /* If ay is tiny, scale both inputs up.  */
> +  if (__glibc_unlikely (ay < TINY_VAL))
> +    {
> +      if (__glibc_unlikely (ax >= ay / EPS))
> +       return ax;
> 
> Why not return ax + ay like in the other special cases? It's more accurate for rounding
> upwards, and it's best to keep the code as similar as possible as the double version.

Ack, I have fix it.
diff mbox series

Patch

diff --git a/sysdeps/ieee754/ldbl-128/e_hypotl.c b/sysdeps/ieee754/ldbl-128/e_hypotl.c
index cd4fdbc4a6..022fa9aaf7 100644
--- a/sysdeps/ieee754/ldbl-128/e_hypotl.c
+++ b/sysdeps/ieee754/ldbl-128/e_hypotl.c
@@ -1,141 +1,107 @@ 
-/* e_hypotl.c -- long double version of e_hypot.c.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/* __ieee754_hypotl(x,y)
- *
- * Method :
- *	If (assume round-to-nearest) z=x*x+y*y
- *	has error less than sqrtl(2)/2 ulp, than
- *	sqrtl(z) has error less than 1 ulp (exercise).
- *
- *	So, compute sqrtl(x*x+y*y) with some care as
- *	follows to get the error below 1 ulp:
- *
- *	Assume x>y>0;
- *	(if possible, set rounding to round-to-nearest)
- *	1. if x > 2y  use
- *		x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
- *	where x1 = x with lower 64 bits cleared, x2 = x-x1; else
- *	2. if x <= 2y use
- *		t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
- *	where t1 = 2x with lower 64 bits cleared, t2 = 2x-t1,
- *	y1= y with lower 64 bits chopped, y2 = y-y1.
- *
- *	NOTE: scaling may be necessary if some argument is too
- *	      large or too tiny
- *
- * Special cases:
- *	hypotl(x,y) is INF if x or y is +INF or -INF; else
- *	hypotl(x,y) is NAN if x or y is NAN.
- *
- * Accuracy:
- *	hypotl(x,y) returns sqrtl(x^2+y^2) with error less
- *	than 1 ulps (units in the last place)
- */
+/* Euclidean distance function.  Long Double/Binary128 version.
+   Copyright (C) 2021 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+/* This implementation is based on 'An Improved Algorithm for hypot(a,b)' by
+   Carlos F. Borges [1] using the MyHypot3 with the following changes:
+
+   - Handle qNaN and sNaN.
+   - Tune the 'widely varying operands' to avoid spurious underflow
+     due the multiplication and fix the return value for upwards
+     rounding mode.
+   - Handle required underflow exception for subnormal results.
+
+   [1] https://arxiv.org/pdf/1904.09481.pdf  */
 
 #include <math.h>
 #include <math_private.h>
 #include <math-underflow.h>
 #include <libm-alias-finite.h>
 
+#define SCALE      L(0x1p-8303)
+#define LARGE_VAL  L(0x1.6a09e667f3bcc908b2fb1366ea95p+8191)
+#define TINY_VAL   L(0x1p-8191)
+#define EPS        L(0x1p-114)
+
+/* Hypot kernel. The inputs must be adjusted so that ax >= ay >= 0
+   and squaring ax, ay and (ax - ay) does not overflow or underflow.  */
+static inline _Float128
+kernel (_Float128 ax, _Float128 ay)
+{
+  _Float128 t1, t2;
+  _Float128 h = sqrtl (ax * ax + ay * ay);
+  if (h <= L(2.0) * ay)
+    {
+      _Float128 delta = h - ay;
+      t1 = ax * (L(2.0) * delta - ax);
+      t2 = (delta - L(2.0) * (ax - ay)) * delta;
+    }
+  else
+    {
+      _Float128 delta = h - ax;
+      t1 = L(2.0) * delta * (ax - L(2.0) * ay);
+      t2 = (L(4.0) * delta - ay) * ay + delta * delta;
+    }
+
+  h -= (t1 + t2) / (L(2.0) * h);
+  return h;
+}
+
 _Float128
 __ieee754_hypotl(_Float128 x, _Float128 y)
 {
-	_Float128 a,b,t1,t2,y1,y2,w;
-	int64_t j,k,ha,hb;
-
-	GET_LDOUBLE_MSW64(ha,x);
-	ha &= 0x7fffffffffffffffLL;
-	GET_LDOUBLE_MSW64(hb,y);
-	hb &= 0x7fffffffffffffffLL;
-	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
-	SET_LDOUBLE_MSW64(a,ha);	/* a <- |a| */
-	SET_LDOUBLE_MSW64(b,hb);	/* b <- |b| */
-	if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */
-	k=0;
-	if(ha > 0x5f3f000000000000LL) {	/* a>2**8000 */
-	   if(ha >= 0x7fff000000000000LL) {	/* Inf or NaN */
-	       uint64_t low;
-	       w = a+b;			/* for sNaN */
-	       if (issignaling (a) || issignaling (b))
-		 return w;
-	       GET_LDOUBLE_LSW64(low,a);
-	       if(((ha&0xffffffffffffLL)|low)==0) w = a;
-	       GET_LDOUBLE_LSW64(low,b);
-	       if(((hb^0x7fff000000000000LL)|low)==0) w = b;
-	       return w;
-	   }
-	   /* scale a and b by 2**-9600 */
-	   ha -= 0x2580000000000000LL;
-	   hb -= 0x2580000000000000LL;	k += 9600;
-	   SET_LDOUBLE_MSW64(a,ha);
-	   SET_LDOUBLE_MSW64(b,hb);
-	}
-	if(hb < 0x20bf000000000000LL) {	/* b < 2**-8000 */
-	    if(hb <= 0x0000ffffffffffffLL) {	/* subnormal b or 0 */
-		uint64_t low;
-		GET_LDOUBLE_LSW64(low,b);
-		if((hb|low)==0) return a;
-		t1=0;
-		SET_LDOUBLE_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */
-		b *= t1;
-		a *= t1;
-		k -= 16382;
-		GET_LDOUBLE_MSW64 (ha, a);
-		GET_LDOUBLE_MSW64 (hb, b);
-		if (hb > ha)
-		  {
-		    t1 = a;
-		    a = b;
-		    b = t1;
-		    j = ha;
-		    ha = hb;
-		    hb = j;
-		  }
-	    } else {		/* scale a and b by 2^9600 */
-		ha += 0x2580000000000000LL;	/* a *= 2^9600 */
-		hb += 0x2580000000000000LL;	/* b *= 2^9600 */
-		k -= 9600;
-		SET_LDOUBLE_MSW64(a,ha);
-		SET_LDOUBLE_MSW64(b,hb);
-	    }
-	}
-    /* medium size a and b */
-	w = a-b;
-	if (w>b) {
-	    t1 = 0;
-	    SET_LDOUBLE_MSW64(t1,ha);
-	    t2 = a-t1;
-	    w  = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
-	} else {
-	    a  = a+a;
-	    y1 = 0;
-	    SET_LDOUBLE_MSW64(y1,hb);
-	    y2 = b - y1;
-	    t1 = 0;
-	    SET_LDOUBLE_MSW64(t1,ha+0x0001000000000000LL);
-	    t2 = a - t1;
-	    w  = sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
-	}
-	if(k!=0) {
-	    uint64_t high;
-	    t1 = 1;
-	    GET_LDOUBLE_MSW64(high,t1);
-	    SET_LDOUBLE_MSW64(t1,high+(k<<48));
-	    w *= t1;
-	    math_check_force_underflow_nonneg (w);
-	    return w;
-	} else return w;
+  if (!isfinite(x) || !isfinite(y))
+    {
+      if ((isinf (x) || isinf (y))
+	  && !issignaling (x) && !issignaling (y))
+	return INFINITY;
+      return x + y;
+    }
+
+  x = fabsl (x);
+  y = fabsl (y);
+
+  _Float128 ax = x < y ? y : x;
+  _Float128 ay = x < y ? x : y;
+
+  /* If ax is huge, scale both inputs down.  */
+  if (__glibc_unlikely (ax > LARGE_VAL))
+    {
+      if (__glibc_unlikely (ay <= ax * EPS))
+	return ax + ay;
+
+      return kernel (ax * SCALE, ay * SCALE) / SCALE;
+    }
+
+  /* If ay is tiny, scale both inputs up.  */
+  if (__glibc_unlikely (ay < TINY_VAL))
+    {
+      if (__glibc_unlikely (ax >= ay / EPS))
+	return ax;
+
+      ax = kernel (ax / SCALE, ay / SCALE) * SCALE;
+      math_check_force_underflow_nonneg (ax);
+      return ax;
+    }
+
+  /* Common case: ax is not huge and ay is not tiny.  */
+  if (__glibc_unlikely (ay <= ax * EPS))
+    return ax + ay;
+
+  return kernel (ax, ay);
 }
 libm_alias_finite (__ieee754_hypotl, __hypotl)