Patchwork [PATCH1/2] Improve performance of sincosf

login
register
mail settings
Submitter Wilco Dijkstra
Date June 29, 2018, 12:20 p.m.
Message ID <DB5PR08MB1030632EF77F675D8FE8FBF3834E0@DB5PR08MB1030.eurprd08.prod.outlook.com>
Download mbox | patch
Permalink /patch/28140/
State New
Headers show

Comments

Wilco Dijkstra - June 29, 2018, 12:20 p.m.
Joseph Myers wrote:

> > +/* New sincosf.  */
>
> Comments describing something as new are never appropriate; they'll 
> always become out of date.

That comment was from early development, the 2nd patch removes it.
I've removed it from patch 1 as well.

>> +extern sincos_t sincosf_table[2] attribute_hidden;
>
> I'd expect this array, and thus all pointers to it, to be const; any 
> writable static / global data is suspect without a clear requirement for 
> it to be writable.

sincos_t itself is cost, however since that seems unclear I've added it to
the uses of sincos_t.

I also needed to change the patch because of recent changes to math barriers,
here is the updated version:


This patch is a complete rewrite of sinf, cosf and sincosf.  The new version
is significantly faster, as well as simple and accurate. 
The worst-case ULP is 0.56072, maximum relative error is 0.5303p-23 over all
4 billion inputs.  In non-nearest rounding modes the error is 1ULP.

The algorithm uses 3 main cases: small inputs which don't need argument
reduction, small inputs which need a simple range reduction and large inputs
requiring complex range reduction.  The code uses approximate integer
comparisons to quickly decide between these cases - on some targets this may
be slow, so this can be configured to use floating point comparisons.

The small range reducer uses a single reduction step to handle values up to
120.0.  It is fastest on targets which support inlined round instructions.

The large range reducer uses integer arithmetic for simplicity.  It does a
32x96 bit multiply to compute a 64-bit modulo result.  This is more than
accurate enough to handle the worst-case cancellation for values close to
an integer multiple of PI/4.  It could be further optimized, however it is
already much faster than necessary.

sincosf throughput gains on Cortex-A72:
* |x| < 0x1p-12 : 1.6x
* |x| < M_PI_4  : 1.7x
* |x| < 2 * M_PI: 1.5x
* |x| < 120.0   : 1.8x
* |x| < Inf     : 2.3x

On a benchmark with significant use of sincosf the overall speedup is >33%.

ChangeLog:
2018-06-29  Wilco Dijkstra  <wdijkstr@arm.com>

	* math/Makefile: Add s_sincosf_data.c.
	* sysdeps/aarch64/fpu/math_private.h (roundtoint): Use round.
	(converttoint): Use lround.
	* sysdeps/ieee754/flt-32/math_config.h (PREFER_FLOAT_COMPARISON): Define.
	* sysdeps/ieee754/flt-32/s_sincosf.h (abstop12): Add new function.
	(sincosf_poly): Likewise.
	(reduce_small): Likewise.
	(reduce_large): Likewise.
	* sysdeps/ieee754/flt-32/s_sincosf.c (sincosf): Rewrite.
	* sysdeps/ieee754/flt-32/s_sincosf_data.c: New file with sincosf data.
--
Joseph Myers - June 29, 2018, 8:25 p.m.
On Fri, 29 Jun 2018, Wilco Dijkstra wrote:

> 	* math/Makefile: Add s_sincosf_data.c.

You then need to add a dummy s_sincosf_data.c for architectures with their 
own sincosf implementation, to avoid wasting space in libm.so with unused 
data from s_sincosf_data.c.  That looks like ia64, m68k, x86_64.  (i686 
always builds in the generic version in case of a lack of SSE2 support, so 
no change is needed there, though I've no idea of the performance merits 
of the SSE2 version for i686 versus building the new generic one with SSE2 
enabled.)

> -/* Hack: only include the large arm_neon.h when needed.  */
> -#ifdef _MATH_CONFIG_H
> -# include <arm_neon.h>
> -
>  /* ACLE intrinsics for frintn and fcvtns instructions.  */
>  # define TOINT_INTRINSICS 1

Since you're removing the #ifdef you also need to remove the space between 
'#' and define (update the preprocessor indentation on what was 
conditional code).

>  #endif
>  
>  #include_next <math_private.h>
> -
> -#endif

So you seem to be moving the #include_next outside the multiple-include 
guards, so probably defeating the multiple-include optimization on this 
header.  Why?

> +#ifndef PREFER_FLOAT_COMPARISON
> +#define PREFER_FLOAT_COMPARISON 0
> +#endif

Missing preprocessor indentation, "# define".

> +/* The constants and polynomials for sine and cosine.  The 2nd entry
> +   computes -cos (x) rather than cos (x) to get negation for free.  */
> +const sincos_t sincosf_table[2] =

This should be __sincosf_table.  I'd have expected you to get 
linknamespace errors from testing your second patch with the name 
sincosf_table, which is in the user's namespace.  (They wouldn't have 
shown up with the first patch because sincosf isn't a standard function, 
but sinf and cosf are standard functions.)

Patch

diff --git a/math/Makefile b/math/Makefile
index 90b3b68916e12d8563ed57772013f4420415928b..e73ceb8d4e8d520c3047f63bf6c6d05856a89186 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -131,7 +131,7 @@  type-double-routines := branred doasin dosincos mpa mpatan2	\
 # float support
 type-float-suffix := f
 type-float-routines := k_rem_pio2f math_errf e_exp2f_data e_logf_data	\
-		       e_log2f_data e_powf_log2_data
+		       e_log2f_data e_powf_log2_data s_sincosf_data
 
 # _Float128 support
 type-float128-suffix := f128
diff --git a/sysdeps/aarch64/fpu/math_private.h b/sysdeps/aarch64/fpu/math_private.h
index fcd02c06543c0841cac9949c7ca052c042b3d85b..161e4cdca265edac2537eebb85a49ec13d65c7d6 100644
--- a/sysdeps/aarch64/fpu/math_private.h
+++ b/sysdeps/aarch64/fpu/math_private.h
@@ -21,6 +21,8 @@ 
 
 #include <fenv.h>
 #include <fpu_control.h>
+#include <stdint.h>
+#include <math.h>
 
 static __always_inline void
 libc_feholdexcept_aarch64 (fenv_t *envp)
@@ -298,26 +300,20 @@  libc_feresetround_noex_aarch64_ctx (struct rm_ctx *ctx)
 #define libc_feresetround_noexf_ctx	libc_feresetround_noex_aarch64_ctx
 #define libc_feresetround_noexl_ctx	libc_feresetround_noex_aarch64_ctx
 
-/* Hack: only include the large arm_neon.h when needed.  */
-#ifdef _MATH_CONFIG_H
-# include <arm_neon.h>
-
 /* ACLE intrinsics for frintn and fcvtns instructions.  */
 # define TOINT_INTRINSICS 1
 
 static inline double_t
 roundtoint (double_t x)
 {
-  return vget_lane_f64 (vrndn_f64 (vld1_f64 (&x)), 0);
+  return round (x);
 }
 
 static inline uint64_t
 converttoint (double_t x)
 {
-  return vcvtnd_s64_f64 (x);
+  return (uint64_t) lround (x);
 }
 #endif
 
 #include_next <math_private.h>
-
-#endif
diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
index c4def9bbc1898d71865572d0084b219034bf3805..b68fc47acc811d9e9844faac1ca71ec3b4d7de41 100644
--- a/sysdeps/ieee754/flt-32/math_config.h
+++ b/sysdeps/ieee754/flt-32/math_config.h
@@ -46,6 +46,9 @@ 
 #ifndef TOINT_SHIFT
 # define TOINT_SHIFT 1
 #endif
+#ifndef PREFER_FLOAT_COMPARISON
+#define PREFER_FLOAT_COMPARISON 0
+#endif
 
 static inline uint32_t
 asuint (float f)
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h
index 35b5eee536ef0ca55e68e3d3a091486dca39d36a..88785aa5ef895bd39b1340ff6f85257669760c63 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.h
+++ b/sysdeps/ieee754/flt-32/s_sincosf.h
@@ -16,6 +16,10 @@ 
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <stdint.h>
+#include <math.h>
+#include "math_config.h"
+
 /* Chebyshev constants for cos, range -PI/4 - PI/4.  */
 static const double C0 = -0x1.ffffffffe98aep-2;
 static const double C1 =  0x1.55555545c50c7p-5;
@@ -153,3 +157,118 @@  reduced_cos (double theta, unsigned int n)
     }
   return sign * cx;
 }
+
+
+/* PI * 2^-64.  */
+static const double pi64 = 0x1.921FB54442D18p-62;
+/* PI / 4.  */
+static const double pio4 = 0x1.921FB54442D18p-1;
+
+typedef struct
+{
+  double sign[4];
+  double hpi_inv, hpi, c0, c1, c2, c3, c4, s1, s2, s3;
+} sincos_t;
+
+extern const sincos_t sincosf_table[2] attribute_hidden;
+
+extern const uint32_t inv_pio4[] attribute_hidden;
+
+/* abstop12 assumes floating point reinterpret is fast by default.
+   If floating point comparisons are faster, define PREFER_FLOAT_COMPARISON.  */
+#if PREFER_FLOAT_COMPARISON
+static inline float
+abstop12 (float x)
+{
+  return fabsf (x);
+}
+#else
+static inline uint32_t
+abstop12 (float x)
+{
+  return (asuint (x) >> 20) & 0x7ff;
+}
+#endif
+
+/* Compute the sine and cosine of inputs X and X2 (X squared), using the
+   polynomial P and store the results in SINP and COSP.  N is the quadrant,
+   if odd the cosine and sine polynomials are swapped.  */
+static inline void
+sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp,
+	      float *cosp)
+{
+  double x3, x4, x5, x6, s, c, c1, c2, s1;
+
+  x4 = x2 * x2;
+  x3 = x2 * x;
+  c2 = p->c3 + x2 * p->c4;
+  s1 = p->s2 + x2 * p->s3;
+
+  /* Swap sin/cos result based on quadrant.  */
+  float *tmp = (n & 1 ? cosp : sinp);
+  cosp = (n & 1 ? sinp : cosp);
+  sinp = tmp;
+
+  c1 = p->c0 + x2 * p->c1;
+  x5 = x3 * x2;
+  x6 = x4 * x2;
+
+  s = x + x3 * p->s1;
+  c = c1 + x4 * p->c2;
+
+  *sinp = s + x5 * s1;
+  *cosp = c + x6 * c2;
+}
+
+/* Fast range reduction using single multiply-subtract.  Return the modulo of
+   X as a value between -PI/4 and PI/4 and store the quadrant in NP.
+   The values for PI/2 and 2/PI are accessed via P.  Since PI/2 as a double
+   is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
+   only 2 multiplies are required and the result is accurate for |X| <= 120.0.
+   Use round/lround if inlined, otherwise convert to int.  To avoid inaccuracies
+   introduced by truncating negative values, compute the quadrant * 2^24.  */
+static inline double
+reduce_fast (double x, const sincos_t *p, int *np)
+{
+  double r;
+#if TOINT_INTRINSICS
+  r = x * p->hpi_inv;
+  *np = converttoint (r);
+  return x - roundtoint (r) * p->hpi;
+#else
+  r = x * p->hpi_inv;
+  int n = ((int32_t)r + 0x800000) >> 24;
+  *np = n;
+  return x - n * p->hpi;
+#endif
+}
+
+/* Reduce the range of XI to a multiple of PI/4 using fast integer arithmetic.
+   XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
+   Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
+   Reduction uses a table of 4/PI with 192 bits of precision.  A 32x96->128 bit
+   multiply computes the exact 2.62-bit fixed-point modulo.  Since the result
+   can have at most 29 leading zeros after the binary point, the double
+   precision result is accurate to 33 bits.  */
+static inline double
+reduce_large (uint32_t xi, int *np)
+{
+  const uint32_t *arr = &inv_pio4[(xi >> 26) & 15];
+  int shift = (xi >> 23) & 7;
+  uint64_t n, res0, res1, res2;
+
+  xi = (xi & 0xffffff) | 0x800000;
+  xi <<= shift;
+
+  res0 = xi * arr[0];
+  res1 = (uint64_t)xi * arr[4];
+  res2 = (uint64_t)xi * arr[8];
+  res0 = (res2 >> 32) | (res0 << 32);
+  res0 += res1;
+
+  n = (res0 + (1ULL << 61)) >> 62;
+  res0 -= n << 62;
+  double x = (int64_t)res0;
+  *np = n;
+  return x * pi64;
+}
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c
index d4a5a1b22c2f9a5cf607cd2d31157dd1dbd14e9a..f8fec38f0fd7c702dfe7e3a6abcc3b19f228cedd 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.c
+++ b/sysdeps/ieee754/flt-32/s_sincosf.c
@@ -1,5 +1,5 @@ 
 /* Compute sine and cosine of argument.
-   Copyright (C) 2017-2018 Free Software Foundation, Inc.
+   Copyright (C) 2018 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
@@ -17,9 +17,11 @@ 
    <http://www.gnu.org/licenses/>.  */
 
 #include <errno.h>
+#include <stdint.h>
 #include <math.h>
-#include <math_private.h>
+#include <math-barriers.h>
 #include <libm-alias-float.h>
+#include "math_config.h"
 #include "s_sincosf.h"
 
 #ifndef SINCOSF
@@ -28,141 +30,72 @@ 
 # define SINCOSF_FUNC SINCOSF
 #endif
 
+/* Fast sincosf implementation.  Worst-case ULP is 0.56072, maximum relative
+   error is 0.5303p-23.  A single-step signed range reduction is used for
+   small values.  Large inputs have their range reduced using fast integer
+   arithmetic.
+*/
 void
-SINCOSF_FUNC (float x, float *sinx, float *cosx)
+SINCOSF_FUNC (float y, float *sinp, float *cosp)
 {
-  double cx;
-  double theta = x;
-  double abstheta = fabs (theta);
-  /* If |x|< Pi/4.  */
-  if (isless (abstheta, M_PI_4))
+  double x = y;
+  double s;
+  int n;
+  const sincos_t *p = &sincosf_table[0];
+
+  if (abstop12 (y) < abstop12 (pio4))
+    {
+      double x2 = x * x;
+
+      if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
+      {
+	/* Force underflow for tiny y.  */
+	if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
+	  math_force_eval ((float)x2);
+	*sinp = y;
+	*cosp = 1.0f;
+	return;
+      }
+
+      sincosf_poly (x, x2, p, 0, sinp, cosp);
+    }
+  else if (abstop12 (y) < abstop12 (120.0f))
     {
-      if (abstheta >= 0x1p-5) /* |x| >= 2^-5.  */
-	{
-	  const double theta2 = theta * theta;
-	  /* Chebyshev polynomial of the form for sin and cos.  */
-	  cx = C3 + theta2 * C4;
-	  cx = C2 + theta2 * cx;
-	  cx = C1 + theta2 * cx;
-	  cx = C0 + theta2 * cx;
-	  cx = 1.0 + theta2 * cx;
-	  *cosx = cx;
-	  cx = S3 + theta2 * S4;
-	  cx = S2 + theta2 * cx;
-	  cx = S1 + theta2 * cx;
-	  cx = S0 + theta2 * cx;
-	  cx = theta + theta * theta2 * cx;
-	  *sinx = cx;
-	}
-      else if (abstheta >= 0x1p-27)     /* |x| >= 2^-27.  */
-	{
-	  /* A simpler Chebyshev approximation is close enough for this range:
-	     for sin: x+x^3*(SS0+x^2*SS1)
-	     for cos: 1.0+x^2*(CC0+x^3*CC1).  */
-	  const double theta2 = theta * theta;
-	  cx = CC0 + theta * theta2 * CC1;
-	  cx = 1.0 + theta2 * cx;
-	  *cosx = cx;
-	  cx = SS0 + theta2 * SS1;
-	  cx = theta + theta * theta2 * cx;
-	  *sinx = cx;
-	}
-      else
-	{
-	  /* Handle some special cases.  */
-	  if (theta)
-	    *sinx = theta - (theta * SMALL);
-	  else
-	    *sinx = theta;
-	  *cosx = 1.0 - abstheta;
-	}
+      x = reduce_fast (x, p, &n);
+
+      /* Setup the signs for sin and cos.  */
+      s = p->sign[n & 3];
+
+      if (n & 2)
+	p = &sincosf_table[1];
+
+      sincosf_poly (x * s, x * x, p, n, sinp, cosp);
     }
-  else                          /* |x| >= Pi/4.  */
+  else if (__glibc_likely (abstop12 (y) < abstop12 (INFINITY)))
     {
-      unsigned int signbit = isless (x, 0);
-      if (isless (abstheta, 9 * M_PI_4))        /* |x| < 9*Pi/4.  */
-	{
-	  /* There are cases where FE_UPWARD rounding mode can
-	     produce a result of abstheta * inv_PI_4 == 9,
-	     where abstheta < 9pi/4, so the domain for
-	     pio2_table must go to 5 (9 / 2 + 1).  */
-	  unsigned int n = (abstheta * inv_PI_4) + 1;
-	  theta = abstheta - pio2_table[n / 2];
-	  *sinx = reduced_sin (theta, n, signbit);
-	  *cosx = reduced_cos (theta, n);
-	}
-      else if (isless (abstheta, INFINITY))
-	{
-	  if (abstheta < 0x1p+23)     /* |x| < 2^23.  */
-	    {
-	      unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
-	      double x = n / 2;
-	      theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
-	      /* Argument reduction needed.  */
-	      *sinx = reduced_sin (theta, n, signbit);
-	      *cosx = reduced_cos (theta, n);
-	    }
-	  else                  /* |x| >= 2^23.  */
-	    {
-	      x = fabsf (x);
-	      int exponent;
-	      GET_FLOAT_WORD (exponent, x);
-	      exponent
-	        = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS;
-	      exponent += 3;
-	      exponent /= 28;
-	      double a = invpio4_table[exponent] * x;
-	      double b = invpio4_table[exponent + 1] * x;
-	      double c = invpio4_table[exponent + 2] * x;
-	      double d = invpio4_table[exponent + 3] * x;
-	      uint64_t l = a;
-	      l &= ~0x7;
-	      a -= l;
-	      double e = a + b;
-	      l = e;
-	      e = a - l;
-	      if (l & 1)
-	        {
-	          e -= 1.0;
-	          e += b;
-	          e += c;
-	          e += d;
-	          e *= M_PI_4;
-		  *sinx = reduced_sin (e, l + 1, signbit);
-		  *cosx = reduced_cos (e, l + 1);
-	        }
-	      else
-		{
-		  e += b;
-		  e += c;
-		  e += d;
-		  if (e <= 1.0)
-		    {
-		      e *= M_PI_4;
-		      *sinx = reduced_sin (e, l + 1, signbit);
-		      *cosx = reduced_cos (e, l + 1);
-		    }
-		  else
-		    {
-		      l++;
-		      e -= 2.0;
-		      e *= M_PI_4;
-		      *sinx = reduced_sin (e, l + 1, signbit);
-		      *cosx = reduced_cos (e, l + 1);
-		    }
-		}
-	    }
-	}
-      else
-	{
-	  int32_t ix;
-	  /* High word of x.  */
-	  GET_FLOAT_WORD (ix, abstheta);
-	  /* sin/cos(Inf or NaN) is NaN.  */
-	  *sinx = *cosx = x - x;
-	  if (ix == 0x7f800000)
-	    __set_errno (EDOM);
-	}
+      uint32_t xi = asuint (y);
+      int sign = xi >> 31;
+
+      x = reduce_large (xi, &n);
+
+      /* Setup signs for sin and cos - include original sign.  */
+      s = p->sign[(n + sign) & 3];
+
+      if ((n + sign) & 2)
+	p = &sincosf_table[1];
+
+      sincosf_poly (x * s, x * x, p, n, sinp, cosp);
+    }
+  else
+    {
+      /* Return NaN if Inf or NaN for both sin and cos.  */
+      *sinp = *cosp = y - y;
+#if WANT_ERRNO
+      /* Needed to set errno for +-Inf, the add is a hack to work
+	 around a gcc register allocation issue: just passing y
+	 affects code generation in the fast path.  */
+      __math_invalidf (y + y);
+#endif
     }
 }
 
diff --git a/sysdeps/ieee754/flt-32/s_sincosf_data.c b/sysdeps/ieee754/flt-32/s_sincosf_data.c
new file mode 100644
index 0000000000000000000000000000000000000000..c3dac5fe2ead99e279df26e5d2bc5a46d8d425d2
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/s_sincosf_data.c
@@ -0,0 +1,74 @@ 
+/* Compute sine and cosine of argument.
+   Copyright (C) 2018 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
+   <http://www.gnu.org/licenses/>.  */
+
+#include <stdint.h>
+#include <math.h>
+#include "math_config.h"
+#include "s_sincosf.h"
+
+/* The constants and polynomials for sine and cosine.  The 2nd entry
+   computes -cos (x) rather than cos (x) to get negation for free.  */
+const sincos_t sincosf_table[2] =
+{
+  {
+    { 1.0, -1.0, -1.0, 1.0 },
+#if TOINT_INTRINSICS
+    0x1.45F306DC9C883p-1,
+#else
+    0x1.45F306DC9C883p+23,
+#endif
+    0x1.921FB54442D18p0,
+    0x1p0,
+    -0x1.ffffffd0c621cp-2,
+    0x1.55553e1068f19p-5,
+    -0x1.6c087e89a359dp-10,
+    0x1.99343027bf8c3p-16,
+    -0x1.555545995a603p-3,
+    0x1.1107605230bc4p-7,
+    -0x1.994eb3774cf24p-13
+  },
+  {
+    { 1.0, -1.0, -1.0, 1.0 },
+#if TOINT_INTRINSICS
+    0x1.45F306DC9C883p-1,
+#else
+    0x1.45F306DC9C883p+23,
+#endif
+    0x1.921FB54442D18p0,
+    -0x1p0,
+    0x1.ffffffd0c621cp-2,
+    -0x1.55553e1068f19p-5,
+    0x1.6c087e89a359dp-10,
+    -0x1.99343027bf8c3p-16,
+    -0x1.555545995a603p-3,
+    0x1.1107605230bc4p-7,
+    -0x1.994eb3774cf24p-13
+  }
+};
+
+/* Table with 4/PI to 192 bit precision.  To avoid unaligned accesses
+   only 8 new bits are added per entry, making the table 4 times larger.  */
+const uint32_t inv_pio4[24] =
+{
+  0xa2,       0xa2f9,	  0xa2f983,   0xa2f9836e,
+  0xf9836e4e, 0x836e4e44, 0x6e4e4415, 0x4e441529,
+  0x441529fc, 0x1529fc27, 0x29fc2757, 0xfc2757d1,
+  0x2757d1f5, 0x57d1f534, 0xd1f534dd, 0xf534ddc0,
+  0x34ddc0db, 0xddc0db62, 0xc0db6295, 0xdb629599,
+  0x6295993c, 0x95993c43, 0x993c4390, 0x3c439041
+};