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 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak nearbyintf = __nearbyintf
  32 
  33 #include "libm.h"
  34 #include <fenv.h>
  35 
  36 float
  37 __nearbyintf(float x)
  38 {
  39         union {
  40                 unsigned i;
  41                 float f;
  42         } xx;
  43 
  44         unsigned hx, sx, i, frac;
  45         int rm;
  46 
  47         xx.f = x;
  48         sx = xx.i & 0x80000000;
  49         hx = xx.i & ~0x80000000;
  50 
  51         /* handle trivial cases */
  52         if (hx >= 0x4b000000) {              /* x is nan, inf, or already integral */
  53                 if (hx > 0x7f800000) /* x is nan */
  54                         return (x * x);         /* + -> * for Cheetah */
  55 
  56                 return (x);
  57         } else if (hx == 0) {                   /* x is zero */
  58                 return (x);
  59         }
  60 
  61         /* get the rounding mode */
  62         rm = fegetround();
  63 
  64         /* flip the sense of directed roundings if x is negative */
  65         if (sx && (rm == FE_UPWARD || rm == FE_DOWNWARD))
  66                 rm = (FE_UPWARD + FE_DOWNWARD) - rm;
  67 
  68         /* handle |x| < 1 */
  69         if (hx < 0x3f800000) {
  70                 if (rm == FE_UPWARD || (rm == FE_TONEAREST && hx > 0x3f000000))
  71                         xx.i = sx | 0x3f800000;
  72                 else
  73                         xx.i = sx;
  74 
  75                 return (xx.f);
  76         }
  77 
  78         /* round x at the integer bit */
  79         i = 1 << (0x96 - (hx >> 23));
  80         frac = hx & (i - 1);
  81 
  82         if (!frac)
  83                 return (x);
  84 
  85         hx &= ~(i - 1);
  86 
  87         if (rm == FE_UPWARD || (rm == FE_TONEAREST && (frac > (i >> 1) ||
  88             ((frac == (i >> 1)) && (hx & i)))))
  89                 xx.i = sx | (hx + i);
  90         else
  91                 xx.i = sx | hx;
  92 
  93         return (xx.f);
  94 }
  95 
  96 #if 0
  97 /*
  98  * Alternate implementations for SPARC, x86, using fp ops.  These may
  99  * be faster depending on how expensive saving and restoring the fp
 100  * modes and status flags is.
 101  */
 102 
 103 #include "libm.h"
 104 #include "fma.h"
 105 
 106 #if defined(__sparc)
 107 float
 108 __nearbyintf(float x)
 109 {
 110         union {
 111                 unsigned i;
 112                 float f;
 113         } xx, yy;
 114 
 115         float z;
 116         unsigned hx, sx, fsr, oldfsr;
 117         int rm;
 118 
 119         xx.f = x;
 120         sx = xx.i & 0x80000000;
 121         hx = xx.i & ~0x80000000;
 122 
 123         /* handle trivial cases */
 124         if (hx >= 0x4b000000)                /* x is nan, inf, or already integral */
 125                 return (x + 0.0f);
 126         else if (hx == 0)               /* x is zero */
 127                 return (x);
 128 
 129         /* save the fsr */
 130         __fenv_getfsr(&oldfsr);
 131 
 132         /* handle |x| < 1 */
 133         if (hx < 0x3f800000) {
 134                 /* flip the sense of directed roundings if x is negative */
 135                 rm = oldfsr >> 30;
 136 
 137                 if (sx)
 138                         rm ^= rm >> 1;
 139 
 140                 if (rm == FSR_RP || (rm == FSR_RN && hx > 0x3f000000))
 141                         xx.i = sx | 0x3f800000;
 142                 else
 143                         xx.i = sx;
 144 
 145                 return (xx.f);
 146         }
 147 
 148         /* clear the inexact trap */
 149         fsr = oldfsr & ~FSR_NXM;
 150         __fenv_setfsr(&fsr);
 151 
 152         /* round x at the integer bit */
 153         yy.i = sx | 0x4b000000;
 154         z = (x + yy.f) - yy.f;
 155 
 156         /* restore the old fsr */
 157         __fenv_setfsr(&oldfsr);
 158 
 159         return (z);
 160 }
 161 #elif defined(__x86)
 162 /* inline template */
 163 extern long double frndint(long double);
 164 
 165 float
 166 __nearbyintf(float x)
 167 {
 168         long double z;
 169         unsigned oldcwsw, cwsw;
 170 
 171         /* save the control and status words, mask the inexact exception */
 172         __fenv_getcwsw(&oldcwsw);
 173         cwsw = oldcwsw | 0x00200000;
 174         __fenv_setcwsw(&cwsw);
 175 
 176         z = frndint((long double)x);
 177 
 178         /*
 179          * restore the control and status words, preserving all but the
 180          * inexact flag
 181          */
 182         __fenv_getcwsw(&cwsw);
 183         oldcwsw |= (cwsw & 0x1f);
 184         __fenv_setcwsw(&oldcwsw);
 185 
 186         /* note: the value of z is representable in single precision */
 187         return (z);
 188 }
 189 #else
 190 #error Unknown architecture
 191 #endif
 192 #endif