lib/libm/src/s_nexttoward.c

   1 /*      $NetBSD: s_nexttoward.c,v 1.2 2013/08/21 13:03:56 martin Exp $  */
   2
   3 /* @(#)s_nextafter.c 5.1 93/09/24 */
   4 /*
   5  * ====================================================
   6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   7  *
   8  * Developed at SunPro, a Sun Microsystems, Inc. business.
   9  * Permission to use, copy, modify, and distribute this
  10  * software is freely granted, provided that this notice
  11  * is preserved.
  12  * ====================================================
  13  */
  14
  15 #include <sys/cdefs.h>
  16 __RCSID("$NetBSD: s_nexttoward.c,v 1.2 2013/08/21 13:03:56 martin Exp $");
  17
  18 /*
  19  * We assume that a long double has a 15-bit exponent.  On systems
  20  * where long double is the same as double, nexttoward() is an alias
  21  * for nextafter(), so we don't use this routine.
  22  */
  23 #include <float.h>
  24
  25 #include <machine/ieee.h>
  26 #include "math.h"
  27 #include "math_private.h"
  28
  29 #if LDBL_MAX_EXP != 0x4000
  30 #error "Unsupported long double format"
  31 #endif
  32
  33 #ifdef LDBL_IMPLICIT_NBIT
  34 #define LDBL_NBIT       0
  35 #endif
  36
  37 /*
  38  * The nexttoward() function is equivalent to nextafter() function,
  39  * except that the second parameter shall have type long double and
  40  * the functions shall return y converted to the type of the function
  41  * if x equals y.
  42  *
  43  * Special cases: XXX
  44  */
  45 double
  46 nexttoward(double x, long double y)
  47 {
  48         union ieee_ext_u uy;
  49         volatile double t;
  50         int32_t hx, ix;
  51         uint32_t lx;
  52
  53         EXTRACT_WORDS(hx, lx, x);
  54         ix = hx & 0x7fffffff;                   /* |x| */
  55         uy.extu_ld = y;
  56
  57         if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) ||
  58             (uy.extu_exp == 0x7fff &&
  59                 ((uy.extu_frach & ~LDBL_NBIT) | uy.extu_fracl) != 0))
  60                 return x+y;                     /* x or y is nan */
  61
  62         if (x == y)
  63                 return (double)y;               /* x=y, return y */
  64
  65         if (x == 0.0) {
  66                 INSERT_WORDS(x, uy.extu_sign<<31, 1);   /* return +-minsubnormal */
  67                 t = x*x;
  68                 if (t == x)
  69                         return t;
  70                 else
  71                         return x;               /* raise underflow flag */
  72         }
  73
  74         if ((hx > 0.0) ^ (x < y)) {             /* x -= ulp */
  75                 if (lx == 0) hx -= 1;
  76                 lx -= 1;
  77         } else {                                /* x += ulp */
  78                 lx += 1;
  79                 if (lx == 0) hx += 1;
  80         }
  81         ix = hx & 0x7ff00000;
  82         if (ix >= 0x7ff00000) return x+x;       /* overflow  */
  83         if (ix <  0x00100000) {                 /* underflow */
  84                 t = x*x;
  85                 if (t != x) {                   /* raise underflow flag */
  86                         INSERT_WORDS(y, hx, lx);
  87                         return y;
  88                 }
  89         }
  90         INSERT_WORDS(x, hx, lx);
  91
  92         return x;
  93 }