overby.me
this repo has no description
1/*
2 * remainderf.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_ilogbf( float x ) __attribute__ ((always_inline));
24static inline int local_ilogbf( float x )
25{
26 union{ float f; uint32_t u;}u = {x};
27
28 u.u &= 0x7fffffff;
29 int32_t exp = u.u >> 23;
30
31 if( __builtin_expect( (uint32_t) exp - 1U >= 254, 0 ) )
32 { // +-denorm, other possibilities: +-0, +-inf, NaN are screened out by caller
33
34 u.u |= 0x3f800000U;
35 u.f -= 1.0f;
36 exp = u.u >> 23;
37
38 return exp - (127+126);
39 }
40
41 return exp - 127;
42}
43
44static inline float local_scalbnf( float x, int n ) __attribute__ ((always_inline));
45static inline float local_scalbnf( float x, int n )
46{
47 //create 2**n in double
48 union{ uint64_t u; double d;} u = { ( (int64_t) n + 1023) << 52 };
49 double d = x;
50
51 //scale by appropriate power of 2
52 d *= u.d;
53
54 //return result
55 return (float) d;
56}
57
58
59
60float remainderf( float x, float y )
61{
62
63 // deal with NaN
64 if( x != x || y != y )
65 return x + y;
66
67 //x = Inf or y = 0, return Invalid per IEEE-754
68 float fabsx = __builtin_fabsf(x);
69 if( fabsx == __builtin_inff() || 0.0f == y )
70 {
71#ifdef ARMLIBM_SET_FLAGS
72 return required_add_float( __builtin_inff(), -__builtin_inff() ); //set invalid
73#else
74 return __builtin_nanf("");
75#endif
76 }
77
78 //handle trivial case
79 float fabsy = __builtin_fabsf(y);
80 if( fabsy == __builtin_inff() || 0.0f == x )
81 {
82#ifdef ARMLIBM_SET_FLAGS
83 required_add_float( fabsx, 0.0f ); //signal underflow (that is, no flag set, but exception occurs if unmasked) if x is denorm
84#endif
85 return x;
86 }
87
88 //we have to work
89 int32_t iquo = 0;
90 int32_t iscx = local_ilogbf(fabsx);
91 int32_t iscy = local_ilogbf(fabsy);
92 int32_t idiff = iscx - iscy;
93 float x1, y1, z;
94 x1 = fabsx;
95 y1 = fabsy;
96 union{ float d[2]; uint32_t u[2];} u = { { x, y } };
97 int32_t sign = u.u[0] ^ u.u[1];
98 sign = sign >> 31;
99
100 if( idiff >= 0 )
101 {
102 int32_t i;
103
104 if( idiff )
105 {
106 y1 = local_scalbnf( y1, -iscy );
107 x1 = local_scalbnf( x1, -iscx );
108 for( i = idiff; i != 0; i-- )
109 {
110 if( x1 >= y1 )
111 {
112 x1 -= y1;
113 iquo += 1;
114 }
115 iquo += iquo;
116 x1 += x1;
117 }
118 x1 = local_scalbnf( x1, iscy );
119 }
120 if( x1 >= fabsy )
121 {
122 x1 -= fabsy;
123 iquo += 1;
124 }
125 }
126
127 if( x1 < 0x1.0p127f )
128 {
129 z = x1 + x1;
130 if( (z > fabsy) || ((z == fabsy) && ((iquo & 1) != 0 )))
131 x1 -= fabsy;
132 }
133 else
134 { // x1 is only this large if y is even bigger, so we can safely divide y by 2
135 float halfy = 0.5f * fabsy;
136 if( (x1 > halfy) || ((x1 == halfy) & ((iquo & 1) != 0 )))
137 x1 -= fabsy;
138 }
139
140 if( x < -0.0f )
141 x1 = -x1;
142
143 return x1;
144}
145