1 /* $NetBSD: dtv_math.c,v 1.5 2011/08/09 01:42:24 jmcneill Exp $ */
2 
3 /*-
4  * Copyright (c) 2011 Alan Barrett <apb@NetBSD.org>
5  * All rights reserved.
6  *
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions
9  * are met:
10  * 1. Redistributions of source code must retain the above copyright
11  *    notice, this list of conditions and the following disclaimer.
12  * 2. Redistributions in binary form must reproduce the above copyright
13  *    notice, this list of conditions and the following disclaimer in the
14  *    documentation and/or other materials provided with the distribution.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
17  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
18  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
19  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
20  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26  * POSSIBILITY OF SUCH DAMAGE.
27  */
28 
29 #include <sys/cdefs.h>
30 __KERNEL_RCSID(0, "$NetBSD: dtv_math.c,v 1.5 2011/08/09 01:42:24 jmcneill Exp $");
31 
32 #include <sys/types.h>
33 #include <sys/bitops.h>
34 #include <sys/module.h>
35 
36 #include <dev/dtv/dtv_math.h>
37 
38 /*
39  * dtv_intlog10 -- return an approximation to log10(x) * 1<<24,
40  * using integer arithmetic.
41  *
42  * As a special case, returns 0 when x == 0.  The mathematical
43  * result is -infinity.
44  *
45  * This function uses 0.5 + x/2 - 1/x as an approximation to
46  * log2(x) for x in the range [1.0, 2.0], and scales the input value
47  * to fit this range.  The resulting error is always better than
48  * 0.2%.
49  *
50  * Here's a table of the desired and actual results, as well
51  * as the absolute and relative errors, for several values of x.
52  *
53  *           x     desired      actual     err_abs err_rel
54  *           0           0           0          +0 +0.00000
55  *           1           0           0          +0 +0.00000
56  *           2     5050445     5050122        -323 -0.00006
57  *           3     8004766     7996348       -8418 -0.00105
58  *           4    10100890    10100887          -3 -0.00000
59  *           5    11726770    11741823      +15053 +0.00128
60  *           6    13055211    13046470       -8741 -0.00067
61  *           7    14178392    14158860      -19532 -0.00138
62  *           8    15151335    15151009        -326 -0.00002
63  *           9    16009532    16028061      +18529 +0.00116
64  *          10    16777216    16792588      +15372 +0.00092
65  *          11    17471670    17475454       +3784 +0.00022
66  *          12    18105656    18097235       -8421 -0.00047
67  *          13    18688868    18672077      -16791 -0.00090
68  *          14    19228837    19209625      -19212 -0.00100
69  *          15    19731537    19717595      -13942 -0.00071
70  *          16    20201781    20201774          -7 -0.00000
71  *          20    21827661    21842710      +15049 +0.00069
72  *          24    23156102    23147357       -8745 -0.00038
73  *          30    24781982    24767717      -14265 -0.00058
74  *          40    26878106    26893475      +15369 +0.00057
75  *          60    29832427    29818482      -13945 -0.00047
76  *         100    33554432    33540809      -13623 -0.00041
77  *        1000    50331648    50325038       -6610 -0.00013
78  *       10000    67108864    67125985      +17121 +0.00026
79  *      100000    83886080    83875492      -10588 -0.00013
80  *     1000000   100663296   100652005      -11291 -0.00011
81  *    10000000   117440512   117458739      +18227 +0.00016
82  *   100000000   134217728   134210175       -7553 -0.00006
83  *  1000000000   150994944   150980258      -14686 -0.00010
84  *  4294967295   161614248   161614192         -56 -0.00000
85  */
86 uint32_t
dtv_intlog10(uint32_t x)87 dtv_intlog10(uint32_t x)
88 {
89 	uint32_t ilog2x;
90 	uint32_t t;
91 	uint32_t t1;
92 
93 	if (__predict_false(x == 0))
94 		return 0;
95 
96 	/*
97 	 * find ilog2x = floor(log2(x)), as an integer in the range [0,31].
98 	 */
99 	ilog2x = ilog2(x);
100 
101 	/*
102 	 * Set "t" to the result of shifting x left or right
103 	 * until the most significant bit that was actually set
104 	 * moves into the 1<<24 position.
105 	 *
106 	 * Now we can think of "t" as representing
107 	 * x / 2**(floor(log2(x))),
108 	 * as a fixed-point value with 8 integer bits and 24 fraction bits.
109 	 *
110 	 * This value is in the semi-closed interval [1.0, 2.0)
111 	 * when interpreting it as a fixed-point number, or in the
112 	 * interval [0x01000000, 0x01ffffff] when examining the
113 	 * underlying uint32_t representation.
114 	 */
115 	t = (ilog2x > 24 ? x >> (ilog2x - 24) : x << (24 - ilog2x));
116 
117 	/*
118 	 * Calculate "t1 = 1 / t" in the 8.24 fixed-point format.
119 	 * This value is in the interval [0.5, 1.0]
120 	 * when interpreting it as a fixed-point number, or in the
121 	 * interval [0x00800000, 0x01000000] when examining the
122 	 * underlying uint32_t representation.
123 	 *
124 	 */
125 	t1 = ((uint64_t)1 << 48) / t;
126 
127 	/*
128 	 * Calculate "t = ilog2x + t/2 - t1 + 0.5" in the 8.24
129 	 * fixed-point format.
130 	 *
131 	 * If x is a power of 2, then t is now exactly equal to log2(x)
132 	 * when interpreting it as a fixed-point number, or exactly
133 	 * log2(x) << 24 when examining the underlying uint32_t
134 	 * representation.
135 	 *
136 	 * If x is not a power of 2, then t is the result of
137 	 * using the function x/2 - 1/x + 0.5 as an approximation for
138 	 * log2(x) for x in the range [1, 2], and scaling both the
139 	 * input and the result by the appropriate number of powers of 2.
140 	 */
141 	t = (ilog2x << 24) + (t >> 1) - t1 + (1 << 23);
142 
143 	/*
144 	 * Multiply t by log10(2) to get the final result.
145 	 *
146 	 * log10(2) is approximately 643/2136  We divide before
147 	 * multiplying to avoid overflow.
148 	 */
149 	return t / 2136 * 643;
150 }
151 
152 #ifdef _KERNEL
153 MODULE(MODULE_CLASS_MISC, dtv_math, NULL);
154 
155 static int
dtv_math_modcmd(modcmd_t cmd,void * opaque)156 dtv_math_modcmd(modcmd_t cmd, void *opaque)
157 {
158 	if (cmd == MODULE_CMD_INIT || cmd == MODULE_CMD_FINI)
159 		return 0;
160 	return ENOTTY;
161 }
162 #endif
163 
164 #ifdef TEST_DTV_MATH
165 /*
166  * To test:
167  *	cc -DTEST_DTV_MATH ./dtv_math.c -lm -o ./a.out && ./a.out
168  */
169 
170 #include <stdio.h>
171 #include <inttypes.h>
172 #include <math.h>
173 
174 int
main(void)175 main(void)
176 {
177 	uint32_t xlist[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
178 			14, 15, 16, 20, 24, 30, 40, 60, 100, 1000, 10000,
179 			100000, 1000000, 10000000, 100000000, 1000000000,
180 			0xffffffff};
181 	int i;
182 
183 	printf("%11s %11s %11s %11s %s\n",
184 		"x", "desired", "actual", "err_abs", "err_rel");
185 	for (i = 0; i < __arraycount(xlist); i++)
186 	{
187 		uint32_t x = xlist[i];
188 		uint32_t desired = (uint32_t)(log10((double)x)
189 						* (double)(1<<24));
190 		uint32_t actual = dtv_intlog10(x);
191 		int32_t err_abs = actual - desired;
192 		double err_rel = (err_abs == 0 ? 0.0
193 				: err_abs / (double)actual);
194 
195 		printf("%11"PRIu32" %11"PRIu32" %11"PRIu32
196 			" %+11"PRId32" %+.5f\n",
197 			x, desired, actual, err_abs, err_rel);
198 	}
199 	return 0;
200 }
201 
202 #endif /* TEST_DTV_MATH */
203