1 /*
2 * rredf.c - trigonometric range reduction function
3 *
4 * Copyright (c) 2009-2018, Arm Limited.
5 * SPDX-License-Identifier: MIT
6 */
7
8 /*
9 * This code is intended to be used as the second half of a range
10 * reducer whose first half is an inline function defined in
11 * rredf.h. Each trig function performs range reduction by invoking
12 * that, which handles the quickest and most common cases inline
13 * before handing off to this function for everything else. Thus a
14 * reasonable compromise is struck between speed and space. (I
15 * hope.) In particular, this approach avoids a function call
16 * overhead in the common case.
17 */
18
19 #include "math_private.h"
20
21 #ifdef __cplusplus
22 extern "C" {
23 #endif /* __cplusplus */
24
25 /*
26 * Input values to this function:
27 * - x is the original user input value, unchanged by the
28 * first-tier reducer in the case where it hands over to us.
29 * - q is still the place where the caller expects us to leave the
30 * quadrant code.
31 * - k is the IEEE bit pattern of x (which it would seem a shame to
32 * recompute given that the first-tier reducer already went to
33 * the effort of extracting it from the VFP). FIXME: in softfp,
34 * on the other hand, it's unconscionably wasteful to replicate
35 * this value into a second register and we should change the
36 * prototype!
37 */
__mathlib_rredf2(float x,int * q,unsigned k)38 float __mathlib_rredf2(float x, int *q, unsigned k)
39 {
40 /*
41 * First, weed out infinities and NaNs, and deal with them by
42 * returning a negative q.
43 */
44 if ((k << 1) >= 0xFF000000) {
45 *q = -1;
46 return x;
47 }
48 /*
49 * We do general range reduction by multiplying by 2/pi, and
50 * retaining the bottom two bits of the integer part and an
51 * initial chunk of the fraction below that. The integer bits
52 * are directly output as *q; the fraction is then multiplied
53 * back up by pi/2 before returning it.
54 *
55 * To get this right, we don't have to multiply by the _whole_
56 * of 2/pi right from the most significant bit downwards:
57 * instead we can discard any bit of 2/pi with a place value
58 * high enough that multiplying it by the LSB of x will yield a
59 * place value higher than 2. Thus we can bound the required
60 * work by a reasonably small constant regardless of the size of
61 * x (unlike, for instance, the IEEE remainder operation).
62 *
63 * At the other end, however, we must take more care: it isn't
64 * adequate just to acquire two integer bits and 24 fraction
65 * bits of (2/pi)x, because if a lot of those fraction bits are
66 * zero then we will suffer significance loss. So we must keep
67 * computing fraction bits as far down as 23 bits below the
68 * _highest set fraction bit_.
69 *
70 * The immediate question, therefore, is what the bound on this
71 * end of the job will be. In other words: what is the smallest
72 * difference between an integer multiple of pi/2 and a
73 * representable IEEE single precision number larger than the
74 * maximum size handled by rredf.h?
75 *
76 * The most difficult cases for each exponent can readily be
77 * found by Tim Peters's modular minimisation algorithm, and are
78 * tabulated in mathlib/tests/directed/rredf.tst. The single
79 * worst case is the IEEE single-precision number 0x6F79BE45,
80 * whose numerical value is in the region of 7.7*10^28; when
81 * reduced mod pi/2, it attains the value 0x30DDEEA9, or about
82 * 0.00000000161. The highest set bit of this value is the one
83 * with place value 2^-30; so its lowest is 2^-53. Hence, to be
84 * sure of having enough fraction bits to output at full single
85 * precision, we must be prepared to collect up to 53 bits of
86 * fraction in addition to our two bits of integer part.
87 *
88 * To begin with, this means we must store the value of 2/pi to
89 * a precision of 128+53 = 181 bits. That's six 32-bit words.
90 * (Hardly a chore, unlike the equivalent problem in double
91 * precision!)
92 */
93 {
94 static const unsigned twooverpi[] = {
95 /* We start with a zero word, because that takes up less
96 * space than the array bounds checking and special-case
97 * handling that would have to occur in its absence. */
98 0,
99 /* 2/pi in hex is 0.a2f9836e... */
100 0xa2f9836e, 0x4e441529, 0xfc2757d1,
101 0xf534ddc0, 0xdb629599, 0x3c439041,
102 /* Again, to avoid array bounds overrun, we store a spare
103 * word at the end. And it would be a shame to fill it
104 * with zeroes when we could use more bits of 2/pi... */
105 0xfe5163ab
106 };
107
108 /*
109 * Multiprecision multiplication of this nature is more
110 * readily done in integers than in VFP, since we can use
111 * UMULL (on CPUs that support it) to multiply 32 by 32 bits
112 * at a time whereas the VFP would only be able to do 12x12
113 * without losing accuracy.
114 *
115 * So extract the mantissa of the input number as a 32-bit
116 * integer.
117 */
118 unsigned mantissa = 0x80000000 | (k << 8);
119
120 /*
121 * Now work out which part of our stored value of 2/pi we're
122 * supposed to be multiplying by.
123 *
124 * Let the IEEE exponent field of x be e. With its bias
125 * removed, (e-127) is the index of the set bit at the top
126 * of 'mantissa' (i.e. that set bit has real place value
127 * 2^(e-127)). So the lowest set bit in 'mantissa', 23 bits
128 * further down, must have place value 2^(e-150).
129 *
130 * We begin taking an interest in the value of 2/pi at the
131 * bit which multiplies by _that_ to give something with
132 * place value at most 2. In other words, the highest bit of
133 * 2/pi we're interested in is the one with place value
134 * 2/(2^(e-150)) = 2^(151-e).
135 *
136 * The bit at the top of the first (zero) word of the above
137 * array has place value 2^31. Hence, the bit we want to put
138 * at the top of the first word we extract from that array
139 * is the one at bit index n, where 31-n = 151-e and hence
140 * n=e-120.
141 */
142 int topbitindex = ((k >> 23) & 0xFF) - 120;
143 int wordindex = topbitindex >> 5;
144 int shiftup = topbitindex & 31;
145 int shiftdown = 32 - shiftup;
146 unsigned word1, word2, word3;
147 if (shiftup) {
148 word1 = (twooverpi[wordindex] << shiftup) | (twooverpi[wordindex+1] >> shiftdown);
149 word2 = (twooverpi[wordindex+1] << shiftup) | (twooverpi[wordindex+2] >> shiftdown);
150 word3 = (twooverpi[wordindex+2] << shiftup) | (twooverpi[wordindex+3] >> shiftdown);
151 } else {
152 word1 = twooverpi[wordindex];
153 word2 = twooverpi[wordindex+1];
154 word3 = twooverpi[wordindex+2];
155 }
156
157 /*
158 * Do the multiplications, and add them together.
159 */
160 unsigned long long mult1 = (unsigned long long)word1 * mantissa;
161 unsigned long long mult2 = (unsigned long long)word2 * mantissa;
162 unsigned long long mult3 = (unsigned long long)word3 * mantissa;
163
164 unsigned /* bottom3 = (unsigned)mult3, */ top3 = (unsigned)(mult3 >> 32);
165 unsigned bottom2 = (unsigned)mult2, top2 = (unsigned)(mult2 >> 32);
166 unsigned bottom1 = (unsigned)mult1, top1 = (unsigned)(mult1 >> 32);
167
168 unsigned out3, out2, out1, carry;
169
170 out3 = top3 + bottom2; carry = (out3 < top3);
171 out2 = top2 + bottom1 + carry; carry = carry ? (out2 <= top2) : (out2 < top2);
172 out1 = top1 + carry;
173
174 /*
175 * The two words we multiplied to get mult1 had their top
176 * bits at (respectively) place values 2^(151-e) and
177 * 2^(e-127). The value of those two bits multiplied
178 * together will have ended up in bit 62 (the
179 * topmost-but-one bit) of mult1, i.e. bit 30 of out1.
180 * Hence, that bit has place value 2^(151-e+e-127) = 2^24.
181 * So the integer value that we want to output as q,
182 * consisting of the bits with place values 2^1 and 2^0,
183 * must be 23 and 24 bits below that, i.e. in bits 7 and 6
184 * of out1.
185 *
186 * Or, at least, it will be once we add 1/2, to round to the
187 * _nearest_ multiple of pi/2 rather than the next one down.
188 */
189 *q = (out1 + (1<<5)) >> 6;
190
191 /*
192 * Now we construct the output fraction, which is most
193 * simply done in the VFP. We just extract three consecutive
194 * bit strings from our chunk of binary data, convert them
195 * to integers, equip each with an appropriate FP exponent,
196 * add them together, and (don't forget) multiply back up by
197 * pi/2. That way we don't have to work out ourselves where
198 * the highest fraction bit ended up.
199 *
200 * Since our displacement from the nearest multiple of pi/2
201 * can be positive or negative, the topmost of these three
202 * values must be arranged with its 2^-1 bit at the very top
203 * of the word, and then treated as a _signed_ integer.
204 */
205 {
206 int i1 = (out1 << 26) | ((out2 >> 19) << 13);
207 unsigned i2 = out2 << 13;
208 unsigned i3 = out3;
209 float f1 = i1, f2 = i2 * (1.0f/524288.0f), f3 = i3 * (1.0f/524288.0f/524288.0f);
210
211 /*
212 * Now f1+f2+f3 is a representation, potentially to
213 * twice double precision, of 2^32 times ((2/pi)*x minus
214 * some integer). So our remaining job is to multiply
215 * back down by (pi/2)*2^-32, and convert back to one
216 * single-precision output number.
217 */
218
219 /* Normalise to a prec-and-a-half representation... */
220 float ftop = CLEARBOTTOMHALF(f1+f2+f3), fbot = f3-((ftop-f1)-f2);
221
222 /* ... and multiply by a prec-and-a-half value of (pi/2)*2^-32. */
223 float ret = (ftop * 0x1.92p-32F) + (ftop * 0x1.fb5444p-44F + fbot * 0x1.921fb6p-32F);
224
225 /* Just before we return, take the input sign into account. */
226 if (k & 0x80000000) {
227 *q = 0x10000000 - *q;
228 ret = -ret;
229 }
230 return ret;
231 }
232 }
233 }
234
235 #ifdef __cplusplus
236 } /* end of extern "C" */
237 #endif /* __cplusplus */
238
239 /* end of rredf.c */
240