sysdeps/libm-ieee754/s_nexttoward.c

   1 /* s_nexttoward.c
   2  * Conversion from s_nextafter.c by Ulrich Drepper, Cygnus Support,
   3  * drepper@cygnus.com.
   4  */
   5
   6 /*
   7  * ====================================================
   8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   9  *
  10  * Developed at SunPro, a Sun Microsystems, Inc. business.
  11  * Permission to use, copy, modify, and distribute this
  12  * software is freely granted, provided that this notice
  13  * is preserved.
  14  * ====================================================
  15  */
  16
  17 #if defined(LIBM_SCCS) && !defined(lint)
  18 static char rcsid[] = "$NetBSD: $";
  19 #endif
  20
  21 /* IEEE functions
  22  *      nextafterx(x,y)
  23  *      return the next machine floating-point number of x in the
  24  *      direction toward y.
  25  *   Special cases:
  26  */
  27
  28 #include "math.h"
  29 #include "math_private.h"
  30
  31 #ifdef __STDC__
  32         double __nexttoward(double x, long double y)
  33 #else
  34         double __nexttoward(x,y)
  35         double x;
  36         long double y;
  37 #endif
  38 {
  39         int32_t hx,ix,iy;
  40         u_int32_t lx,hy,ly,esy;
  41
  42         EXTRACT_WORDS(hx,lx,x);
  43         GET_LDOUBLE_WORDS(esy,hy,ly,y);
  44         ix = hx&0x7fffffff;             /* |x| */
  45         iy = esy&0x7fff;                /* |y| */
  46
  47         if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */
  48            ((iy>=0x7fff)&&((hy|ly)|-(hy|ly))!=0))           /* y is nan */
  49            return x+y;
  50         if((long double) x==y) return x;        /* x=y, return x */
  51         if((ix|lx)==0) {                        /* x == 0 */
  52             double x2;
  53             INSERT_WORDS(x,esy&0x8000?0x80000000:0,1);/* return +-minsub */
  54             x2 = x*x;
  55             if(x2==x) return x2; else return x; /* raise underflow flag */
  56         }
  57         if(hx>=0) {                             /* x > 0 */
  58             if (esy>=0x8000||((ix>>20)&0x7ff)>iy
  59                 || (((ix>>20)&0x7ff)==iy
  60                     && (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
  61                         || (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
  62                             && (lx<<11)>ly)))) {        /* x > y, x -= ulp */
  63                 if(lx==0) hx -= 1;
  64                 lx -= 1;
  65             } else {                            /* x < y, x += ulp */
  66                 lx += 1;
  67                 if(lx==0) hx += 1;
  68             }
  69         } else {                                /* x < 0 */
  70             if (esy<0x8000||((ix>>20)&0x7ff)>iy
  71                 || (((ix>>20)&0x7ff)==iy
  72                     && (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
  73                         || (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
  74                             && (lx<<11)>ly))))  {/* x < y, x -= ulp */
  75                 if(lx==0) hx -= 1;
  76                 lx -= 1;
  77             } else {                            /* x > y, x += ulp */
  78                 lx += 1;
  79                 if(lx==0) hx += 1;
  80             }
  81         }
  82         hy = hx&0x7ff00000;
  83         if(hy>=0x7ff00000) return x+x;  /* overflow  */
  84         if(hy<0x00100000) {             /* underflow */
  85             double x2 = x*x;
  86             if(x2!=x) {         /* raise underflow flag */
  87                 INSERT_WORDS(x2,hx,lx);
  88                 return x2;
  89             }
  90         }
  91         INSERT_WORDS(x,hx,lx);
  92         return x;
  93 }
  94 weak_alias (__nexttoward, nexttoward)
  95 #ifdef NO_LONG_DOUBLE
  96 strong_alias (__nexttoward, __nexttowardl)
  97 weak_alias (__nexttoward, nexttowardl)
  98 #endif