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

   1 #include "libm.h"
   2
   3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
   4 long double nextafterl(long double x, long double y)
   5 {
   6         return nextafter(x, y);
   7 }
   8 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
   9 #define MSB ((uint64_t)1<<63)
  10 long double nextafterl(long double x, long double y)
  11 {
  12         union ldshape ux, uy;
  13
  14         if (isnan(x) || isnan(y))
  15                 return x + y;
  16         if (x == y)
  17                 return y;
  18         ux.value = x;
  19         if (x == 0) {
  20                 uy.value = y;
  21                 ux.bits.m = 1;
  22                 ux.bits.sign = uy.bits.sign;
  23         } else if (x < y ^ ux.bits.sign) {
  24                 ux.bits.m++;
  25                 if ((ux.bits.m & ~MSB) == 0) {
  26                         ux.bits.m = MSB;
  27                         ux.bits.exp++;
  28                 }
  29         } else {
  30                 if ((ux.bits.m & ~MSB) == 0) {
  31                         ux.bits.exp--;
  32                         if (ux.bits.exp)
  33                                 ux.bits.m = 0;
  34                 }
  35                 ux.bits.m--;
  36         }
  37         /* raise overflow if ux.value is infinite and x is finite */
  38         if (ux.bits.exp == 0x7fff)
  39                 return x + x;
  40         /* raise underflow if ux.value is subnormal or zero */
  41         if (ux.bits.exp == 0) {
  42                 volatile float z;
  43                 z = x*x + ux.value*ux.value;
  44         }
  45         return ux.value;
  46 }
  47 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
  48 long double nextafterl(long double x, long double y)
  49 {
  50         union ldshape ux, uy;
  51
  52         if (isnan(x) || isnan(y))
  53                 return x + y;
  54         if (x == y)
  55                 return y;
  56         ux.value = x;
  57         if (x == 0) {
  58                 uy.value = y;
  59                 ux.bits.mlo = 1;
  60                 ux.bits.sign = uy.bits.sign;
  61         } else if (x < y ^ ux.bits.sign) {
  62                 ux.bits.mlo++;
  63                 if (ux.bits.mlo == 0) {
  64                         ux.bits.mhi++;
  65                         if (ux.bits.mhi == 0)
  66                                 ux.bits.exp++;
  67                 }
  68         } else {
  69                 if (ux.bits.mlo == 0) {
  70                         if (ux.bits.mhi == 0)
  71                                 ux.bits.exp--;
  72                         ux.bits.mhi--;
  73                 }
  74                 ux.bits.mlo--;
  75         }
  76         /* raise overflow if ux.value is infinite and x is finite */
  77         if (ux.bits.exp == 0x7fff)
  78                 return x + x;
  79         /* raise underflow if ux.value is subnormal or zero */
  80         if (ux.bits.exp == 0) {
  81                 volatile float z;
  82                 z = x*x + ux.value*ux.value;
  83         }
  84         return ux.value;
  85 }
  86 #endif