lib/msun/src/s_nexttowardf.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 #include <sys/cdefs.h>
  13 #include <float.h>
  14
  15 #include "fpmath.h"
  16 #include "math.h"
  17 #include "math_private.h"
  18
  19 #define LDBL_INFNAN_EXP (LDBL_MAX_EXP * 2 - 1)
  20
  21 float
  22 nexttowardf(float x, long double y)
  23 {
  24         union IEEEl2bits uy;
  25         volatile float t;
  26         int32_t hx,ix;
  27
  28         GET_FLOAT_WORD(hx,x);
  29         ix = hx&0x7fffffff;             /* |x| */
  30         uy.e = y;
  31
  32         if((ix>0x7f800000) ||
  33            (uy.bits.exp == LDBL_INFNAN_EXP &&
  34             ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
  35            return x+y;  /* x or y is nan */
  36         if(x==y) return (float)y;               /* x=y, return y */
  37         if(ix==0) {                             /* x == 0 */
  38             SET_FLOAT_WORD(x,(uy.bits.sign<<31)|1);/* return +-minsubnormal */
  39             t = x*x;
  40             if(t==x) return t; else return x;   /* raise underflow flag */
  41         }
  42         if(hx>=0 ^ x < y)                       /* x -= ulp */
  43             hx -= 1;
  44         else                                    /* x += ulp */
  45             hx += 1;
  46         ix = hx&0x7f800000;
  47         if(ix>=0x7f800000) return x+x;  /* overflow  */
  48         if(ix<0x00800000) {             /* underflow */
  49             t = x*x;
  50             if(t!=x) {          /* raise underflow flag */
  51                 SET_FLOAT_WORD(x,hx);
  52                 return x;
  53             }
  54         }
  55         SET_FLOAT_WORD(x,hx);
  56         return x;
  57 }