Source-Changes-HG archive
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index][Old Index]
[src/trunk]: src/sys/dev/dtv Instead of returning an integer between 0 and 31...
details: https://anonhg.NetBSD.org/src/rev/5d110ed47af6
branches: trunk
changeset: 767354:5d110ed47af6
user: apb <apb%NetBSD.org@localhost>
date: Sat Jul 16 22:30:26 2011 +0000
description:
Instead of returning an integer between 0 and 31 (scaled by 1<<24), use
(0.5 + x/2 - 1/x) as an approximation to log2(x) for x from 1 to 2, and
scale the input to fit this range. Now the error is always less than
0.2%.
Also add an test program, hidden behind #ifdef TEST_DTV_MATH, to print
a table of expected and actual results, and the errors.
diffstat:
sys/dev/dtv/dtv_math.c | 169 ++++++++++++++++++++++++++++++++++++++++--------
1 files changed, 140 insertions(+), 29 deletions(-)
diffs (199 lines):
diff -r 4442557e0b01 -r 5d110ed47af6 sys/dev/dtv/dtv_math.c
--- a/sys/dev/dtv/dtv_math.c Sat Jul 16 22:16:59 2011 +0000
+++ b/sys/dev/dtv/dtv_math.c Sat Jul 16 22:30:26 2011 +0000
@@ -1,4 +1,4 @@
-/* $NetBSD: dtv_math.c,v 1.2 2011/07/16 16:13:13 jmcneill Exp $ */
+/* $NetBSD: dtv_math.c,v 1.3 2011/07/16 22:30:26 apb Exp $ */
/*-
* Copyright (c) 2011 Alan Barrett <apb%NetBSD.org@localhost>
@@ -27,50 +27,161 @@
*/
#include <sys/cdefs.h>
-__KERNEL_RCSID(0, "$NetBSD: dtv_math.c,v 1.2 2011/07/16 16:13:13 jmcneill Exp $");
+__KERNEL_RCSID(0, "$NetBSD: dtv_math.c,v 1.3 2011/07/16 22:30:26 apb Exp $");
#include <sys/types.h>
#include <sys/bitops.h>
-#include <dev/dtv/dtvif.h>
-
-
-#define LOG10_2_x24 5050445 /* floor(log10(2.0) * 2**24 */
-
/*
* dtv_intlog10 -- return an approximation to log10(x) * 1<<24,
* using integer arithmetic.
*
- * As a special case, returns 0 when x == 0.
+ * As a special case, returns 0 when x == 0. The mathematical
+ * result is -infinity.
*
- * Results should be approximately as follows, bearing in
- * mind that this function returns only an approximation
- * to the exact results.
+ * This function uses 0.5 + x/2 - 1/x as an approximation to
+ * log2(x) for x in the range [1.0, 2.0], and scales the input value
+ * to fit this range. The resulting error is always better than
+ * 0.2%.
+ *
+ * Here's a table of the desired and actual results, as well
+ * as the absolute and relative errors, for several values of x.
*
- * dtv_intlog10(0) = 0 (special case; the mathematical value is undefined)
- * dtv_intlog10(1) = 0
- * dtv_intlog10(2) = 5050445 (approx 0.30102999 * 2**24)
- * dtv_intlog10(10) = 16777216 (1.0 * 2**24)
- * dtv_intlog10(100) = 33554432 (2.0 * 2**24)
- * dtv_intlog10(1000) = 50331648 (3.0 * 2**24)
- * dtv_intlog10(10000) = 67108864 (4.0 * 2**24)
- * dtv_intlog10(100000) = 83886080 (5.0 * 2**24)
- * dtv_intlog10(1000000) = 100663296 (6.0 * 2**24)
- * dtv_intlog10(10000000) = 117440512 (7.0 * 2**24)
- * dtv_intlog10(100000000) = 134217728 (8.0 * 2**24)
- * dtv_intlog10(1000000000) = 150994944 (9.0 * 2**24)
- * dtv_intlog10(4294967295) = 161614248 (approx 9.63295986 * 2**24)
+ * x desired actual err_abs err_rel
+ * 0 0 0 +0 +0.00000
+ * 1 0 0 +0 +0.00000
+ * 2 5050445 5050122 -323 -0.00006
+ * 3 8004766 7996348 -8418 -0.00105
+ * 4 10100890 10100887 -3 -0.00000
+ * 5 11726770 11741823 +15053 +0.00128
+ * 6 13055211 13046470 -8741 -0.00067
+ * 7 14178392 14158860 -19532 -0.00138
+ * 8 15151335 15151009 -326 -0.00002
+ * 9 16009532 16028061 +18529 +0.00116
+ * 10 16777216 16792588 +15372 +0.00092
+ * 11 17471670 17475454 +3784 +0.00022
+ * 12 18105656 18097235 -8421 -0.00047
+ * 13 18688868 18672077 -16791 -0.00090
+ * 14 19228837 19209625 -19212 -0.00100
+ * 15 19731537 19717595 -13942 -0.00071
+ * 16 20201781 20201774 -7 -0.00000
+ * 20 21827661 21842710 +15049 +0.00069
+ * 24 23156102 23147357 -8745 -0.00038
+ * 30 24781982 24767717 -14265 -0.00058
+ * 40 26878106 26893475 +15369 +0.00057
+ * 60 29832427 29818482 -13945 -0.00047
+ * 100 33554432 33540809 -13623 -0.00041
+ * 1000 50331648 50325038 -6610 -0.00013
+ * 10000 67108864 67125985 +17121 +0.00026
+ * 100000 83886080 83875492 -10588 -0.00013
+ * 1000000 100663296 100652005 -11291 -0.00011
+ * 10000000 117440512 117458739 +18227 +0.00016
+ * 100000000 134217728 134210175 -7553 -0.00006
+ * 1000000000 150994944 150980258 -14686 -0.00010
+ * 4294967295 161614248 161614192 -56 -0.00000
*/
uint32_t
dtv_intlog10(uint32_t x)
{
+ uint32_t ilog2x;
+ uint32_t t;
+ uint32_t t1;
+
if (__predict_false(x == 0))
return 0;
+
/*
- * all we do is find log2(x), as an integer between 0 and 31,
- * and scale it. Thus, there are only 32 values that this
- * function will ever return. To do a better job, we would
- * need a lookup table and interpolation.
+ * find ilog2x = floor(log2(x)), as an integer in the range [0,31].
+ */
+ ilog2x = ilog2(x);
+
+ /*
+ * Set "t" to the result of shifting x left or right
+ * until the most significant bit that was actually set
+ * moves into the 1<<24 position.
+ *
+ * Now we can think of "t" as representing
+ * x / 2**(floor(log2(x))),
+ * as a fixed-point value with 8 integer bits and 24 fraction bits.
+ *
+ * This value is in the semi-closed interval [1.0, 2.0)
+ * when interpreting it as a fixed-point number, or in the
+ * interval [0x01000000, 0x01ffffff] when examining the
+ * underlying uint32_t representation.
+ */
+ t = (ilog2x > 24 ? x >> (ilog2x - 24) : x << (24 - ilog2x));
+
+ /*
+ * Calculate "t1 = 1 / t" in the 8.24 fixed-point format.
+ * This value is in the interval [0.5, 1.0]
+ * when interpreting it as a fixed-point number, or in the
+ * interval [0x00800000, 0x01000000] when examining the
+ * underlying uint32_t representation.
+ *
+ */
+ t1 = ((uint64_t)1 << 48) / t;
+
+ /*
+ * Calculate "t = ilog2x + t/2 - t1 + 0.5" in the 8.24
+ * fixed-point format.
+ *
+ * If x is a power of 2, then t is now exactly equal to log2(x)
+ * when interpreting it as a fixed-point number, or exactly
+ * log2(x) << 24 when examining the underlying uint32_t
+ * representation.
+ *
+ * If x is not a power of 2, then t is the result of
+ * using the function x/2 - 1/x + 0.5 as an approximation for
+ * log2(x) for x in the range [1, 2], and scaling both the
+ * input and the result by the appropriate number of powers of 2.
*/
- return (uint32_t)(LOG10_2_x24) * (uint32_t)ilog2(x);
+ t = (ilog2x << 24) + (t >> 1) - t1 + (1 << 23);
+
+ /*
+ * Multiply t by log10(2) to get the final result.
+ *
+ * log10(2) is approximately 643/2136 We divide before
+ * multiplying to avoid overflow.
+ */
+ return t / 2136 * 643;
}
+
+#ifdef TEST_DTV_MATH
+/*
+ * To test:
+ * cc -DTEST_DTV_MATH ./dtv_math.c -lm -o ./a.out && ./a.out
+ */
+
+#include <stdio.h>
+#include <inttypes.h>
+#include <math.h>
+
+int
+main(void)
+{
+ uint32_t xlist[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
+ 14, 15, 16, 20, 24, 30, 40, 60, 100, 1000, 10000,
+ 100000, 1000000, 10000000, 100000000, 1000000000,
+ 0xffffffff};
+ int i;
+
+ printf("%11s %11s %11s %11s %s\n",
+ "x", "desired", "actual", "err_abs", "err_rel");
+ for (i = 0; i < __arraycount(xlist); i++)
+ {
+ uint32_t x = xlist[i];
+ uint32_t desired = (uint32_t)(log10((double)x)
+ * (double)(1<<24));
+ uint32_t actual = dtv_intlog10(x);
+ int32_t err_abs = actual - desired;
+ double err_rel = (err_abs == 0 ? 0.0
+ : err_abs / (double)actual);
+
+ printf("%11"PRIu32" %11"PRIu32" %11"PRIu32
+ " %+11"PRId32" %+.5f\n",
+ x, desired, actual, err_abs, err_rel);
+ }
+ return 0;
+}
+
+#endif /* TEST_DTV_MATH */
Home |
Main Index |
Thread Index |
Old Index