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 #ifdef _NEED_FLOAT64
21
22 __float64
nexttoward(__float64 x,long double y)23 nexttoward(__float64 x, long double y)
24 {
25 int32_t hx,ix,iy;
26 u_int32_t lx,hy,ly,esy;
27
28 EXTRACT_WORDS(hx,lx,x);
29 GET_LDOUBLE_WORDS(esy,hy,ly,y);
30 ix = hx&0x7fffffff; /* |x| */
31 iy = esy&0x7fff; /* |y| */
32 hy &= 0x7fffffff; /* mask off leading 1 */
33
34 if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */
35 ((iy>=0x7fff)&&(hy|ly)!=0)) /* y is nan */
36 return (long double)x+y;
37 if((long double) x==y) return y; /* x=y, return y */
38 if((ix|lx)==0) { /* x == 0 */
39 INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
40 force_eval_float64(opt_barrier_float64(x)*x);
41 return x;
42 }
43 if(hx>=0) { /* x > 0 */
44 if ((long double) x > y) { /* x > y, x -= ulp */
45 if(lx==0) hx -= 1;
46 lx -= 1;
47 } else { /* x < y, x += ulp */
48 lx += 1;
49 if(lx==0) hx += 1;
50 }
51 } else { /* x < 0 */
52 if ((long double) x < y) { /* x < y, x -= ulp */
53 if(lx==0) hx -= 1;
54 lx -= 1;
55 } else { /* x > y, x += ulp */
56 lx += 1;
57 if(lx==0) hx += 1;
58 }
59 }
60 hy = hx&0x7ff00000;
61 if(hy>=0x7ff00000)
62 return __math_oflow(hx<0);
63 INSERT_WORDS(x,hx,lx);
64 if(hy<0x00100000)
65 return __math_denorm(x);
66 return x;
67 }
68
69 _MATH_ALIAS_d_dl(nexttoward)
70
71 #endif /* _NEED_FLOAT64 */
72
