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