lib/libc/sparc64/fpu/fpu_implode.c

   1 /*
   2  * Copyright (c) 1992, 1993
   3  *      The Regents of the University of California.  All rights reserved.
   4  *
   5  * This software was developed by the Computer Systems Engineering group
   6  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
   7  * contributed to Berkeley.
   8  *
   9  * All advertising materials mentioning features or use of this software
  10  * must display the following acknowledgement:
  11  *      This product includes software developed by the University of
  12  *      California, Lawrence Berkeley Laboratory.
  13  *
  14  * Redistribution and use in source and binary forms, with or without
  15  * modification, are permitted provided that the following conditions
  16  * are met:
  17  * 1. Redistributions of source code must retain the above copyright
  18  *    notice, this list of conditions and the following disclaimer.
  19  * 2. Redistributions in binary form must reproduce the above copyright
  20  *    notice, this list of conditions and the following disclaimer in the
  21  *    documentation and/or other materials provided with the distribution.
  22  * 4. Neither the name of the University nor the names of its contributors
  23  *    may be used to endorse or promote products derived from this software
  24  *    without specific prior written permission.
  25  *
  26  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  27  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  28  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  29  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  30  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  31  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  32  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  33  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  34  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  35  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  36  * SUCH DAMAGE.
  37  *
  38  *      @(#)fpu_implode.c       8.1 (Berkeley) 6/11/93
  39  *      $NetBSD: fpu_implode.c,v 1.8 2001/08/26 05:44:46 eeh Exp $
  40  */
  41
  42 #include <sys/cdefs.h>
  43 __FBSDID("$FreeBSD$");
  44
  45 /*
  46  * FPU subroutines: `implode' internal format numbers into the machine's
  47  * `packed binary' format.
  48  */
  49
  50 #include <sys/param.h>
  51
  52 #include <machine/frame.h>
  53 #include <machine/fp.h>
  54 #include <machine/fsr.h>
  55 #include <machine/ieee.h>
  56 #include <machine/instr.h>
  57
  58 #include "fpu_arith.h"
  59 #include "fpu_emu.h"
  60 #include "fpu_extern.h"
  61 #include "__sparc_utrap_private.h"
  62
  63 static int fpround(struct fpemu *, struct fpn *);
  64 static int toinf(struct fpemu *, int);
  65
  66 /*
  67  * Round a number (algorithm from Motorola MC68882 manual, modified for
  68  * our internal format).  Set inexact exception if rounding is required.
  69  * Return true iff we rounded up.
  70  *
  71  * After rounding, we discard the guard and round bits by shifting right
  72  * 2 bits (a la fpu_shr(), but we do not bother with fp->fp_sticky).
  73  * This saves effort later.
  74  *
  75  * Note that we may leave the value 2.0 in fp->fp_mant; it is the caller's
  76  * responsibility to fix this if necessary.
  77  */
  78 static int
  79 fpround(struct fpemu *fe, struct fpn *fp)
  80 {
  81         u_int m0, m1, m2, m3;
  82         int gr, s;
  83
  84         m0 = fp->fp_mant[0];
  85         m1 = fp->fp_mant[1];
  86         m2 = fp->fp_mant[2];
  87         m3 = fp->fp_mant[3];
  88         gr = m3 & 3;
  89         s = fp->fp_sticky;
  90
  91         /* mant >>= FP_NG */
  92         m3 = (m3 >> FP_NG) | (m2 << (32 - FP_NG));
  93         m2 = (m2 >> FP_NG) | (m1 << (32 - FP_NG));
  94         m1 = (m1 >> FP_NG) | (m0 << (32 - FP_NG));
  95         m0 >>= FP_NG;
  96
  97         if ((gr | s) == 0)      /* result is exact: no rounding needed */
  98                 goto rounddown;
  99
 100         fe->fe_cx |= FSR_NX;    /* inexact */
 101
 102         /* Go to rounddown to round down; break to round up. */
 103         switch (FSR_GET_RD(fe->fe_fsr)) {
 104         case FSR_RD_N:
 105         default:
 106                 /*
 107                  * Round only if guard is set (gr & 2).  If guard is set,
 108                  * but round & sticky both clear, then we want to round
 109                  * but have a tie, so round to even, i.e., add 1 iff odd.
 110                  */
 111                 if ((gr & 2) == 0)
 112                         goto rounddown;
 113                 if ((gr & 1) || fp->fp_sticky || (m3 & 1))
 114                         break;
 115                 goto rounddown;
 116
 117         case FSR_RD_Z:
 118                 /* Round towards zero, i.e., down. */
 119                 goto rounddown;
 120
 121         case FSR_RD_NINF:
 122                 /* Round towards -Inf: up if negative, down if positive. */
 123                 if (fp->fp_sign)
 124                         break;
 125                 goto rounddown;
 126
 127         case FSR_RD_PINF:
 128                 /* Round towards +Inf: up if positive, down otherwise. */
 129                 if (!fp->fp_sign)
 130                         break;
 131                 goto rounddown;
 132         }
 133
 134         /* Bump low bit of mantissa, with carry. */
 135         FPU_ADDS(m3, m3, 1);
 136         FPU_ADDCS(m2, m2, 0);
 137         FPU_ADDCS(m1, m1, 0);
 138         FPU_ADDC(m0, m0, 0);
 139         fp->fp_mant[0] = m0;
 140         fp->fp_mant[1] = m1;
 141         fp->fp_mant[2] = m2;
 142         fp->fp_mant[3] = m3;
 143         return (1);
 144
 145 rounddown:
 146         fp->fp_mant[0] = m0;
 147         fp->fp_mant[1] = m1;
 148         fp->fp_mant[2] = m2;
 149         fp->fp_mant[3] = m3;
 150         return (0);
 151 }
 152
 153 /*
 154  * For overflow: return true if overflow is to go to +/-Inf, according
 155  * to the sign of the overflowing result.  If false, overflow is to go
 156  * to the largest magnitude value instead.
 157  */
 158 static int
 159 toinf(struct fpemu *fe, int sign)
 160 {
 161         int inf;
 162
 163         /* look at rounding direction */
 164         switch (FSR_GET_RD(fe->fe_fsr)) {
 165         default:
 166         case FSR_RD_N:          /* the nearest value is always Inf */
 167                 inf = 1;
 168                 break;
 169
 170         case FSR_RD_Z:          /* toward 0 => never towards Inf */
 171                 inf = 0;
 172                 break;
 173
 174         case FSR_RD_PINF:       /* toward +Inf iff positive */
 175                 inf = sign == 0;
 176                 break;
 177
 178         case FSR_RD_NINF:       /* toward -Inf iff negative */
 179                 inf = sign;
 180                 break;
 181         }
 182         return (inf);
 183 }
 184
 185 /*
 186  * fpn -> int (int value returned as return value).
 187  *
 188  * N.B.: this conversion always rounds towards zero (this is a peculiarity
 189  * of the SPARC instruction set).
 190  */
 191 u_int
 192 __fpu_ftoi(fe, fp)
 193         struct fpemu *fe;
 194         struct fpn *fp;
 195 {
 196         u_int i;
 197         int sign, exp;
 198
 199         sign = fp->fp_sign;
 200         switch (fp->fp_class) {
 201
 202         case FPC_ZERO:
 203                 return (0);
 204
 205         case FPC_NUM:
 206                 /*
 207                  * If exp >= 2^32, overflow.  Otherwise shift value right
 208                  * into last mantissa word (this will not exceed 0xffffffff),
 209                  * shifting any guard and round bits out into the sticky
 210                  * bit.  Then ``round'' towards zero, i.e., just set an
 211                  * inexact exception if sticky is set (see round()).
 212                  * If the result is > 0x80000000, or is positive and equals
 213                  * 0x80000000, overflow; otherwise the last fraction word
 214                  * is the result.
 215                  */
 216                 if ((exp = fp->fp_exp) >= 32)
 217                         break;
 218                 /* NB: the following includes exp < 0 cases */
 219                 if (__fpu_shr(fp, FP_NMANT - 1 - exp) != 0)
 220                         fe->fe_cx |= FSR_NX;
 221                 i = fp->fp_mant[3];
 222                 if (i >= ((u_int)0x80000000 + sign))
 223                         break;
 224                 return (sign ? -i : i);
 225
 226         default:                /* Inf, qNaN, sNaN */
 227                 break;
 228         }
 229         /* overflow: replace any inexact exception with invalid */
 230         fe->fe_cx = (fe->fe_cx & ~FSR_NX) | FSR_NV;
 231         return (0x7fffffff + sign);
 232 }
 233
 234 /*
 235  * fpn -> extended int (high bits of int value returned as return value).
 236  *
 237  * N.B.: this conversion always rounds towards zero (this is a peculiarity
 238  * of the SPARC instruction set).
 239  */
 240 u_int
 241 __fpu_ftox(fe, fp, res)
 242         struct fpemu *fe;
 243         struct fpn *fp;
 244         u_int *res;
 245 {
 246         u_int64_t i;
 247         int sign, exp;
 248
 249         sign = fp->fp_sign;
 250         switch (fp->fp_class) {
 251
 252         case FPC_ZERO:
 253                 res[1] = 0;
 254                 return (0);
 255
 256         case FPC_NUM:
 257                 /*
 258                  * If exp >= 2^64, overflow.  Otherwise shift value
 259                  * right into last mantissa word (this will not exceed
 260                  * 0xffffffffffffffff), shifting any guard and round
 261                  * bits out into the sticky bit.  Then ``round'' towards
 262                  * zero, i.e., just set an inexact exception if sticky
 263                  * is set (see round()).
 264                  * If the result is > 0x8000000000000000, or is positive
 265                  * and equals 0x8000000000000000, overflow; otherwise
 266                  * the last fraction word is the result.
 267                  */
 268                 if ((exp = fp->fp_exp) >= 64)
 269                         break;
 270                 /* NB: the following includes exp < 0 cases */
 271                 if (__fpu_shr(fp, FP_NMANT - 1 - exp) != 0)
 272                         fe->fe_cx |= FSR_NX;
 273                 i = ((u_int64_t)fp->fp_mant[2]<<32)|fp->fp_mant[3];
 274                 if (i >= ((u_int64_t)0x8000000000000000LL + sign))
 275                         break;
 276                 if (sign)
 277                         i = -i;
 278                 res[1] = (int)i;
 279                 return (i >> 32);
 280
 281         default:                /* Inf, qNaN, sNaN */
 282                 break;
 283         }
 284         /* overflow: replace any inexact exception with invalid */
 285         fe->fe_cx = (fe->fe_cx & ~FSR_NX) | FSR_NV;
 286         return (0x7fffffffffffffffLL + sign);
 287 }
 288
 289 /*
 290  * fpn -> single (32 bit single returned as return value).
 291  * We assume <= 29 bits in a single-precision fraction (1.f part).
 292  */
 293 u_int
 294 __fpu_ftos(fe, fp)
 295         struct fpemu *fe;
 296         struct fpn *fp;
 297 {
 298         u_int sign = fp->fp_sign << 31;
 299         int exp;
 300
 301 #define SNG_EXP(e)      ((e) << SNG_FRACBITS)   /* makes e an exponent */
 302 #define SNG_MASK        (SNG_EXP(1) - 1)        /* mask for fraction */
 303
 304         /* Take care of non-numbers first. */
 305         if (ISNAN(fp)) {
 306                 /*
 307                  * Preserve upper bits of NaN, per SPARC V8 appendix N.
 308                  * Note that fp->fp_mant[0] has the quiet bit set,
 309                  * even if it is classified as a signalling NaN.
 310                  */
 311                 (void) __fpu_shr(fp, FP_NMANT - 1 - SNG_FRACBITS);
 312                 exp = SNG_EXP_INFNAN;
 313                 goto done;
 314         }
 315         if (ISINF(fp))
 316                 return (sign | SNG_EXP(SNG_EXP_INFNAN));
 317         if (ISZERO(fp))
 318                 return (sign);
 319
 320         /*
 321          * Normals (including subnormals).  Drop all the fraction bits
 322          * (including the explicit ``implied'' 1 bit) down into the
 323          * single-precision range.  If the number is subnormal, move
 324          * the ``implied'' 1 into the explicit range as well, and shift
 325          * right to introduce leading zeroes.  Rounding then acts
 326          * differently for normals and subnormals: the largest subnormal
 327          * may round to the smallest normal (1.0 x 2^minexp), or may
 328          * remain subnormal.  In the latter case, signal an underflow
 329          * if the result was inexact or if underflow traps are enabled.
 330          *
 331          * Rounding a normal, on the other hand, always produces another
 332          * normal (although either way the result might be too big for
 333          * single precision, and cause an overflow).  If rounding a
 334          * normal produces 2.0 in the fraction, we need not adjust that
 335          * fraction at all, since both 1.0 and 2.0 are zero under the
 336          * fraction mask.
 337          *
 338          * Note that the guard and round bits vanish from the number after
 339          * rounding.
 340          */
 341         if ((exp = fp->fp_exp + SNG_EXP_BIAS) <= 0) {   /* subnormal */
 342                 /* -NG for g,r; -SNG_FRACBITS-exp for fraction */
 343                 (void) __fpu_shr(fp, FP_NMANT - FP_NG - SNG_FRACBITS - exp);
 344                 if (fpround(fe, fp) && fp->fp_mant[3] == SNG_EXP(1))
 345                         return (sign | SNG_EXP(1) | 0);
 346                 if ((fe->fe_cx & FSR_NX) ||
 347                     (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT)))
 348                         fe->fe_cx |= FSR_UF;
 349                 return (sign | SNG_EXP(0) | fp->fp_mant[3]);
 350         }
 351         /* -FP_NG for g,r; -1 for implied 1; -SNG_FRACBITS for fraction */
 352         (void) __fpu_shr(fp, FP_NMANT - FP_NG - 1 - SNG_FRACBITS);
 353 #ifdef DIAGNOSTIC
 354         if ((fp->fp_mant[3] & SNG_EXP(1 << FP_NG)) == 0)
 355                 __utrap_panic("fpu_ftos");
 356 #endif
 357         if (fpround(fe, fp) && fp->fp_mant[3] == SNG_EXP(2))
 358                 exp++;
 359         if (exp >= SNG_EXP_INFNAN) {
 360                 /* overflow to inf or to max single */
 361                 fe->fe_cx |= FSR_OF | FSR_NX;
 362                 if (toinf(fe, sign))
 363                         return (sign | SNG_EXP(SNG_EXP_INFNAN));
 364                 return (sign | SNG_EXP(SNG_EXP_INFNAN - 1) | SNG_MASK);
 365         }
 366 done:
 367         /* phew, made it */
 368         return (sign | SNG_EXP(exp) | (fp->fp_mant[3] & SNG_MASK));
 369 }
 370
 371 /*
 372  * fpn -> double (32 bit high-order result returned; 32-bit low order result
 373  * left in res[1]).  Assumes <= 61 bits in double precision fraction.
 374  *
 375  * This code mimics fpu_ftos; see it for comments.
 376  */
 377 u_int
 378 __fpu_ftod(fe, fp, res)
 379         struct fpemu *fe;
 380         struct fpn *fp;
 381         u_int *res;
 382 {
 383         u_int sign = fp->fp_sign << 31;
 384         int exp;
 385
 386 #define DBL_EXP(e)      ((e) << (DBL_FRACBITS & 31))
 387 #define DBL_MASK        (DBL_EXP(1) - 1)
 388
 389         if (ISNAN(fp)) {
 390                 (void) __fpu_shr(fp, FP_NMANT - 1 - DBL_FRACBITS);
 391                 exp = DBL_EXP_INFNAN;
 392                 goto done;
 393         }
 394         if (ISINF(fp)) {
 395                 sign |= DBL_EXP(DBL_EXP_INFNAN);
 396                 goto zero;
 397         }
 398         if (ISZERO(fp)) {
 399 zero:           res[1] = 0;
 400                 return (sign);
 401         }
 402
 403         if ((exp = fp->fp_exp + DBL_EXP_BIAS) <= 0) {
 404                 (void) __fpu_shr(fp, FP_NMANT - FP_NG - DBL_FRACBITS - exp);
 405                 if (fpround(fe, fp) && fp->fp_mant[2] == DBL_EXP(1)) {
 406                         res[1] = 0;
 407                         return (sign | DBL_EXP(1) | 0);
 408                 }
 409                 if ((fe->fe_cx & FSR_NX) ||
 410                     (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT)))
 411                         fe->fe_cx |= FSR_UF;
 412                 exp = 0;
 413                 goto done;
 414         }
 415         (void) __fpu_shr(fp, FP_NMANT - FP_NG - 1 - DBL_FRACBITS);
 416         if (fpround(fe, fp) && fp->fp_mant[2] == DBL_EXP(2))
 417                 exp++;
 418         if (exp >= DBL_EXP_INFNAN) {
 419                 fe->fe_cx |= FSR_OF | FSR_NX;
 420                 if (toinf(fe, sign)) {
 421                         res[1] = 0;
 422                         return (sign | DBL_EXP(DBL_EXP_INFNAN) | 0);
 423                 }
 424                 res[1] = ~0;
 425                 return (sign | DBL_EXP(DBL_EXP_INFNAN) | DBL_MASK);
 426         }
 427 done:
 428         res[1] = fp->fp_mant[3];
 429         return (sign | DBL_EXP(exp) | (fp->fp_mant[2] & DBL_MASK));
 430 }
 431
 432 /*
 433  * fpn -> extended (32 bit high-order result returned; low-order fraction
 434  * words left in res[1]..res[3]).  Like ftod, which is like ftos ... but
 435  * our internal format *is* extended precision, plus 2 bits for guard/round,
 436  * so we can avoid a small bit of work.
 437  */
 438 u_int
 439 __fpu_ftoq(fe, fp, res)
 440         struct fpemu *fe;
 441         struct fpn *fp;
 442         u_int *res;
 443 {
 444         u_int sign = fp->fp_sign << 31;
 445         int exp;
 446
 447 #define EXT_EXP(e)      ((e) << (EXT_FRACBITS & 31))
 448 #define EXT_MASK        (EXT_EXP(1) - 1)
 449
 450         if (ISNAN(fp)) {
 451                 (void) __fpu_shr(fp, 2);        /* since we are not rounding */
 452                 exp = EXT_EXP_INFNAN;
 453                 goto done;
 454         }
 455         if (ISINF(fp)) {
 456                 sign |= EXT_EXP(EXT_EXP_INFNAN);
 457                 goto zero;
 458         }
 459         if (ISZERO(fp)) {
 460 zero:           res[1] = res[2] = res[3] = 0;
 461                 return (sign);
 462         }
 463
 464         if ((exp = fp->fp_exp + EXT_EXP_BIAS) <= 0) {
 465                 (void) __fpu_shr(fp, FP_NMANT - FP_NG - EXT_FRACBITS - exp);
 466                 if (fpround(fe, fp) && fp->fp_mant[0] == EXT_EXP(1)) {
 467                         res[1] = res[2] = res[3] = 0;
 468                         return (sign | EXT_EXP(1) | 0);
 469                 }
 470                 if ((fe->fe_cx & FSR_NX) ||
 471                     (fe->fe_fsr & (FSR_UF << FSR_TEM_SHIFT)))
 472                         fe->fe_cx |= FSR_UF;
 473                 exp = 0;
 474                 goto done;
 475         }
 476         /* Since internal == extended, no need to shift here. */
 477         if (fpround(fe, fp) && fp->fp_mant[0] == EXT_EXP(2))
 478                 exp++;
 479         if (exp >= EXT_EXP_INFNAN) {
 480                 fe->fe_cx |= FSR_OF | FSR_NX;
 481                 if (toinf(fe, sign)) {
 482                         res[1] = res[2] = res[3] = 0;
 483                         return (sign | EXT_EXP(EXT_EXP_INFNAN) | 0);
 484                 }
 485                 res[1] = res[2] = res[3] = ~0;
 486                 return (sign | EXT_EXP(EXT_EXP_INFNAN) | EXT_MASK);
 487         }
 488 done:
 489         res[1] = fp->fp_mant[1];
 490         res[2] = fp->fp_mant[2];
 491         res[3] = fp->fp_mant[3];
 492         return (sign | EXT_EXP(exp) | (fp->fp_mant[0] & EXT_MASK));
 493 }
 494
 495 /*
 496  * Implode an fpn, writing the result into the given space.
 497  */
 498 void
 499 __fpu_implode(fe, fp, type, space)
 500         struct fpemu *fe;
 501         struct fpn *fp;
 502         int type;
 503         u_int *space;
 504 {
 505
 506         switch (type) {
 507
 508         case FTYPE_LNG:
 509                 space[0] = __fpu_ftox(fe, fp, space);
 510                 break;
 511
 512         case FTYPE_INT:
 513                 space[0] = __fpu_ftoi(fe, fp);
 514                 break;
 515
 516         case FTYPE_SNG:
 517                 space[0] = __fpu_ftos(fe, fp);
 518                 break;
 519
 520         case FTYPE_DBL:
 521                 space[0] = __fpu_ftod(fe, fp, space);
 522                 break;
 523
 524         case FTYPE_EXT:
 525                 /* funky rounding precision options ?? */
 526                 space[0] = __fpu_ftoq(fe, fp, space);
 527                 break;
 528
 529         default:
 530                 __utrap_panic("fpu_implode");
 531         }
 532         DPRINTF(FPE_REG, ("fpu_implode: %x %x %x %x\n",
 533                 space[0], space[1], space[2], space[3]));
 534 }