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