lib/libm/src/s_nexttowardf.c

   1 /*      $NetBSD: s_nexttowardf.c,v 1.3 2013/02/09 23:14:44 christos Exp $       */
   2
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14 #include <sys/cdefs.h>
  15 #if 0
  16 __FBSDID("$FreeBSD: src/lib/msun/src/s_nexttowardf.c,v 1.3 2011/02/10 07:38:38 das Exp $");
  17 #else
  18 __RCSID("$NetBSD: s_nexttowardf.c,v 1.3 2013/02/09 23:14:44 christos Exp $");
  19 #endif
  20
  21 #include <string.h>
  22 #include <float.h>
  23 #include <machine/ieee.h>
  24
  25 #include "math.h"
  26 #include "math_private.h"
  27
  28 #ifdef EXT_EXP_INFNAN
  29 float
  30 nexttowardf(float x, long double y)
  31 {
  32         volatile float t;
  33         int32_t hx,ix;
  34         union ieee_ext_u uy;
  35
  36         GET_FLOAT_WORD(hx,x);
  37         ix = hx&0x7fffffff;             /* |x| */
  38
  39         memset(&uy, 0, sizeof(uy));
  40         uy.extu_ld = y;
  41         uy.extu_ext.ext_frach &= ~0x80000000;
  42
  43         if((ix>0x7f800000) ||
  44            (uy.extu_ext.ext_exp == EXT_EXP_INFNAN &&
  45             (uy.extu_ext.ext_frach | uy.extu_ext.ext_fracl) != 0))
  46            return x+y;  /* x or y is nan */
  47         if(x==y) return (float)y;               /* x=y, return y */
  48         if(ix==0) {                             /* x == 0 */
  49             SET_FLOAT_WORD(x,(uy.extu_ext.ext_sign<<31)|1);/* return +-minsubnormal */
  50             t = x*x;
  51             if(t==x) return t; else return x;   /* raise underflow flag */
  52         }
  53         if((hx >= 0) ^ (x < y))                 /* x -= ulp */
  54             hx -= 1;
  55         else                                    /* x += ulp */
  56             hx += 1;
  57         ix = hx&0x7f800000;
  58         if(ix>=0x7f800000) return x+x;  /* overflow  */
  59         if(ix<0x00800000) {             /* underflow */
  60             t = x*x;
  61             if(t!=x) {          /* raise underflow flag */
  62                 SET_FLOAT_WORD(x,hx);
  63                 return x;
  64             }
  65         }
  66         SET_FLOAT_WORD(x,hx);
  67         return x;
  68 }
  69 #endif