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 (at) 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 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 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 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
Indexes created Tue Sep 30 11:09:46 GMT 2025