overby.me
this repo has no description
1/*
2 * remquof.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 remquof( float x, float y, int *quo )
61{
62
63 *quo = 0; // Init quotient result
64
65 // deal with NaN
66 if( x != x || y != y )
67 return x + y;
68
69 //x = Inf or y = 0, return Invalid per IEEE-754
70 float fabsx = __builtin_fabsf(x);
71 if( fabsx == __builtin_inff() || 0.0f == y )
72 {
73#ifdef ARMLIBM_SET_FLAGS
74 return required_add_float( __builtin_inff(), -__builtin_inff() ); //set invalid
75#else
76 return __builtin_nanf("");
77#endif
78 }
79
80 //handle trivial case
81 float fabsy = __builtin_fabsf(y);
82 if( fabsy == __builtin_inff() || 0.0f == x )
83 {
84#ifdef ARMLIBM_SET_FLAGS
85 required_add_float( fabsx, 0.0f ); //signal underflow (that is, no flag set, but exception occurs if unmasked) if x is denorm
86#endif
87 return x;
88 }
89
90 //we have to work
91 int32_t iquo = 0;
92 int32_t iscx = local_ilogbf(fabsx);
93 int32_t iscy = local_ilogbf(fabsy);
94 int32_t idiff = iscx - iscy;
95 float x1, y1, z;
96 x1 = fabsx;
97 y1 = fabsy;
98 union{ float d[2]; uint32_t u[2];} u = { { x, y } };
99 int32_t sign = u.u[0] ^ u.u[1];
100 sign = sign >> 31;
101
102 if( idiff >= 0 )
103 {
104 int32_t i;
105
106 if( idiff )
107 {
108 y1 = local_scalbnf( y1, -iscy );
109 x1 = local_scalbnf( x1, -iscx );
110 for( i = idiff; i != 0; i-- )
111 {
112 if( x1 >= y1 )
113 {
114 x1 -= y1;
115 iquo += 1;
116 }
117 iquo += iquo;
118 x1 += x1;
119 }
120 x1 = local_scalbnf( x1, iscy );
121 }
122 if( x1 >= fabsy )
123 {
124 x1 -= fabsy;
125 iquo += 1;
126 }
127 }
128
129 if( x1 < 0x1.0p127f )
130 {
131 z = x1 + x1;
132 if( (z > fabsy) || ((z == fabsy) && ((iquo & 1) != 0 )))
133 {
134 x1 -= fabsy;
135 iquo += 1;
136 }
137 }
138 else
139 { // x1 is only this large if y is even bigger, so we can safely divide y by 2
140 float halfy = 0.5f * fabsy;
141 if( (x1 > halfy) || ((x1 == halfy) & ((iquo & 1) != 0 )))
142 {
143 x1 -= fabsy;
144 iquo += 1;
145 }
146 }
147
148 if( x < -0.0f )
149 x1 = -x1;
150 iquo &= 0x7f;
151
152 *quo = (iquo ^ sign) - sign;
153
154 return x1;
155}
156