fpu_log.c revision 1.6
1 /* $NetBSD: fpu_log.c,v 1.6 1999/05/30 20:17:48 briggs Exp $ */ 2 3 /* 4 * Copyright (c) 1995 Ken Nakata 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 * 3. Neither the name of the author nor the names of its contributors 16 * may be used to endorse or promote products derived from this software 17 * without specific prior written permission. 18 * 19 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 29 * SUCH DAMAGE. 30 * 31 * @(#)fpu_log.c 10/8/95 32 */ 33 34 #include <sys/types.h> 35 #include <sys/systm.h> 36 37 #include "fpu_emulate.h" 38 39 static u_int logA6[] = { 0x3FC2499A, 0xB5E4040B }; 40 static u_int logA5[] = { 0xBFC555B5, 0x848CB7DB }; 41 static u_int logA4[] = { 0x3FC99999, 0x987D8730 }; 42 static u_int logA3[] = { 0xBFCFFFFF, 0xFF6F7E97 }; 43 static u_int logA2[] = { 0x3FD55555, 0x555555A4 }; 44 static u_int logA1[] = { 0xBFE00000, 0x00000008 }; 45 46 static u_int logB5[] = { 0x3F175496, 0xADD7DAD6 }; 47 static u_int logB4[] = { 0x3F3C71C2, 0xFE80C7E0 }; 48 static u_int logB3[] = { 0x3F624924, 0x928BCCFF }; 49 static u_int logB2[] = { 0x3F899999, 0x999995EC }; 50 static u_int logB1[] = { 0x3FB55555, 0x55555555 }; 51 52 /* sfpn = shortened fp number; can represent only positive numbers */ 53 static struct sfpn { 54 int sp_exp; 55 u_int sp_m0, sp_m1; 56 } logtbl[] = { 57 { 0x3FFE - 0x3fff, 0xFE03F80FU, 0xE03F80FEU }, 58 { 0x3FF7 - 0x3fff, 0xFF015358U, 0x833C47E2U }, 59 { 0x3FFE - 0x3fff, 0xFA232CF2U, 0x52138AC0U }, 60 { 0x3FF9 - 0x3fff, 0xBDC8D83EU, 0xAD88D549U }, 61 { 0x3FFE - 0x3fff, 0xF6603D98U, 0x0F6603DAU }, 62 { 0x3FFA - 0x3fff, 0x9CF43DCFU, 0xF5EAFD48U }, 63 { 0x3FFE - 0x3fff, 0xF2B9D648U, 0x0F2B9D65U }, 64 { 0x3FFA - 0x3fff, 0xDA16EB88U, 0xCB8DF614U }, 65 { 0x3FFE - 0x3fff, 0xEF2EB71FU, 0xC4345238U }, 66 { 0x3FFB - 0x3fff, 0x8B29B775U, 0x1BD70743U }, 67 { 0x3FFE - 0x3fff, 0xEBBDB2A5U, 0xC1619C8CU }, 68 { 0x3FFB - 0x3fff, 0xA8D839F8U, 0x30C1FB49U }, 69 { 0x3FFE - 0x3fff, 0xE865AC7BU, 0x7603A197U }, 70 { 0x3FFB - 0x3fff, 0xC61A2EB1U, 0x8CD907ADU }, 71 { 0x3FFE - 0x3fff, 0xE525982AU, 0xF70C880EU }, 72 { 0x3FFB - 0x3fff, 0xE2F2A47AU, 0xDE3A18AFU }, 73 { 0x3FFE - 0x3fff, 0xE1FC780EU, 0x1FC780E2U }, 74 { 0x3FFB - 0x3fff, 0xFF64898EU, 0xDF55D551U }, 75 { 0x3FFE - 0x3fff, 0xDEE95C4CU, 0xA037BA57U }, 76 { 0x3FFC - 0x3fff, 0x8DB956A9U, 0x7B3D0148U }, 77 { 0x3FFE - 0x3fff, 0xDBEB61EEU, 0xD19C5958U }, 78 { 0x3FFC - 0x3fff, 0x9B8FE100U, 0xF47BA1DEU }, 79 { 0x3FFE - 0x3fff, 0xD901B203U, 0x6406C80EU }, 80 { 0x3FFC - 0x3fff, 0xA9372F1DU, 0x0DA1BD17U }, 81 { 0x3FFE - 0x3fff, 0xD62B80D6U, 0x2B80D62CU }, 82 { 0x3FFC - 0x3fff, 0xB6B07F38U, 0xCE90E46BU }, 83 { 0x3FFE - 0x3fff, 0xD3680D36U, 0x80D3680DU }, 84 { 0x3FFC - 0x3fff, 0xC3FD0329U, 0x06488481U }, 85 { 0x3FFE - 0x3fff, 0xD0B69FCBU, 0xD2580D0BU }, 86 { 0x3FFC - 0x3fff, 0xD11DE0FFU, 0x15AB18CAU }, 87 { 0x3FFE - 0x3fff, 0xCE168A77U, 0x25080CE1U }, 88 { 0x3FFC - 0x3fff, 0xDE1433A1U, 0x6C66B150U }, 89 { 0x3FFE - 0x3fff, 0xCB8727C0U, 0x65C393E0U }, 90 { 0x3FFC - 0x3fff, 0xEAE10B5AU, 0x7DDC8ADDU }, 91 { 0x3FFE - 0x3fff, 0xC907DA4EU, 0x871146ADU }, 92 { 0x3FFC - 0x3fff, 0xF7856E5EU, 0xE2C9B291U }, 93 { 0x3FFE - 0x3fff, 0xC6980C69U, 0x80C6980CU }, 94 { 0x3FFD - 0x3fff, 0x82012CA5U, 0xA68206D7U }, 95 { 0x3FFE - 0x3fff, 0xC4372F85U, 0x5D824CA6U }, 96 { 0x3FFD - 0x3fff, 0x882C5FCDU, 0x7256A8C5U }, 97 { 0x3FFE - 0x3fff, 0xC1E4BBD5U, 0x95F6E947U }, 98 { 0x3FFD - 0x3fff, 0x8E44C60BU, 0x4CCFD7DEU }, 99 { 0x3FFE - 0x3fff, 0xBFA02FE8U, 0x0BFA02FFU }, 100 { 0x3FFD - 0x3fff, 0x944AD09EU, 0xF4351AF6U }, 101 { 0x3FFE - 0x3fff, 0xBD691047U, 0x07661AA3U }, 102 { 0x3FFD - 0x3fff, 0x9A3EECD4U, 0xC3EAA6B2U }, 103 { 0x3FFE - 0x3fff, 0xBB3EE721U, 0xA54D880CU }, 104 { 0x3FFD - 0x3fff, 0xA0218434U, 0x353F1DE8U }, 105 { 0x3FFE - 0x3fff, 0xB92143FAU, 0x36F5E02EU }, 106 { 0x3FFD - 0x3fff, 0xA5F2FCABU, 0xBBC506DAU }, 107 { 0x3FFE - 0x3fff, 0xB70FBB5AU, 0x19BE3659U }, 108 { 0x3FFD - 0x3fff, 0xABB3B8BAU, 0x2AD362A5U }, 109 { 0x3FFE - 0x3fff, 0xB509E68AU, 0x9B94821FU }, 110 { 0x3FFD - 0x3fff, 0xB1641795U, 0xCE3CA97BU }, 111 { 0x3FFE - 0x3fff, 0xB30F6352U, 0x8917C80BU }, 112 { 0x3FFD - 0x3fff, 0xB7047551U, 0x5D0F1C61U }, 113 { 0x3FFE - 0x3fff, 0xB11FD3B8U, 0x0B11FD3CU }, 114 { 0x3FFD - 0x3fff, 0xBC952AFEU, 0xEA3D13E1U }, 115 { 0x3FFE - 0x3fff, 0xAF3ADDC6U, 0x80AF3ADEU }, 116 { 0x3FFD - 0x3fff, 0xC2168ED0U, 0xF458BA4AU }, 117 { 0x3FFE - 0x3fff, 0xAD602B58U, 0x0AD602B6U }, 118 { 0x3FFD - 0x3fff, 0xC788F439U, 0xB3163BF1U }, 119 { 0x3FFE - 0x3fff, 0xAB8F69E2U, 0x8359CD11U }, 120 { 0x3FFD - 0x3fff, 0xCCECAC08U, 0xBF04565DU }, 121 { 0x3FFE - 0x3fff, 0xA9C84A47U, 0xA07F5638U }, 122 { 0x3FFD - 0x3fff, 0xD2420487U, 0x2DD85160U }, 123 { 0x3FFE - 0x3fff, 0xA80A80A8U, 0x0A80A80BU }, 124 { 0x3FFD - 0x3fff, 0xD7894992U, 0x3BC3588AU }, 125 { 0x3FFE - 0x3fff, 0xA655C439U, 0x2D7B73A8U }, 126 { 0x3FFD - 0x3fff, 0xDCC2C4B4U, 0x9887DACCU }, 127 { 0x3FFE - 0x3fff, 0xA4A9CF1DU, 0x96833751U }, 128 { 0x3FFD - 0x3fff, 0xE1EEBD3EU, 0x6D6A6B9EU }, 129 { 0x3FFE - 0x3fff, 0xA3065E3FU, 0xAE7CD0E0U }, 130 { 0x3FFD - 0x3fff, 0xE70D785CU, 0x2F9F5BDCU }, 131 { 0x3FFE - 0x3fff, 0xA16B312EU, 0xA8FC377DU }, 132 { 0x3FFD - 0x3fff, 0xEC1F392CU, 0x5179F283U }, 133 { 0x3FFE - 0x3fff, 0x9FD809FDU, 0x809FD80AU }, 134 { 0x3FFD - 0x3fff, 0xF12440D3U, 0xE36130E6U }, 135 { 0x3FFE - 0x3fff, 0x9E4CAD23U, 0xDD5F3A20U }, 136 { 0x3FFD - 0x3fff, 0xF61CCE92U, 0x346600BBU }, 137 { 0x3FFE - 0x3fff, 0x9CC8E160U, 0xC3FB19B9U }, 138 { 0x3FFD - 0x3fff, 0xFB091FD3U, 0x8145630AU }, 139 { 0x3FFE - 0x3fff, 0x9B4C6F9EU, 0xF03A3CAAU }, 140 { 0x3FFD - 0x3fff, 0xFFE97042U, 0xBFA4C2ADU }, 141 { 0x3FFE - 0x3fff, 0x99D722DAU, 0xBDE58F06U }, 142 { 0x3FFE - 0x3fff, 0x825EFCEDU, 0x49369330U }, 143 { 0x3FFE - 0x3fff, 0x9868C809U, 0x868C8098U }, 144 { 0x3FFE - 0x3fff, 0x84C37A7AU, 0xB9A905C9U }, 145 { 0x3FFE - 0x3fff, 0x97012E02U, 0x5C04B809U }, 146 { 0x3FFE - 0x3fff, 0x87224C2EU, 0x8E645FB7U }, 147 { 0x3FFE - 0x3fff, 0x95A02568U, 0x095A0257U }, 148 { 0x3FFE - 0x3fff, 0x897B8CACU, 0x9F7DE298U }, 149 { 0x3FFE - 0x3fff, 0x94458094U, 0x45809446U }, 150 { 0x3FFE - 0x3fff, 0x8BCF55DEU, 0xC4CD05FEU }, 151 { 0x3FFE - 0x3fff, 0x92F11384U, 0x0497889CU }, 152 { 0x3FFE - 0x3fff, 0x8E1DC0FBU, 0x89E125E5U }, 153 { 0x3FFE - 0x3fff, 0x91A2B3C4U, 0xD5E6F809U }, 154 { 0x3FFE - 0x3fff, 0x9066E68CU, 0x955B6C9BU }, 155 { 0x3FFE - 0x3fff, 0x905A3863U, 0x3E06C43BU }, 156 { 0x3FFE - 0x3fff, 0x92AADE74U, 0xC7BE59E0U }, 157 { 0x3FFE - 0x3fff, 0x8F1779D9U, 0xFDC3A219U }, 158 { 0x3FFE - 0x3fff, 0x94E9BFF6U, 0x15845643U }, 159 { 0x3FFE - 0x3fff, 0x8DDA5202U, 0x37694809U }, 160 { 0x3FFE - 0x3fff, 0x9723A1B7U, 0x20134203U }, 161 { 0x3FFE - 0x3fff, 0x8CA29C04U, 0x6514E023U }, 162 { 0x3FFE - 0x3fff, 0x995899C8U, 0x90EB8990U }, 163 { 0x3FFE - 0x3fff, 0x8B70344AU, 0x139BC75AU }, 164 { 0x3FFE - 0x3fff, 0x9B88BDAAU, 0x3A3DAE2FU }, 165 { 0x3FFE - 0x3fff, 0x8A42F870U, 0x5669DB46U }, 166 { 0x3FFE - 0x3fff, 0x9DB4224FU, 0xFFE1157CU }, 167 { 0x3FFE - 0x3fff, 0x891AC73AU, 0xE9819B50U }, 168 { 0x3FFE - 0x3fff, 0x9FDADC26U, 0x8B7A12DAU }, 169 { 0x3FFE - 0x3fff, 0x87F78087U, 0xF78087F8U }, 170 { 0x3FFE - 0x3fff, 0xA1FCFF17U, 0xCE733BD4U }, 171 { 0x3FFE - 0x3fff, 0x86D90544U, 0x7A34ACC6U }, 172 { 0x3FFE - 0x3fff, 0xA41A9E8FU, 0x5446FB9FU }, 173 { 0x3FFE - 0x3fff, 0x85BF3761U, 0x2CEE3C9BU }, 174 { 0x3FFE - 0x3fff, 0xA633CD7EU, 0x6771CD8BU }, 175 { 0x3FFE - 0x3fff, 0x84A9F9C8U, 0x084A9F9DU }, 176 { 0x3FFE - 0x3fff, 0xA8489E60U, 0x0B435A5EU }, 177 { 0x3FFE - 0x3fff, 0x83993052U, 0x3FBE3368U }, 178 { 0x3FFE - 0x3fff, 0xAA59233CU, 0xCCA4BD49U }, 179 { 0x3FFE - 0x3fff, 0x828CBFBEU, 0xB9A020A3U }, 180 { 0x3FFE - 0x3fff, 0xAC656DAEU, 0x6BCC4985U }, 181 { 0x3FFE - 0x3fff, 0x81848DA8U, 0xFAF0D277U }, 182 { 0x3FFE - 0x3fff, 0xAE6D8EE3U, 0x60BB2468U }, 183 { 0x3FFE - 0x3fff, 0x80808080U, 0x80808081U }, 184 { 0x3FFE - 0x3fff, 0xB07197A2U, 0x3C46C654U }, 185 }; 186 187 static struct fpn *__fpu_logn __P((struct fpemu *fe)); 188 189 /* 190 * natural log - algorithm taken from Motorola FPSP, 191 * except this doesn't bother to check for invalid input. 192 */ 193 static struct fpn * 194 __fpu_logn(fe) 195 struct fpemu *fe; 196 { 197 static struct fpn X, F, U, V, W, KLOG2; 198 struct fpn *d; 199 int i, k; 200 201 CPYFPN(&X, &fe->fe_f2); 202 203 /* see if |X-1| < 1/16 approx. */ 204 if ((-1 == X.fp_exp && (0xf07d0000U >> (31 - FP_LG)) <= X.fp_mant[0]) || 205 (0 == X.fp_exp && X.fp_mant[0] <= (0x88410000U >> (31 - FP_LG)))) { 206 /* log near 1 */ 207 #if FPE_DEBUG 208 printf("__fpu_logn: log near 1\n"); 209 #endif 210 211 fpu_const(&fe->fe_f1, 0x32); 212 /* X+1 */ 213 d = fpu_add(fe); 214 CPYFPN(&V, d); 215 216 CPYFPN(&fe->fe_f1, &X); 217 fpu_const(&fe->fe_f2, 0x32); /* 1.0 */ 218 fe->fe_f2.fp_sign = 1; /* -1.0 */ 219 /* X-1 */ 220 d = fpu_add(fe); 221 CPYFPN(&fe->fe_f1, d); 222 /* 2(X-1) */ 223 fe->fe_f1.fp_exp++; /* *= 2 */ 224 CPYFPN(&fe->fe_f2, &V); 225 /* U=2(X-1)/(X+1) */ 226 d = fpu_div(fe); 227 CPYFPN(&U, d); 228 CPYFPN(&fe->fe_f1, d); 229 CPYFPN(&fe->fe_f2, d); 230 /* V=U*U */ 231 d = fpu_mul(fe); 232 CPYFPN(&V, d); 233 CPYFPN(&fe->fe_f1, d); 234 CPYFPN(&fe->fe_f2, d); 235 /* W=V*V */ 236 d = fpu_mul(fe); 237 CPYFPN(&W, d); 238 239 /* calculate U+U*V*([B1+W*(B3+W*B5)]+[V*(B2+W*B4)]) */ 240 241 /* B1+W*(B3+W*B5) part */ 242 CPYFPN(&fe->fe_f1, d); 243 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB5); 244 /* W*B5 */ 245 d = fpu_mul(fe); 246 CPYFPN(&fe->fe_f1, d); 247 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB3); 248 /* B3+W*B5 */ 249 d = fpu_add(fe); 250 CPYFPN(&fe->fe_f1, d); 251 CPYFPN(&fe->fe_f2, &W); 252 /* W*(B3+W*B5) */ 253 d = fpu_mul(fe); 254 CPYFPN(&fe->fe_f1, d); 255 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB1); 256 /* B1+W*(B3+W*B5) */ 257 d = fpu_add(fe); 258 CPYFPN(&X, d); 259 260 /* [V*(B2+W*B4)] part */ 261 CPYFPN(&fe->fe_f1, &W); 262 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB4); 263 /* W*B4 */ 264 d = fpu_mul(fe); 265 CPYFPN(&fe->fe_f1, d); 266 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB2); 267 /* B2+W*B4 */ 268 d = fpu_add(fe); 269 CPYFPN(&fe->fe_f1, d); 270 CPYFPN(&fe->fe_f2, &V); 271 /* V*(B2+W*B4) */ 272 d = fpu_mul(fe); 273 CPYFPN(&fe->fe_f1, d); 274 CPYFPN(&fe->fe_f2, &X); 275 /* B1+W*(B3+W*B5)+V*(B2+W*B4) */ 276 d = fpu_add(fe); 277 CPYFPN(&fe->fe_f1, d); 278 CPYFPN(&fe->fe_f2, &V); 279 /* V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */ 280 d = fpu_mul(fe); 281 CPYFPN(&fe->fe_f1, d); 282 CPYFPN(&fe->fe_f2, &U); 283 /* U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */ 284 d = fpu_mul(fe); 285 CPYFPN(&fe->fe_f1, d); 286 CPYFPN(&fe->fe_f2, &U); 287 /* U+U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */ 288 d = fpu_add(fe); 289 } else /* the usual case */ { 290 #if FPE_DEBUG 291 printf("__fpu_logn: the usual case. X=(%d,%08x,%08x...)\n", 292 X.fp_exp, X.fp_mant[0], X.fp_mant[1]); 293 #endif 294 295 k = X.fp_exp; 296 /* X <- Y */ 297 X.fp_exp = fe->fe_f2.fp_exp = 0; 298 299 /* get the most significant 7 bits of X */ 300 F.fp_class = FPC_NUM; 301 F.fp_sign = 0; 302 F.fp_exp = X.fp_exp; 303 F.fp_mant[0] = X.fp_mant[0] & (0xfe000000U >> (31 - FP_LG)); 304 F.fp_mant[0] |= (0x01000000U >> (31 - FP_LG)); 305 F.fp_mant[1] = F.fp_mant[2] = 0; 306 F.fp_sticky = 0; 307 308 #if FPE_DEBUG 309 printf("__fpu_logn: X=Y*2^k=(%d,%08x,%08x...)*2^%d\n", 310 fe->fe_f2.fp_exp, fe->fe_f2.fp_mant[0], 311 fe->fe_f2.fp_mant[1], k); 312 printf("__fpu_logn: F=(%d,%08x,%08x...)\n", 313 F.fp_exp, F.fp_mant[0], F.fp_mant[1]); 314 #endif 315 316 /* index to the table */ 317 i = ((F.fp_mant[0] << (7 - FP_LG)) | 318 (F.fp_mant[1] >> (32 - (7 - FP_LG)))) & 0x7e; 319 320 #if FPE_DEBUG 321 printf("__fpu_logn: index to logtbl i=%d(%x)\n", i, i); 322 #endif 323 324 CPYFPN(&fe->fe_f1, &F); 325 /* -F */ 326 fe->fe_f1.fp_sign = 1; 327 /* Y-F */ 328 d = fpu_add(fe); 329 CPYFPN(&fe->fe_f1, d); 330 331 /* fe_f2 = 1/F */ 332 fe->fe_f2.fp_class = FPC_NUM; 333 fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[2] = 0; 334 fe->fe_f2.fp_exp = logtbl[i].sp_exp; 335 fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG)); 336 fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) | 337 (logtbl[i].sp_m1 >> (31 - FP_LG)); 338 fe->fe_f2.fp_mant[2] = (u_int)(logtbl[i].sp_m1 << (FP_LG + 1)); 339 340 #if FPE_DEBUG 341 printf("__fpu_logn: 1/F=(%d,%08x,%08x...)\n", fe->fe_f2.fp_exp, 342 fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]); 343 #endif 344 345 /* U = (Y-F) * (1/F) */ 346 d = fpu_mul(fe); 347 CPYFPN(&U, d); 348 349 /* KLOG2 = K * ln(2) */ 350 /* fe_f1 == (fpn)k */ 351 fpu_explode(fe, &fe->fe_f1, FTYPE_LNG, &k); 352 (void)fpu_const(&fe->fe_f2, 0x30 /* ln(2) */); 353 #if FPE_DEBUG 354 printf("__fpu_logn: fp(k)=(%d,%08x,%08x...)\n", fe->fe_f1.fp_exp, 355 fe->fe_f1.fp_mant[0], fe->fe_f1.fp_mant[1]); 356 printf("__fpu_logn: ln(2)=(%d,%08x,%08x...)\n", fe->fe_f2.fp_exp, 357 fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]); 358 #endif 359 /* K * LOGOF2 */ 360 d = fpu_mul(fe); 361 CPYFPN(&KLOG2, d); 362 363 /* V=U*U */ 364 CPYFPN(&fe->fe_f1, &U); 365 CPYFPN(&fe->fe_f2, &U); 366 d = fpu_mul(fe); 367 CPYFPN(&V, d); 368 369 /* 370 * approximation of LOG(1+U) by 371 * (U+V*(A1+V*(A3+V*A5)))+(U*V*(A2+V*(A4+V*A6))) 372 */ 373 374 /* (U+V*(A1+V*(A3+V*A5))) part */ 375 CPYFPN(&fe->fe_f1, d); 376 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA5); 377 /* V*A5 */ 378 d = fpu_mul(fe); 379 380 CPYFPN(&fe->fe_f1, d); 381 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA3); 382 /* A3+V*A5 */ 383 d = fpu_add(fe); 384 385 CPYFPN(&fe->fe_f1, d); 386 CPYFPN(&fe->fe_f2, &V); 387 /* V*(A3+V*A5) */ 388 d = fpu_mul(fe); 389 390 CPYFPN(&fe->fe_f1, d); 391 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA1); 392 /* A1+V*(A3+V*A5) */ 393 d = fpu_add(fe); 394 395 CPYFPN(&fe->fe_f1, d); 396 CPYFPN(&fe->fe_f2, &V); 397 /* V*(A1+V*(A3+V*A5)) */ 398 d = fpu_mul(fe); 399 400 CPYFPN(&fe->fe_f1, d); 401 CPYFPN(&fe->fe_f2, &U); 402 /* U+V*(A1+V*(A3+V*A5)) */ 403 d = fpu_add(fe); 404 405 CPYFPN(&X, d); 406 407 /* (U*V*(A2+V*(A4+V*A6))) part */ 408 CPYFPN(&fe->fe_f1, &V); 409 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA6); 410 /* V*A6 */ 411 d = fpu_mul(fe); 412 CPYFPN(&fe->fe_f1, d); 413 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA4); 414 /* A4+V*A6 */ 415 d = fpu_add(fe); 416 CPYFPN(&fe->fe_f1, d); 417 CPYFPN(&fe->fe_f2, &V); 418 /* V*(A4+V*A6) */ 419 d = fpu_mul(fe); 420 CPYFPN(&fe->fe_f1, d); 421 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA2); 422 /* A2+V*(A4+V*A6) */ 423 d = fpu_add(fe); 424 CPYFPN(&fe->fe_f1, d); 425 CPYFPN(&fe->fe_f2, &V); 426 /* V*(A2+V*(A4+V*A6)) */ 427 d = fpu_mul(fe); 428 CPYFPN(&fe->fe_f1, d); 429 CPYFPN(&fe->fe_f2, &U); 430 /* U*V*(A2+V*(A4+V*A6)) */ 431 d = fpu_mul(fe); 432 CPYFPN(&fe->fe_f1, d); 433 i++; 434 /* fe_f2 = logtbl[i+1] (== LOG(F)) */ 435 fe->fe_f2.fp_class = FPC_NUM; 436 fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[2] = 0; 437 fe->fe_f2.fp_exp = logtbl[i].sp_exp; 438 fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG)); 439 fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) | 440 (logtbl[i].sp_m1 >> (31 - FP_LG)); 441 fe->fe_f2.fp_mant[2] = (logtbl[i].sp_m1 << (FP_LG + 1)); 442 443 #if FPE_DEBUG 444 printf("__fpu_logn: ln(F)=(%d,%08x,%08x,...)\n", fe->fe_f2.fp_exp, 445 fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]); 446 #endif 447 448 /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */ 449 d = fpu_add(fe); 450 CPYFPN(&fe->fe_f1, d); 451 CPYFPN(&fe->fe_f2, &X); 452 /* LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */ 453 d = fpu_add(fe); 454 455 #if FPE_DEBUG 456 printf("__fpu_logn: ln(Y)=(%c,%d,%08x,%08x,%08x)\n", 457 d->fp_sign ? '-' : '+', d->fp_exp, 458 d->fp_mant[0], d->fp_mant[1], d->fp_mant[2]); 459 #endif 460 461 CPYFPN(&fe->fe_f1, d); 462 CPYFPN(&fe->fe_f2, &KLOG2); 463 /* K*LOGOF2+LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */ 464 d = fpu_add(fe); 465 } 466 467 return d; 468 } 469 470 struct fpn * 471 fpu_log10(fe) 472 struct fpemu *fe; 473 { 474 struct fpn *fp = &fe->fe_f2; 475 u_int fpsr; 476 477 fpsr = fe->fe_fpsr & ~FPSR_EXCP; /* clear all exceptions */ 478 479 if (fp->fp_class >= FPC_NUM) { 480 if (fp->fp_sign) { /* negative number or Inf */ 481 fp = fpu_newnan(fe); 482 fpsr |= FPSR_OPERR; 483 } else if (fp->fp_class == FPC_NUM) { 484 /* the real work here */ 485 fp = __fpu_logn(fe); 486 if (fp != &fe->fe_f1) 487 CPYFPN(&fe->fe_f1, fp); 488 (void)fpu_const(&fe->fe_f2, 0x31 /* ln(10) */); 489 fp = fpu_div(fe); 490 } /* else if fp == +Inf, return +Inf */ 491 } else if (fp->fp_class == FPC_ZERO) { 492 /* return -Inf */ 493 fp->fp_class = FPC_INF; 494 fp->fp_sign = 1; 495 fpsr |= FPSR_DZ; 496 } else if (fp->fp_class == FPC_SNAN) { 497 fpsr |= FPSR_SNAN; 498 fp = fpu_newnan(fe); 499 } else { 500 fp = fpu_newnan(fe); 501 } 502 503 fe->fe_fpsr = fpsr; 504 505 return fp; 506 } 507 508 struct fpn * 509 fpu_log2(fe) 510 struct fpemu *fe; 511 { 512 struct fpn *fp = &fe->fe_f2; 513 u_int fpsr; 514 515 fpsr = fe->fe_fpsr & ~FPSR_EXCP; /* clear all exceptions */ 516 517 if (fp->fp_class >= FPC_NUM) { 518 if (fp->fp_sign) { /* negative number or Inf */ 519 fp = fpu_newnan(fe); 520 fpsr |= FPSR_OPERR; 521 } else if (fp->fp_class == FPC_NUM) { 522 /* the real work here */ 523 if (fp->fp_mant[0] == FP_1 && fp->fp_mant[1] == 0 && 524 fp->fp_mant[2] == 0) { 525 /* fp == 2.0 ^ exp <--> log2(fp) == exp */ 526 fpu_explode(fe, &fe->fe_f3, FTYPE_LNG, &fp->fp_exp); 527 fp = &fe->fe_f3; 528 } else { 529 fp = __fpu_logn(fe); 530 if (fp != &fe->fe_f1) 531 CPYFPN(&fe->fe_f1, fp); 532 (void)fpu_const(&fe->fe_f2, 0x30 /* ln(2) */); 533 fp = fpu_div(fe); 534 } 535 } /* else if fp == +Inf, return +Inf */ 536 } else if (fp->fp_class == FPC_ZERO) { 537 /* return -Inf */ 538 fp->fp_class = FPC_INF; 539 fp->fp_sign = 1; 540 fpsr |= FPSR_DZ; 541 } else if (fp->fp_class == FPC_SNAN) { 542 fpsr |= FPSR_SNAN; 543 fp = fpu_newnan(fe); 544 } else { 545 fp = fpu_newnan(fe); 546 } 547 548 fe->fe_fpsr = fpsr; 549 return fp; 550 } 551 552 struct fpn * 553 fpu_logn(fe) 554 struct fpemu *fe; 555 { 556 struct fpn *fp = &fe->fe_f2; 557 u_int fpsr; 558 559 fpsr = fe->fe_fpsr & ~FPSR_EXCP; /* clear all exceptions */ 560 561 if (fp->fp_class >= FPC_NUM) { 562 if (fp->fp_sign) { /* negative number or Inf */ 563 fp = fpu_newnan(fe); 564 fpsr |= FPSR_OPERR; 565 } else if (fp->fp_class == FPC_NUM) { 566 /* the real work here */ 567 fp = __fpu_logn(fe); 568 } /* else if fp == +Inf, return +Inf */ 569 } else if (fp->fp_class == FPC_ZERO) { 570 /* return -Inf */ 571 fp->fp_class = FPC_INF; 572 fp->fp_sign = 1; 573 fpsr |= FPSR_DZ; 574 } else if (fp->fp_class == FPC_SNAN) { 575 fpsr |= FPSR_SNAN; 576 fp = fpu_newnan(fe); 577 } else { 578 fp = fpu_newnan(fe); 579 } 580 581 fe->fe_fpsr = fpsr; 582 583 return fp; 584 } 585 586 struct fpn * 587 fpu_lognp1(fe) 588 struct fpemu *fe; 589 { 590 struct fpn *fp; 591 592 /* build a 1.0 */ 593 fp = fpu_const(&fe->fe_f1, 0x32); /* get 1.0 */ 594 /* fp = 1.0 + f2 */ 595 fp = fpu_add(fe); 596 597 /* copy the result to the src opr */ 598 if (&fe->fe_f2 != fp) 599 CPYFPN(&fe->fe_f2, fp); 600 601 return fpu_logn(fe); 602 } 603
Indexes created Sun Sep 28 07:09:56 GMT 2025