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