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 #include <sys/cdefs.h>
  16 __RCSID("$NetBSD: s_nexttoward.c,v 1.1 2010/09/15 16:12:05 christos 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 /*
  34  * The nexttoward() function is equivalent to nextafter() function,
  35  * except that the second parameter shall have type long double and
  36  * the functions shall return y converted to the type of the function
  37  * if x equals y.
  38  *
  39  * Special cases: XXX
  40  */
  41 double
  42 nexttoward(double x, long double y)
  43 {
  44         union ieee_ext_u uy;
  45         volatile double t;
  46         int32_t hx, ix;
  47         uint32_t lx;
  48
  49         EXTRACT_WORDS(hx, lx, x);
  50         ix = hx & 0x7fffffff;                   /* |x| */
  51         uy.extu_ld = y;
  52
  53         if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) ||
  54             (uy.extu_exp == 0x7fff &&
  55                 ((uy.extu_frach & ~LDBL_NBIT) | uy.extu_fracl) != 0))
  56                 return x+y;                     /* x or y is nan */
  57
  58         if (x == y)
  59                 return (double)y;               /* x=y, return y */
  60
  61         if (x == 0.0) {
  62                 INSERT_WORDS(x, uy.extu_sign<<31, 1);   /* return +-minsubnormal */
  63                 t = x*x;
  64                 if (t == x)
  65                         return t;
  66                 else
  67                         return x;               /* raise underflow flag */
  68         }
  69
  70         if ((hx > 0.0) ^ (x < y)) {             /* x -= ulp */
  71                 if (lx == 0) hx -= 1;
  72                 lx -= 1;
  73         } else {                                /* x += ulp */
  74                 lx += 1;
  75                 if (lx == 0) hx += 1;
  76         }
  77         ix = hx & 0x7ff00000;
  78         if (ix >= 0x7ff00000) return x+x;       /* overflow  */
  79         if (ix <  0x00100000) {                 /* underflow */
  80                 t = x*x;
  81                 if (t != x) {                   /* raise underflow flag */
  82                         INSERT_WORDS(y, hx, lx);
  83                         return y;
  84                 }
  85         }
  86         INSERT_WORDS(x, hx, lx);
  87
  88         return x;
  89 }