sys/platform/pc32/gnu/fpemul/poly_l2.c

   1 /*
   2  *  poly_l2.c
   3  *
   4  * Compute the base 2 log of a FPU_REG, using a polynomial approximation.
   5  *
   6  *
   7  * Copyright (C) 1992,1993,1994
   8  *                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
   9  *                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au
  10  * All rights reserved.
  11  *
  12  * This copyright notice covers the redistribution and use of the
  13  * FPU emulator developed by W. Metzenthen. It covers only its use
  14  * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
  15  * use is not permitted under this copyright.
  16  *
  17  * Redistribution and use in source and binary forms, with or without
  18  * modification, are permitted provided that the following conditions
  19  * are met:
  20  * 1. Redistributions of source code must retain the above copyright
  21  *    notice, this list of conditions and the following disclaimer.
  22  * 2. Redistributions in binary form must include information specifying
  23  *    that source code for the emulator is freely available and include
  24  *    either:
  25  *      a) an offer to provide the source code for a nominal distribution
  26  *         fee, or
  27  *      b) list at least two alternative methods whereby the source
  28  *         can be obtained, e.g. a publically accessible bulletin board
  29  *         and an anonymous ftp site from which the software can be
  30  *         downloaded.
  31  * 3. All advertising materials specifically mentioning features or use of
  32  *    this emulator must acknowledge that it was developed by W. Metzenthen.
  33  * 4. The name of W. Metzenthen may not be used to endorse or promote
  34  *    products derived from this software without specific prior written
  35  *    permission.
  36  *
  37  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
  38  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
  39  * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL
  40  * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  41  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  42  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  43  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  44  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  45  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  46  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  47  *
  48  *
  49  * The purpose of this copyright, based upon the Berkeley copyright, is to
  50  * ensure that the covered software remains freely available to everyone.
  51  *
  52  * The software (with necessary differences) is also available, but under
  53  * the terms of the GNU copyleft, for the Linux operating system and for
  54  * the djgpp ms-dos extender.
  55  *
  56  * W. Metzenthen   June 1994.
  57  *
  58  *
  59  * $FreeBSD: src/sys/gnu/i386/fpemul/poly_l2.c,v 1.10 1999/08/28 00:42:53 peter Exp $
  60  * $DragonFly: src/sys/platform/pc32/gnu/fpemul/Attic/poly_l2.c,v 1.3 2003/08/07 21:17:20 dillon Exp $
  61  *
  62  */
  63
  64
  65 #include "reg_constant.h"
  66 #include "control_w.h"
  67
  68
  69
  70 #define HIPOWER 9
  71 static unsigned short lterms[HIPOWER][4] =
  72 {
  73         /* Ideal computation with these coeffs gives about 64.6 bit rel
  74          * accuracy. */
  75         {0xe177, 0xb82f, 0x7652, 0x7154},
  76         {0xee0f, 0xe80f, 0x2770, 0x7b1c},
  77         {0x0fc0, 0xbe87, 0xb143, 0x49dd},
  78         {0x78b9, 0xdadd, 0xec54, 0x34c2},
  79         {0x003a, 0x5de9, 0x628b, 0x2909},
  80         {0x5588, 0xed16, 0x4abf, 0x2193},
  81         {0xb461, 0x85f7, 0x347a, 0x1c6a},
  82         {0x0975, 0x87b3, 0xd5bf, 0x1876},
  83         {0xe85c, 0xcec9, 0x84e7, 0x187d}
  84 };
  85
  86
  87
  88
  89 /*--- poly_l2() -------------------------------------------------------------+
  90  |   Base 2 logarithm by a polynomial approximation.                         |
  91  +---------------------------------------------------------------------------*/
  92 void
  93 poly_l2(FPU_REG * arg, FPU_REG * result)
  94 {
  95         short   exponent;
  96         char    zero;           /* flag for an Xx == 0 */
  97         unsigned short bits, shift;
  98         long long Xsq;
  99         FPU_REG accum, denom, num, Xx;
 100
 101
 102         exponent = arg->exp - EXP_BIAS;
 103
 104         accum.tag = TW_Valid;   /* set the tags to Valid */
 105
 106         if (arg->sigh > (unsigned) 0xb504f334) {
 107                 /* This is good enough for the computation of the polynomial
 108                  * sum, but actually results in a loss of precision for the
 109                  * computation of Xx. This will matter only if exponent
 110                  * becomes zero. */
 111                 exponent++;
 112                 accum.sign = 1; /* sign to negative */
 113                 num.exp = EXP_BIAS;     /* needed to prevent errors in div
 114                                          * routine */
 115                 reg_u_div(&CONST_1, arg, &num, FULL_PRECISION);
 116         } else {
 117                 accum.sign = 0; /* set the sign to positive */
 118                 num.sigl = arg->sigl;   /* copy the mantissa */
 119                 num.sigh = arg->sigh;
 120         }
 121
 122
 123         /* shift num left, lose the ms bit */
 124         num.sigh <<= 1;
 125         if (num.sigl & 0x80000000)
 126                 num.sigh |= 1;
 127         num.sigl <<= 1;
 128
 129         denom.sigl = num.sigl;
 130         denom.sigh = num.sigh;
 131         poly_div4((long long *) &(denom.sigl));
 132         denom.sigh += 0x80000000;       /* set the msb */
 133         Xx.exp = EXP_BIAS;      /* needed to prevent errors in div routine */
 134         reg_u_div(&num, &denom, &Xx, FULL_PRECISION);
 135
 136         zero = !(Xx.sigh | Xx.sigl);
 137
 138         mul64((long long *) &Xx.sigl, (long long *) &Xx.sigl, &Xsq);
 139         poly_div16(&Xsq);
 140
 141         accum.exp = -1;         /* exponent of accum */
 142
 143         /* Do the basic fixed point polynomial evaluation */
 144         polynomial((unsigned *) &accum.sigl, (unsigned *) &Xsq, lterms, HIPOWER - 1);
 145
 146         if (!exponent) {
 147                 /* If the exponent is zero, then we would lose precision by
 148                  * sticking to fixed point computation here */
 149                 /* We need to re-compute Xx because of loss of precision. */
 150                 FPU_REG lXx;
 151                 char    sign;
 152
 153                 sign = accum.sign;
 154                 accum.sign = 0;
 155
 156                 /* make accum compatible and normalize */
 157                 accum.exp = EXP_BIAS + accum.exp;
 158                 normalize(&accum);
 159
 160                 if (zero) {
 161                         reg_move(&CONST_Z, result);
 162                 } else {
 163                         /* we need to re-compute lXx to better accuracy */
 164                         num.tag = TW_Valid;     /* set the tags to Vaild */
 165                         num.sign = 0;   /* set the sign to positive */
 166                         num.exp = EXP_BIAS - 1;
 167                         if (sign) {
 168                                 /* The argument is of the form 1-x */
 169                                 /* Use  1-1/(1-x) = x/(1-x) */
 170                                 *((long long *) &num.sigl) = -*((long long *) &(arg->sigl));
 171                                 normalize(&num);
 172                                 reg_div(&num, arg, &num, FULL_PRECISION);
 173                         } else {
 174                                 normalize(&num);
 175                         }
 176
 177                         denom.tag = TW_Valid;   /* set the tags to Valid */
 178                         denom.sign = SIGN_POS;  /* set the sign to positive */
 179                         denom.exp = EXP_BIAS;
 180
 181                         reg_div(&num, &denom, &lXx, FULL_PRECISION);
 182
 183                         reg_u_mul(&lXx, &accum, &accum, FULL_PRECISION);
 184
 185                         reg_u_add(&lXx, &accum, result, FULL_PRECISION);
 186
 187                         normalize(result);
 188                 }
 189
 190                 result->sign = sign;
 191                 return;
 192         }
 193         mul64((long long *) &accum.sigl,
 194             (long long *) &Xx.sigl, (long long *) &accum.sigl);
 195
 196         *((long long *) (&accum.sigl)) += *((long long *) (&Xx.sigl));
 197
 198         if (Xx.sigh > accum.sigh) {
 199                 /* There was an overflow */
 200
 201                 poly_div2((long long *) &accum.sigl);
 202                 accum.sigh |= 0x80000000;
 203                 accum.exp++;
 204         }
 205         /* When we add the exponent to the accum result later, we will require
 206          * that their signs are the same. Here we ensure that this is so. */
 207         if (exponent && ((exponent < 0) ^ (accum.sign))) {
 208                 /* signs are different */
 209
 210                 accum.sign = !accum.sign;
 211
 212                 /* An exceptional case is when accum is zero */
 213                 if (accum.sigl | accum.sigh) {
 214                         /* find 1-accum */
 215                         /* Shift to get exponent == 0 */
 216                         if (accum.exp < 0) {
 217                                 poly_div2((long long *) &accum.sigl);
 218                                 accum.exp++;
 219                         }
 220                         /* Just negate, but throw away the sign */
 221                         *((long long *) &(accum.sigl)) = -*((long long *) &(accum.sigl));
 222                         if (exponent < 0)
 223                                 exponent++;
 224                         else
 225                                 exponent--;
 226                 }
 227         }
 228         shift = exponent >= 0 ? exponent : -exponent;
 229         bits = 0;
 230         if (shift) {
 231                 if (accum.exp) {
 232                         accum.exp++;
 233                         poly_div2((long long *) &accum.sigl);
 234                 }
 235                 while (shift) {
 236                         poly_div2((long long *) &accum.sigl);
 237                         if (shift & 1)
 238                                 accum.sigh |= 0x80000000;
 239                         shift >>= 1;
 240                         bits++;
 241                 }
 242         }
 243         /* Convert to 64 bit signed-compatible */
 244         accum.exp += bits + EXP_BIAS - 1;
 245
 246         reg_move(&accum, result);
 247         normalize(result);
 248
 249         return;
 250 }
 251
 252
 253 /*--- poly_l2p1() -----------------------------------------------------------+
 254  |   Base 2 logarithm by a polynomial approximation.                         |
 255  |   log2(x+1)                                                               |
 256  +---------------------------------------------------------------------------*/
 257 int
 258 poly_l2p1(FPU_REG * arg, FPU_REG * result)
 259 {
 260         char    sign = 0;
 261         long long Xsq;
 262         FPU_REG arg_pl1, denom, accum, local_arg, poly_arg;
 263
 264
 265         sign = arg->sign;
 266
 267         reg_add(arg, &CONST_1, &arg_pl1, FULL_PRECISION);
 268
 269         if ((arg_pl1.sign) | (arg_pl1.tag)) {   /* We need a valid positive
 270                                                  * number! */
 271                 return 1;
 272         }
 273         reg_add(&CONST_1, &arg_pl1, &denom, FULL_PRECISION);
 274         reg_div(arg, &denom, &local_arg, FULL_PRECISION);
 275         local_arg.sign = 0;     /* Make the sign positive */
 276
 277         /* Now we need to check that  |local_arg| is less than 3-2*sqrt(2) =
 278          * 0.17157.. = .0xafb0ccc0 * 2^-2 */
 279
 280         if (local_arg.exp >= EXP_BIAS - 3) {
 281                 if ((local_arg.exp > EXP_BIAS - 3) ||
 282                     (local_arg.sigh > (unsigned) 0xafb0ccc0)) {
 283                         /* The argument is large */
 284                         poly_l2(&arg_pl1, result);
 285                         return 0;
 286                 }
 287         }
 288         /* Make a copy of local_arg */
 289         reg_move(&local_arg, &poly_arg);
 290
 291         /* Get poly_arg bits aligned as required */
 292         shrx((unsigned *) &(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));
 293
 294         mul64((long long *) &(poly_arg.sigl), (long long *) &(poly_arg.sigl), &Xsq);
 295         poly_div16(&Xsq);
 296
 297         /* Do the basic fixed point polynomial evaluation */
 298         polynomial((u_int *) &accum.sigl, (unsigned *) &Xsq, lterms, HIPOWER - 1);
 299
 300         accum.tag = TW_Valid;   /* set the tags to Valid */
 301         accum.sign = SIGN_POS;  /* and make accum positive */
 302
 303         /* make accum compatible and normalize */
 304         accum.exp = EXP_BIAS - 1;
 305         normalize(&accum);
 306
 307         reg_u_mul(&local_arg, &accum, &accum, FULL_PRECISION);
 308
 309         reg_u_add(&local_arg, &accum, result, FULL_PRECISION);
 310
 311         /* Multiply the result by 2 */
 312         result->exp++;
 313
 314         result->sign = sign;
 315
 316         return 0;
 317 }