1/* $NetBSD: dtv_math.c,v 1.4 2011/07/16 22:41:59 apb 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.4 2011/07/16 22:41:59 apb 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 */
86uint32_t
87dtv_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
153MODULE(MODULE_CLASS_MISC, dtv_math, NULL);
154
155static int
156dtv_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
174int
175main(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