1 /* @(#)s_nextafter.c 5.1 93/09/24 */
2 /*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13 /* IEEE functions
14 * nexttoward(x,y)
15 * return the next machine floating-point number of x in the
16 * direction toward y.
17 * Special cases:
18 */
19
20 #include <math.h>
21 #include <float.h>
22
23 #include "math_private.h"
24
25 double
nexttoward(double x,long double y)26 nexttoward(double x, long double y)
27 {
28 int32_t hx,ix;
29 int64_t hy,iy;
30 u_int32_t lx;
31 u_int64_t ly;
32
33 EXTRACT_WORDS(hx,lx,x);
34 GET_LDOUBLE_WORDS64(hy,ly,y);
35 ix = hx&0x7fffffff; /* |x| */
36 iy = hy&0x7fffffffffffffffLL; /* |y| */
37
38 if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */
39 ((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0))
40 /* y is nan */
41 return x+y;
42 if((long double) x==y) return y; /* x=y, return y */
43 if((ix|lx)==0) { /* x == 0 */
44 volatile double u;
45 INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */
46 u = x;
47 u = u * u; /* raise underflow flag */
48 return x;
49 }
50 if(hx>=0) { /* x > 0 */
51 if (hy<0||(ix>>20)>(iy>>48)-0x3c00
52 || ((ix>>20)==(iy>>48)-0x3c00
53 && (((((int64_t)hx)<<28)|(lx>>4))>(hy&0x0000ffffffffffffLL)
54 || (((((int64_t)hx)<<28)|(lx>>4))==(hy&0x0000ffffffffffffLL)
55 && (lx&0xf)>(ly>>60))))) { /* x > y, x -= ulp */
56 if(lx==0) hx -= 1;
57 lx -= 1;
58 } else { /* x < y, x += ulp */
59 lx += 1;
60 if(lx==0) hx += 1;
61 }
62 } else { /* x < 0 */
63 if (hy>=0||(ix>>20)>(iy>>48)-0x3c00
64 || ((ix>>20)==(iy>>48)-0x3c00
65 && (((((int64_t)hx)<<28)|(lx>>4))>(hy&0x0000ffffffffffffLL)
66 || (((((int64_t)hx)<<28)|(lx>>4))==(hy&0x0000ffffffffffffLL)
67 && (lx&0xf)>(ly>>60))))) { /* x < y, x -= ulp */
68 if(lx==0) hx -= 1;
69 lx -= 1;
70 } else { /* x > y, x += ulp */
71 lx += 1;
72 if(lx==0) hx += 1;
73 }
74 }
75 hy = hx&0x7ff00000;
76 if(hy>=0x7ff00000) {
77 x = x+x; /* overflow */
78 return x;
79 }
80 if(hy<0x00100000) {
81 volatile double u = x*x; /* underflow */
82 }
83 INSERT_WORDS(x,hx,lx);
84 return x;
85 }
86