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
14typedef int int32_t __attribute__ ((__mode__ (  __SI__ ))) ;
15typedef unsigned int u_int32_t __attribute__ ((__mode__ (  __SI__ ))) ;
16
17typedef union
18{
19  float value;
20  u_int32_t word;
21} ieee_float_shape_type;
22
23float __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