summaryrefslogtreecommitdiff
path: root/sys/dev/dtv/dtv_math.c
blob: 45a715f30b16d713f55909056ccfa8ac196efc30 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
/* $NetBSD: dtv_math.c,v 1.5 2011/08/09 01:42:24 jmcneill Exp $ */

/*-
 * Copyright (c) 2011 Alan Barrett <apb@NetBSD.org>
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include <sys/cdefs.h>
__KERNEL_RCSID(0, "$NetBSD: dtv_math.c,v 1.5 2011/08/09 01:42:24 jmcneill Exp $");

#include <sys/types.h>
#include <sys/bitops.h>
#include <sys/module.h>

#include <dev/dtv/dtv_math.h>

/*
 * dtv_intlog10 -- return an approximation to log10(x) * 1<<24,
 * using integer arithmetic.
 *
 * As a special case, returns 0 when x == 0.  The mathematical
 * result is -infinity.
 *
 * 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.
 *
 *           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;

	/*
	 * 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.
	 */
	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 _KERNEL
MODULE(MODULE_CLASS_MISC, dtv_math, NULL);

static int
dtv_math_modcmd(modcmd_t cmd, void *opaque)
{
	if (cmd == MODULE_CMD_INIT || cmd == MODULE_CMD_FINI)
		return 0;
	return ENOTTY;
}
#endif

#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 */