usr/src/lib/libm/common/m9x/nearbyintf.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29
  30 #pragma weak nearbyintf = __nearbyintf
  31
  32 #include "libm.h"
  33 #include <fenv.h>
  34
  35 float
  36 __nearbyintf(float x) {
  37         union {
  38                 unsigned i;
  39                 float f;
  40         } xx;
  41         unsigned hx, sx, i, frac;
  42         int rm;
  43
  44         xx.f = x;
  45         sx = xx.i & 0x80000000;
  46         hx = xx.i & ~0x80000000;
  47
  48         /* handle trivial cases */
  49         if (hx >= 0x4b000000) { /* x is nan, inf, or already integral */
  50                 if (hx > 0x7f800000)    /* x is nan */
  51                         return (x * x);         /* + -> * for Cheetah */
  52                 return (x);
  53         } else if (hx == 0)             /* x is zero */
  54                 return (x);
  55
  56         /* get the rounding mode */
  57         rm = fegetround();
  58
  59         /* flip the sense of directed roundings if x is negative */
  60         if (sx && (rm == FE_UPWARD || rm == FE_DOWNWARD))
  61                 rm = (FE_UPWARD + FE_DOWNWARD) - rm;
  62
  63         /* handle |x| < 1 */
  64         if (hx < 0x3f800000) {
  65                 if (rm == FE_UPWARD || (rm == FE_TONEAREST && hx > 0x3f000000))
  66                         xx.i = sx | 0x3f800000;
  67                 else
  68                         xx.i = sx;
  69                 return (xx.f);
  70         }
  71
  72         /* round x at the integer bit */
  73         i = 1 << (0x96 - (hx >> 23));
  74         frac = hx & (i - 1);
  75         if (!frac)
  76                 return (x);
  77
  78         hx &= ~(i - 1);
  79         if (rm == FE_UPWARD || (rm == FE_TONEAREST && (frac > (i >> 1) ||
  80                 ((frac == (i >> 1)) && (hx & i)))))
  81                 xx.i = sx | (hx + i);
  82         else
  83                 xx.i = sx | hx;
  84         return (xx.f);
  85 }
  86
  87 #if 0
  88
  89 /*
  90  * Alternate implementations for SPARC, x86, using fp ops.  These may
  91  * be faster depending on how expensive saving and restoring the fp
  92  * modes and status flags is.
  93  */
  94
  95 #include "libm.h"
  96 #include "fma.h"
  97
  98 #if defined(__sparc)
  99
 100 float
 101 __nearbyintf(float x) {
 102         union {
 103                 unsigned i;
 104                 float f;
 105         } xx, yy;
 106         float z;
 107         unsigned hx, sx, fsr, oldfsr;
 108         int rm;
 109
 110         xx.f = x;
 111         sx = xx.i & 0x80000000;
 112         hx = xx.i & ~0x80000000;
 113
 114         /* handle trivial cases */
 115         if (hx >= 0x4b000000)   /* x is nan, inf, or already integral */
 116                 return (x + 0.0f);
 117         else if (hx == 0)       /* x is zero */
 118                 return (x);
 119
 120         /* save the fsr */
 121         __fenv_getfsr(&oldfsr);
 122
 123         /* handle |x| < 1 */
 124         if (hx < 0x3f800000) {
 125                 /* flip the sense of directed roundings if x is negative */
 126                 rm = oldfsr >> 30;
 127                 if (sx)
 128                         rm ^= rm >> 1;
 129                 if (rm == FSR_RP || (rm == FSR_RN && hx > 0x3f000000))
 130                         xx.i = sx | 0x3f800000;
 131                 else
 132                         xx.i = sx;
 133                 return (xx.f);
 134         }
 135
 136         /* clear the inexact trap */
 137         fsr = oldfsr & ~FSR_NXM;
 138         __fenv_setfsr(&fsr);
 139
 140         /* round x at the integer bit */
 141         yy.i = sx | 0x4b000000;
 142         z = (x + yy.f) - yy.f;
 143
 144         /* restore the old fsr */
 145         __fenv_setfsr(&oldfsr);
 146
 147         return (z);
 148 }
 149
 150 #elif defined(__x86)
 151
 152 /* inline template */
 153 extern long double frndint(long double);
 154
 155 float
 156 __nearbyintf(float x) {
 157         long double z;
 158         unsigned oldcwsw, cwsw;
 159
 160         /* save the control and status words, mask the inexact exception */
 161         __fenv_getcwsw(&oldcwsw);
 162         cwsw = oldcwsw | 0x00200000;
 163         __fenv_setcwsw(&cwsw);
 164
 165         z = frndint((long double) x);
 166
 167         /*
 168          * restore the control and status words, preserving all but the
 169          * inexact flag
 170          */
 171         __fenv_getcwsw(&cwsw);
 172         oldcwsw |= (cwsw & 0x1f);
 173         __fenv_setcwsw(&oldcwsw);
 174
 175         /* note: the value of z is representable in single precision */
 176         return (z);
 177 }
 178
 179 #else
 180 #error Unknown architecture
 181 #endif
 182
 183 #endif