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