1 /*
2  * inspired by glibc-2.0.6/sysdeps/libm-ieee754/s_nextafterf.c
3  *
4  * gcc -O2 -S -DOP=+ gives faddp %st(1),%st
5  * gcc -O2 -S -DOP=* gives fmulp %st(1),%st
6  * gcc -O2 -S -DOP=- gives fsubrp %st(1),%st
7  * gcc -O2 -S -DOP=/ gives fdivrp %st(1),%st
8  */
9 
10 #ifndef OP
11 #define OP *
12 #endif
13 
14 typedef int int32_t __attribute__ ((__mode__ (  __SI__ ))) ;
15 typedef unsigned int u_int32_t __attribute__ ((__mode__ (  __SI__ ))) ;
16 
17 typedef union
18 {
19   float value;
20   u_int32_t word;
21 } ieee_float_shape_type;
22 
__nextafterf(float x,float y)23 float __nextafterf(float x, float y)
24 {
25  int32_t hx,hy,ix,iy;
26 
27  {
28   ieee_float_shape_type gf_u;
29   gf_u.value = x;
30   hx = gf_u.word;
31  }
32  {
33   ieee_float_shape_type gf_u;
34   gf_u.value = y;
35   hy = gf_u.word;
36  }
37  ix = hx&0x7fffffff;
38  iy = hy&0x7fffffff;
39 
40  if ( ix > 0x7f800000 || iy > 0x7f800000 )
41     return x+y;
42  if (x == y) return x;
43  if (ix == 0)
44    {
45     {
46      ieee_float_shape_type sf_u;
47      sf_u.word = (hy&0x80000000) | 1;
48      x = sf_u.value;
49     }
50     y = x*x;
51     if (y == x) return y; else return x;
52    }
53  if (hx >= 0)
54    {
55     if (hx > hy)
56        hx -= 1;
57     else
58        hx += 1;
59    }
60  else
61    {
62     if (hy >= 0 || hx > hy)
63        hx -= 1;
64     else
65        hx += 1;
66    }
67  hy = hx & 0x7f800000;
68  if (hy >= 0x7f800000)
69     return x+x;
70  if (hy < 0x00800000)
71    {
72     y = x OP x;
73     if (y != x)
74       {
75        ieee_float_shape_type sf_u;
76        sf_u.word = hx;
77        y = sf_u.value;
78        return y;
79       }
80    }
81  {
82   ieee_float_shape_type sf_u;
83   sf_u.word = hx;
84   x = sf_u.value;
85  }
86  return x;
87 }
88 
89 
90