nsz Git - musl/blob - src/math/nexttoward.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_nexttoward.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
  13 #include "libm.h"
  14
  15 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  16 double nexttoward(double x, long double y)
  17 {
  18         return nextafter(x, y);
  19 }
  20 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
  21 double nexttoward(double x, long double y)
  22 {
  23         union IEEEl2bits uy;
  24         volatile double t;
  25         int32_t hx,ix;
  26         uint32_t lx;
  27
  28         EXTRACT_WORDS(hx, lx, x);
  29         ix = hx & 0x7fffffff;
  30         uy.e = y;
  31
  32         if ((ix >= 0x7ff00000 && ((ix-0x7ff00000)|lx) != 0) ||
  33             (uy.bits.exp == 0x7fff && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
  34                 return x + y;  /* x or y is nan */
  35         if (x == y)
  36                 return (double)y;
  37         if (x == 0.0) {
  38                 INSERT_WORDS(x, uy.bits.sign<<31, 1);  /* return +-minsubnormal */
  39                 /* raise underflow */
  40                 t = x * x;
  41                 if (t == x)
  42                         return t;
  43                 return x;
  44         }
  45         if (hx > 0.0 ^ x < y) {  /* x -= ulp */
  46                 if (lx == 0)
  47                         hx--;
  48                 lx--;
  49         } else {                 /* x += ulp */
  50                 lx++;
  51                 if (lx == 0)
  52                         hx++;
  53         }
  54         ix = hx & 0x7ff00000;
  55         if (ix >= 0x7ff00000)   /* overflow  */
  56                 return x + x;
  57         if (ix < 0x00100000) {  /* underflow */
  58                 /* raise underflow flag */
  59                 t = x * x;
  60                 if (t != x) {
  61                         INSERT_WORDS(x, hx, lx);
  62                         return x;
  63                 }
  64         }
  65         INSERT_WORDS(x, hx, lx);
  66         return x;
  67 }
  68 #endif