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 #include <stdio.h>
14
15 float
nexttowardf(float x,long double y)16 nexttowardf(float x, long double y)
17 {
18 int32_t hx,ix;
19 int64_t hy,iy;
20 u_int64_t ly;
21
22 GET_FLOAT_WORD(hx,x);
23 GET_LDOUBLE_WORDS64(hy,ly,y);
24 ix = hx&0x7fffffff; /* |x| */
25 iy = hy&0x7fffffffffffffffLL; /* |y| */
26
27 if(ix>0x7f800000) { /* x is nan */
28 force_eval_long_double(y+y);
29 return x + x;
30 }
31 if((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0) { /* y is nan */
32 return (float) (y + y);
33 }
34 if((long double) x==y) return y; /* x=y, return y */
35 if(ix==0) { /* x == 0 */
36 SET_FLOAT_WORD(x,(u_int32_t)((hy>>32)&0x80000000)|1);/* return +-minsub*/
37 force_eval_float(x*x);
38 return x;
39 }
40 if(hx>=0) { /* x > 0 */
41 if((long double) x > y) { /* x > y, x -= ulp */
42 hx -= 1;
43 } else { /* x < y, x += ulp */
44 hx += 1;
45 }
46 } else { /* x < 0 */
47 if((long double) x < y) { /* x < y, x -= ulp */
48 hx -= 1;
49 } else { /* x > y, x += ulp */
50 hx += 1;
51 }
52 }
53 ix = hx&0x7f800000;
54 if(ix>=0x7f800000)
55 return __math_oflowf(hx<0);
56 SET_FLOAT_WORD(x,hx);
57 if(ix<0x00800000)
58 return __math_denormf(x);
59 return x;
60 }
61
