overby.me
this repo has no description
1/*
2 * remainder.c
3 * cLibm
4 *
5 * Written by Jon Okada, started on December 7th, 1992.
6 * Modified by Paul Finlayson (PAF) for MathLib v2.
7 * Modified by A. Sazegari (ali) for MathLib v3.
8 * Modified and ported by Robert A. Murley (ram) for Mac OS X.
9 * Modified for cLibm, fixed a few edge cases, rewrote local_ funcs by Ian Ollmann.
10 *
11 * Copyright 2007 Apple Inc. All rights reserved.
12 *
13 */
14
15 #include <math.h>
16 #include <float.h>
17 #include <stdint.h>
18
19#ifdef ARMLIBM_SET_FLAGS
20 #include "required_arithmetic.h"
21#endif
22
23static inline int local_ilogb( double x ) __attribute__ ((always_inline));
24static inline int local_ilogb( double x )
25{
26 union{ double d; uint64_t u;}u = {x};
27
28 u.u &= 0x7fffffffffffffffULL;
29 int32_t exp = (int32_t) (u.u >> 52);
30
31 if( __builtin_expect( (uint32_t) exp - 1U >= 2046, 0 ) )
32 { // +-denorm, the other possibilities: +0 +-inf, NaN are screend out by caller
33 u.u |= 0x3ff0000000000000ULL;
34 u.d -= 1.0;
35 exp = (int32_t) (u.u >> 52);
36
37 return exp - (1023+1022);
38 }
39
40 return exp - 1023;
41}
42
43static inline double local_scalbn( double x, int n );
44static inline double local_scalbn( double x, int n )
45{
46 union{ uint64_t u; double d;} u;
47
48 uint32_t absn = n >> (8*sizeof(n)-1);
49 absn = (n ^ absn) - absn;
50
51 if( absn > 1022 )
52 {
53 // step = copysign( 1022, n )
54 int signMask = n >> (8 * sizeof( int ) - 1);
55 int step = (1022 ^ signMask) - signMask;
56
57 u.u = ( (int64_t) step + 1023ULL) << 52;
58
59 if( n < 0 )
60 {
61 do
62 {
63 x *= u.d;
64 n -= step;
65 }while( n < -1022 );
66 }
67 else
68 {
69 do
70 {
71 x *= u.d;
72 n -= step;
73 }while( n > 1022 );
74 }
75 }
76
77 //create 2**n in double
78 u.u = ( (int64_t) n + 1023ULL) << 52;
79
80 //scale by appropriate power of 2
81 x *= u.d;
82
83 //return result
84 return x;
85}
86
87
88double remainder( double x, double y )
89{
90
91 // deal with NaN
92 if( x != x || y != y )
93 return x + y;
94
95 //x = Inf or y = 0, return Invalid per IEEE-754
96 double fabsx = __builtin_fabs(x);
97 if( fabsx == __builtin_inf() || 0.0 == y )
98 {
99#ifdef ARMLIBM_SET_FLAGS
100 return required_add_double( __builtin_inf(), -__builtin_inf() ); //set invalid
101#else
102 return __builtin_nan("");
103#endif
104 }
105
106 //handle trivial case
107 double fabsy = __builtin_fabs(y);
108 if( fabsy == __builtin_inf() || 0.0 == x )
109 {
110#ifdef ARMLIBM_SET_FLAGS
111 required_add_double( fabsx, 0.0 ); // signal underflow (that is, no flag set,
112 // but exception occurs if unmasked) if x is denorm
113#endif
114 return x;
115 }
116
117 //we have to work
118 int32_t iquo = 0;
119 int32_t iscx = local_ilogb(fabsx);
120 int32_t iscy = local_ilogb(fabsy);
121 int32_t idiff = iscx - iscy;
122 double x1, y1, z;
123 x1 = fabsx;
124 y1 = fabsy;
125 if( idiff >= 0 )
126 {
127 int32_t i;
128
129 if( idiff )
130 {
131 y1 = local_scalbn( y1, -iscy );
132 x1 = local_scalbn( x1, -iscx );
133 for( i = idiff; i != 0; i-- )
134 {
135 if( x1 >= y1 )
136 {
137 x1 -= y1;
138 iquo += 1;
139 }
140 iquo += iquo;
141 x1 += x1;
142 }
143 x1 = local_scalbn( x1, iscy );
144 }
145 if( x1 >= fabsy )
146 {
147 x1 -= fabsy;
148 iquo += 1;
149 }
150 }
151
152 if( x1 < 0x1.0p1022 )
153 {
154 z = x1 + x1;
155 if( (z > fabsy) || ((z == fabsy) & ((iquo & 1) != 0 )))
156 x1 -= fabsy;
157 }
158 else
159 {
160 // x1 is only this large if y is even bigger, so we can safely divide y by 2
161 double halfy = 0.5*fabsy;
162 if( (x1 > halfy) || ((x1 == halfy) & ((iquo & 1) != 0 )))
163 x1 -= fabsy;
164 }
165
166
167 if( x < 0 )
168 x1 = -x1;
169
170 return x1;
171}
172