Lossy Floating point: libdroundfast

TurboPFor-Integer-Compression/lib/ext/libdroundfast.c

@@ -0,0 +1,61 @@
+/*
+ * Copyright (c) 2019, CNES.
+ *
+ * This source code is licensed under MIT-style license (found in the
+ * COPYING file in the root directory of this source tree).
+ */
+//https://github.com/CNES/Digit_Rounding/blob/master/libdround/src/libdroundfast.c
+
+#include <math.h>
+
+#define LOG2_10		3.321928095		// log2(10)
+#define LOG10_2		0.301029996		// log10(2)
+
+#define SIGN(x)		( (x<0) ? -1 : 1 )
+
+const float TABLE[5][2] = {
+  {0.6, -LOG10_2},
+  {0.7,-0.221848749},
+  {0.8,-0.154901959},
+  {0.9,-0.096910013},
+  {1.0,-0.045757490},
+};
+
+/*
+ * Round the float value keeping nsd significant digits.
+ * Fast method that does not uses log10() function.
+ */
+double droundFast(double v, int nsd)
+{
+	// compute the number of digits before the decimal point of the input floating-point value v
+	// The value v is interpreted as v = 10^d + eps = 2^e + m
+	// with 0 <= m < 0.5
+	int e;
+	double m = frexp(v, &e);	// return the binary exponent e of the input value v = 2^e + m with 0 <= m < 0.5
+
+	// =============
+	// --- tabulated method ---
+	// tabulate the LOG10(m)
+	int i = 0;
+	while (TABLE[i][0] < m)
+	  {
+	    i++;
+	  }
+	float log10m = TABLE[i][1];
+	
+	// --- low precision method ---
+	// float log10m = -LOG10_2;
+	// =============
+
+	// convert the binary exponent to a number of digits: d = floor(e*log10(2) + log10(m)) + 1
+	int d = (int) floor(e*LOG10_2 + log10m) + 1;
+
+	// compute the power of the quantization step: q = 2^p
+	int p = (int) floor(LOG2_10 * (d - nsd));
+	// compute quantization step: q = 2^p
+	double q = ldexp(1, p);
+
+	// apply the quantization step depending on the bias
+	return SIGN(v) * (floor(fabs(v) / q) + 0.5) * q;
+}
+