1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
/* s_nextafterf.c -- float version of s_nextafter.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "math.h"
#include "math_private.h"
#ifndef math_opt_barrier
# define math_opt_barrier(x) ({ __typeof (x) __x = x; __asm ("" : "+m" (__x)); __x; })
# define math_force_eval(x) __asm __volatile ("" : : "m" (x))
#endif
float nextafterf(float x, float y)
{
int32_t hx, hy, ix, iy;
GET_FLOAT_WORD(hx, x);
GET_FLOAT_WORD(hy, y);
ix = hx & 0x7fffffff; /* |x| */
iy = hy & 0x7fffffff; /* |y| */
/* x is nan or y is nan? */
if ((ix > 0x7f800000) || (iy > 0x7f800000))
return x + y;
if (x == y)
return y;
if (ix == 0) { /* x == 0? */
float u;
/* return +-minsubnormal */
SET_FLOAT_WORD(x, (hy & 0x80000000) | 1);
u = math_opt_barrier(x);
u = u * u;
math_force_eval(u); /* raise underflow flag */
return x;
}
if (hx >= 0) { /* x > 0 */
if (hx > hy) { /* x > y: x -= ulp */
hx -= 1;
} else { /* x < y: x += ulp */
hx += 1;
}
} else { /* x < 0 */
if (hy >= 0 || hx > hy) { /* x < y: x -= ulp */
hx -= 1;
} else { /* x > y: x += ulp */
hx += 1;
}
}
hy = hx & 0x7f800000;
if (hy >= 0x7f800000) {
x = x + x; /* overflow */
//?? if (FLT_EVAL_METHOD != 0)
// asm ("" : "+m"(x));
return x; /* overflow */
}
if (hy < 0x00800000) {
float u = x * x; /* underflow */
math_force_eval(u); /* raise underflow flag */
}
SET_FLOAT_WORD(x, hx);
return x;
}
#if 0
/* "testprog N a b"
* calculates a = nextafterf(a, b) and prints a as float
* and as raw bytes; repeats it N times.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(int argc, char **argv)
{
int cnt, i;
float a, b;
cnt = atoi(argv[1]);
a = strtod(argv[2], NULL);
b = strtod(argv[3], NULL);
while (cnt-- > 0) {
for (i = 0; i < sizeof(a); i++) {
unsigned char c = ((char*)(&a))[i];
printf("%x%x", (c >> 4), (c & 0xf));
}
printf(" %f\n", a);
a = nextafterf(a, b);
}
return 0;
}
#endif
|
