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