sys/arch/m68k/fpe/fpu_log.c

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