overby.me
this repo has no description
1/*
2 * nextafter.c
3 * cLibm
4 *
5 * Created by Ian Ollmann on 6/14/07.
6 * Copyright 2007 Apple Inc. All rights reserved.
7 *
8 */
9
10#include <math.h>
11#include <stdint.h>
12
13#ifdef ARMLIBM_SET_FLAGS
14#include <fenv.h>
15#include "required_arithmetic.h"
16#pragma STDC FENV_ACCESS ON
17
18double nextafter( double x, double y )
19{
20 union{ double d; uint64_t u;} ux = { x };
21 uint64_t step = 1;
22
23 if( y != y || x != x)
24 return x + y;
25
26 if( y <= x ) // a subtraction here would risk Inf-Inf
27 {
28 if( y == x )
29 return y; //make sure sign of 0 is correct
30
31 step = -1LL;
32 }
33
34 //correct for the sign
35 int64_t signMask = (int64_t) ux.u >> 63;
36 step = (step ^ signMask) - signMask;
37
38 uint64_t absux = ux.u & 0x7fffffffffffffffULL;
39
40 if( absux - 0x0010000000000001ULL >= 0x7fefffffffffffffULL - 0x0010000000000001ULL )
41 { //0, Inf, max finite, min normal, denorm
42 //Nan is handled above and won't occur here
43 //Inf can just fall through. We are only here if y!=x.
44
45 if( absux == 0ULL ) // zero
46 {
47 ux.d = y;
48 required_multiply_double( 0x1.0p-1000, 0x1.0p-1000 );
49 ux.u = (ux.u & 0x8000000000000000ULL) + 1U;
50 return ux.d;
51 }
52 else if( absux == 0x7fefffffffffffffULL ) // max finite
53 {
54 ux.u += step;
55
56 //if infinity is the result, set some flags
57 if( 0 == (ux.u & 2ULL) )
58 {
59 required_add_double( x, 1.0 ); //set inexact
60 required_multiply_double( x, x ); //set overflow
61 }
62
63 return ux.d;
64 }
65
66 ux.u += step;
67 if( 0ULL == (ux.u & 0x7ff0000000000000ULL))
68 required_multiply_double( 0x1.0p-1000, 0x1.0p-1000 );
69
70 return ux.d;
71 }
72
73 ux.u += step;
74
75 return ux.d;
76}
77
78#else
79
80double nextafter( double x, double y )
81{
82 union{ double d; uint64_t u;} ux = { x };
83 uint64_t step = 1;
84
85 if( y != y || x != x)
86 return x + y;
87
88 if( y <= x ) // a subtraction here would risk Inf-Inf
89 {
90 if( y == x )
91 return y; //make sure sign of 0 is correct
92
93 step = -1LL;
94 }
95
96 //correct for the sign
97 int64_t signMask = (int64_t) ux.u >> 63;
98 step = (step ^ signMask) - signMask;
99
100 uint64_t absux = ux.u & 0x7fffffffffffffffULL;
101
102 if( absux == 0ULL ) // zero
103 {
104 ux.d = y;
105 ux.u = (ux.u & 0x8000000000000000ULL) + 1U;
106 return ux.d;
107 }
108
109 ux.u += step;
110 return ux.d;
111}
112
113#endif // ARMLIBM_SET_FLAGS