nsz Git - libm/blob - src/math/nextafter.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_nextafter.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /* IEEE functions
  13  *      nextafter(x,y)
  14  *      return the next machine floating-point number of x in the
  15  *      direction toward y.
  16  *   Special cases:
  17  */
  18
  19 #include "libm.h"
  20
  21 double nextafter(double x, double y)
  22 {
  23         volatile double t;
  24         int32_t hx,hy,ix,iy;
  25         uint32_t lx,ly;
  26
  27         EXTRACT_WORDS(hx, lx, x);
  28         EXTRACT_WORDS(hy, ly, y);
  29         ix = hx & 0x7fffffff;  /* |x| */
  30         iy = hy & 0x7fffffff;  /* |y| */
  31
  32         if ((ix >= 0x7ff00000 && (ix-0x7ff00000)|lx) != 0 ||   /* x is nan */
  33             (iy >= 0x7ff00000 && (iy-0x7ff00000)|ly) != 0)     /* y is nan */
  34                 return x+y;
  35         if (x == y)          /* x == y */
  36                 return y;
  37         if ((ix|lx) == 0) {  /* x == 0 */
  38                 INSERT_WORDS(x, hy&0x80000000, 1);  /* return +-minsubnormal */
  39                 /* raise underflow flag */
  40                 t = x*x;
  41                 if (t == x)
  42                         return t;
  43                 return x;
  44         }
  45         if (hx >= 0) {  /* x > 0 */
  46                 if (hx > hy || (hx == hy && lx > ly)) {  /* x > y, x -= ulp */
  47                         if (lx == 0)
  48                                 hx -= 1;
  49                         lx -= 1;
  50                 } else {                                 /* x < y, x += ulp */
  51                         lx += 1;
  52                         if (lx == 0)
  53                                 hx += 1;
  54                 }
  55         } else {        /* x < 0 */
  56                 if (hy >= 0 || hx > hy || (hx == hy && lx > ly)) { /* x < y, x -= ulp */
  57                         if (lx == 0)
  58                                 hx -= 1;
  59                         lx -= 1;
  60                 } else {                                 /* x > y, x += ulp */
  61                         lx += 1;
  62                         if (lx == 0)
  63                                 hx += 1;
  64                 }
  65         }
  66         hy = hx & 0x7ff00000;
  67         if (hy >= 0x7ff00000)  /* overflow  */
  68                 return x+x;
  69         if (hy < 0x00100000) { /* underflow */
  70                 /* raise underflow flag */
  71                 t = x*x;
  72                 if (t != x) {
  73                         INSERT_WORDS(y, hx, lx);
  74                         return y;
  75                 }
  76         }
  77         INSERT_WORDS(x, hx, lx);
  78         return x;
  79 }