lib/libm/src/s_nextafterl.c

   1 /*      $NetBSD: s_nextafterl.c,v 1.2 2010/09/17 20:39:39 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_nextafterl.c,v 1.2 2010/09/17 20:39:39 christos Exp $");
  17
  18 #include <float.h>
  19 #include <math.h>
  20 #include <machine/ieee.h>
  21
  22 #ifdef EXT_EXP_INFNAN
  23 #if LDBL_MAX_EXP != 0x4000
  24 #error "Unsupported long double format"
  25 #endif
  26
  27 /*
  28  * IEEE functions
  29  *      nextafterl(x,y)
  30  *      return the next machine floating-point number of x in the
  31  *      direction toward y.
  32  *   Special cases:
  33  *      If x == y, y shall be returned
  34  *      If x or y is NaN, a NaN shall be returned
  35  */
  36 long double
  37 nextafterl(long double x, long double y)
  38 {
  39         volatile long double t;
  40         union ieee_ext_u ux, uy;
  41
  42         ux.extu_ld = x;
  43         uy.extu_ld = y;
  44
  45         if ((ux.extu_exp == EXT_EXP_NAN &&
  46                 ((ux.extu_frach &~ LDBL_NBIT)|ux.extu_fracl) != 0) ||
  47             (uy.extu_exp == EXT_EXP_NAN &&
  48                 ((uy.extu_frach &~ LDBL_NBIT)|uy.extu_fracl) != 0))
  49                 return x+y;                     /* x or y is nan */
  50
  51         if (x == y) return y;                   /* x=y, return y */
  52
  53         if (x == 0.0) {
  54                 ux.extu_frach = 0;              /* return +-minsubnormal */
  55                 ux.extu_fracl = 1;
  56                 ux.extu_sign = uy.extu_sign;
  57                 t = ux.extu_ld * ux.extu_ld;
  58                 if (t == ux.extu_ld)
  59                         return t;
  60                 else
  61                         return ux.extu_ld;      /* raise underflow flag */
  62         }
  63
  64         if ((x>0.0) ^ (x<y)) {                  /* x -= ulp */
  65                 if (ux.extu_fracl == 0) {
  66                         if ((ux.extu_frach & ~LDBL_NBIT) == 0)
  67                                 ux.extu_exp -= 1;
  68                         ux.extu_frach = (ux.extu_frach - 1) |
  69                                         (ux.extu_frach & LDBL_NBIT);
  70                 }
  71                 ux.extu_fracl -= 1;
  72         } else {                                /* x += ulp */
  73                 ux.extu_fracl += 1;
  74                 if (ux.extu_fracl == 0) {
  75                         ux.extu_frach = (ux.extu_frach + 1) |
  76                                         (ux.extu_frach & LDBL_NBIT);
  77                         if ((ux.extu_frach & ~LDBL_NBIT) == 0)
  78                                 ux.extu_exp += 1;
  79                 }
  80         }
  81
  82         if (ux.extu_exp == EXT_EXP_INF)
  83                 return x+x;                     /* overflow  */
  84
  85         if (ux.extu_exp == 0) {                 /* underflow */
  86                 mask_nbit_l(ux);
  87                 t = ux.extu_ld * ux.extu_ld;
  88                 if (t != ux.extu_ld)            /* raise underflow flag */
  89                         return ux.extu_ld;
  90         }
  91
  92         return ux.extu_ld;
  93 }
  94 #endif