sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c

   1 /* s_nextafterl.c -- long double version of s_nextafter.c.
   2  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
   3  */
   4
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #if defined(LIBM_SCCS) && !defined(lint)
  17 static char rcsid[] = "$NetBSD: $";
  18 #endif
  19
  20 /* IEEE functions
  21  *      nextafterl(x,y)
  22  *      return the next machine floating-point number of x in the
  23  *      direction toward y.
  24  *   Special cases:
  25  */
  26
  27 #include <math.h>
  28 #include <math_private.h>
  29 #include <math_ldbl_opt.h>
  30
  31 long double __nextafterl(long double x, long double y)
  32 {
  33         int64_t hx, hy, ihx, ihy, lx;
  34         double xhi, xlo, yhi;
  35
  36         ldbl_unpack (x, &xhi, &xlo);
  37         EXTRACT_WORDS64 (hx, xhi);
  38         EXTRACT_WORDS64 (lx, xlo);
  39         yhi = ldbl_high (y);
  40         EXTRACT_WORDS64 (hy, yhi);
  41         ihx = hx&0x7fffffffffffffffLL;          /* |hx| */
  42         ihy = hy&0x7fffffffffffffffLL;          /* |hy| */
  43
  44         if((ihx>0x7ff0000000000000LL) ||        /* x is nan */
  45            (ihy>0x7ff0000000000000LL))          /* y is nan */
  46             return x+y; /* signal the nan */
  47         if(x==y)
  48             return y;           /* x=y, return y */
  49         if(ihx == 0) {                          /* x == 0 */
  50             long double u;                      /* return +-minsubnormal */
  51             hy = (hy & 0x8000000000000000ULL) | 1;
  52             INSERT_WORDS64 (yhi, hy);
  53             x = yhi;
  54             u = math_opt_barrier (x);
  55             u = u * u;
  56             math_force_eval (u);                /* raise underflow flag */
  57             return x;
  58         }
  59
  60         long double u;
  61         if(x > y) {     /* x > y, x -= ulp */
  62             /* This isn't the largest magnitude correctly rounded
  63                long double as you can see from the lowest mantissa
  64                bit being zero.  It is however the largest magnitude
  65                long double with a 106 bit mantissa, and nextafterl
  66                is insane with variable precision.  So to make
  67                nextafterl sane we assume 106 bit precision.  */
  68             if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL)) {
  69               u = x+x;  /* overflow, return -inf */
  70               math_force_eval (u);
  71               return y;
  72             }
  73             if (hx >= 0x7ff0000000000000LL) {
  74               u = 0x1.fffffffffffff7ffffffffffff8p+1023L;
  75               return u;
  76             }
  77             if(ihx <= 0x0360000000000000LL) {  /* x <= LDBL_MIN */
  78               u = math_opt_barrier (x);
  79               x -= __LDBL_DENORM_MIN__;
  80               if (ihx < 0x0360000000000000LL
  81                   || (hx > 0 && lx <= 0)
  82                   || (hx < 0 && lx > 1)) {
  83                 u = u * u;
  84                 math_force_eval (u);            /* raise underflow flag */
  85               }
  86               return x;
  87             }
  88             /* If the high double is an exact power of two and the low
  89                double is the opposite sign, then 1ulp is one less than
  90                what we might determine from the high double.  Similarly
  91                if X is an exact power of two, and positive, because
  92                making it a little smaller will result in the exponent
  93                decreasing by one and normalisation of the mantissa.   */
  94             if ((hx & 0x000fffffffffffffLL) == 0
  95                 && ((lx != 0 && (hx ^ lx) < 0)
  96                     || (lx == 0 && hx >= 0)))
  97               ihx -= 1LL << 52;
  98             if (ihx < (106LL << 52)) { /* ulp will denormal */
  99               INSERT_WORDS64 (yhi, ihx & (0x7ffLL<<52));
 100               u = yhi * 0x1p-105;
 101             } else {
 102               INSERT_WORDS64 (yhi, (ihx & (0x7ffLL<<52))-(105LL<<52));
 103               u = yhi;
 104             }
 105             return x - u;
 106         } else {                                /* x < y, x += ulp */
 107             if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL)) {
 108               u = x+x;  /* overflow, return +inf */
 109               math_force_eval (u);
 110               return y;
 111             }
 112             if ((uint64_t) hx >= 0xfff0000000000000ULL) {
 113               u = -0x1.fffffffffffff7ffffffffffff8p+1023L;
 114               return u;
 115             }
 116             if(ihx <= 0x0360000000000000LL) {  /* x <= LDBL_MIN */
 117               u = math_opt_barrier (x);
 118               x += __LDBL_DENORM_MIN__;
 119               if (ihx < 0x0360000000000000LL
 120                   || (hx > 0 && lx < 0 && lx != 0x8000000000000001LL)
 121                   || (hx < 0 && lx >= 0)) {
 122                 u = u * u;
 123                 math_force_eval (u);            /* raise underflow flag */
 124               }
 125               if (x == 0.0L)    /* handle negative __LDBL_DENORM_MIN__ case */
 126                 x = -0.0L;
 127               return x;
 128             }
 129             /* If the high double is an exact power of two and the low
 130                double is the opposite sign, then 1ulp is one less than
 131                what we might determine from the high double.  Similarly
 132                if X is an exact power of two, and negative, because
 133                making it a little larger will result in the exponent
 134                decreasing by one and normalisation of the mantissa.   */
 135             if ((hx & 0x000fffffffffffffLL) == 0
 136                 && ((lx != 0 && (hx ^ lx) < 0)
 137                     || (lx == 0 && hx < 0)))
 138               ihx -= 1LL << 52;
 139             if (ihx < (106LL << 52)) { /* ulp will denormal */
 140               INSERT_WORDS64 (yhi, ihx & (0x7ffLL<<52));
 141               u = yhi * 0x1p-105;
 142             } else {
 143               INSERT_WORDS64 (yhi, (ihx & (0x7ffLL<<52))-(105LL<<52));
 144               u = yhi;
 145             }
 146             return x + u;
 147         }
 148 }
 149 strong_alias (__nextafterl, __nexttowardl)
 150 long_double_symbol (libm, __nextafterl, nextafterl);
 151 long_double_symbol (libm, __nexttowardl, nexttowardl);