lib/libm/src/s_nexttoward.c

   1 /*      $NetBSD: s_nexttoward.c,v 1.1 2010/09/15 16:12:05 christos 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 /*
  16  * We assume that a long double has a 15-bit exponent.  On systems
  17  * where long double is the same as double, nexttoward() is an alias
  18  * for nextafter(), so we don't use this routine.
  19  */
  20 #include <float.h>
  21
  22 #include <machine/ieee.h>
  23 #include <math.h>
  24 #include "math_private.h"
  25
  26 #if LDBL_MAX_EXP != 0x4000
  27 #error "Unsupported long double format"
  28 #endif
  29
  30 /*
  31  * The nexttoward() function is equivalent to nextafter() function,
  32  * except that the second parameter shall have type long double and
  33  * the functions shall return y converted to the type of the function
  34  * if x equals y.
  35  *
  36  * Special cases: XXX
  37  */
  38 double
  39 nexttoward(double x, long double y)
  40 {
  41         union ieee_ext_u uy;
  42         volatile double t;
  43         int32_t hx, ix;
  44         uint32_t lx;
  45
  46         EXTRACT_WORDS(hx, lx, x);
  47         ix = hx & 0x7fffffff;                   /* |x| */
  48         uy.extu_ld = y;
  49
  50         if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) ||
  51             (uy.extu_exp == 0x7fff &&
  52                 ((uy.extu_frach & ~LDBL_NBIT) | uy.extu_fracl) != 0))
  53                 return x+y;                     /* x or y is nan */
  54
  55         if (x == y)
  56                 return (double)y;               /* x=y, return y */
  57
  58         if (x == 0.0) {
  59                 INSERT_WORDS(x, uy.extu_sign<<31, 1);   /* return +-minsubnormal */
  60                 t = x*x;
  61                 if (t == x)
  62                         return t;
  63                 else
  64                         return x;               /* raise underflow flag */
  65         }
  66
  67         if ((hx > 0.0) ^ (x < y)) {             /* x -= ulp */
  68                 if (lx == 0) hx -= 1;
  69                 lx -= 1;
  70         } else {                                /* x += ulp */
  71                 lx += 1;
  72                 if (lx == 0) hx += 1;
  73         }
  74         ix = hx & 0x7ff00000;
  75         if (ix >= 0x7ff00000) return x+x;       /* overflow  */
  76         if (ix <  0x00100000) {                 /* underflow */
  77                 t = x*x;
  78                 if (t != x) {                   /* raise underflow flag */
  79                         INSERT_WORDS(y, hx, lx);
  80                         return y;
  81                 }
  82         }
  83         INSERT_WORDS(x, hx, lx);
  84
  85         return x;
  86 }