?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Minor cleanups throughout.
[glibc/history.git] / soft-fp / op-2.h
blob5c9bce4c1792ee760c8d1e4f1ef52b524563d8ee
1 /* Software floating-point emulation.
2 Basic two-word fraction declaration and manipulation.
3 Copyright (C) 1997,1998,1999,2006 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
5 Contributed by Richard Henderson (rth@cygnus.com),
6 Jakub Jelinek (jj@ultra.linux.cz),
7 David S. Miller (davem@redhat.com) and
8 Peter Maydell (pmaydell@chiark.greenend.org.uk).
9
10 The GNU C Library is free software; you can redistribute it and/or
11 modify it under the terms of the GNU Lesser General Public
12 License as published by the Free Software Foundation; either
13 version 2.1 of the License, or (at your option) any later version.
14
15 In addition to the permissions in the GNU Lesser General Public
16 License, the Free Software Foundation gives you unlimited
17 permission to link the compiled version of this file into
18 combinations with other programs, and to distribute those
19 combinations without any restriction coming from the use of this
20 file. (The Lesser General Public License restrictions do apply in
21 other respects; for example, they cover modification of the file,
22 and distribution when not linked into a combine executable.)
23
24 The GNU C Library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
28
29 You should have received a copy of the GNU Lesser General Public
30 License along with the GNU C Library; if not, write to the Free
31 Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
32 MA 02110-1301, USA. */
33
34 #define _FP_FRAC_DECL_2(X) _FP_W_TYPE X##_f0, X##_f1
35 #define _FP_FRAC_COPY_2(D,S) (D##_f0 = S##_f0, D##_f1 = S##_f1)
36 #define _FP_FRAC_SET_2(X,I) __FP_FRAC_SET_2(X, I)
37 #define _FP_FRAC_HIGH_2(X) (X##_f1)
38 #define _FP_FRAC_LOW_2(X) (X##_f0)
39 #define _FP_FRAC_WORD_2(X,w) (X##_f##w)
40
41 #define _FP_FRAC_SLL_2(X,N) \
42 (void)(((N) < _FP_W_TYPE_SIZE) \
43 ? ({ \
44 if (__builtin_constant_p(N) && (N) == 1) \
45 { \
46 X##_f1 = X##_f1 + X##_f1 + (((_FP_WS_TYPE)(X##_f0)) < 0); \
47 X##_f0 += X##_f0; \
48 } \
49 else \
50 { \
51 X##_f1 = X##_f1 << (N) | X##_f0 >> (_FP_W_TYPE_SIZE - (N)); \
52 X##_f0 <<= (N); \
53 } \
54 0; \
55 }) \
56 : ({ \
57 X##_f1 = X##_f0 << ((N) - _FP_W_TYPE_SIZE); \
58 X##_f0 = 0; \
59 }))
60
61
62 #define _FP_FRAC_SRL_2(X,N) \
63 (void)(((N) < _FP_W_TYPE_SIZE) \
64 ? ({ \
65 X##_f0 = X##_f0 >> (N) | X##_f1 << (_FP_W_TYPE_SIZE - (N)); \
66 X##_f1 >>= (N); \
67 }) \
68 : ({ \
69 X##_f0 = X##_f1 >> ((N) - _FP_W_TYPE_SIZE); \
70 X##_f1 = 0; \
71 }))
72
73 /* Right shift with sticky-lsb. */
74 #define _FP_FRAC_SRST_2(X,S, N,sz) \
75 (void)(((N) < _FP_W_TYPE_SIZE) \
76 ? ({ \
77 S = (__builtin_constant_p(N) && (N) == 1 \
78 ? X##_f0 & 1 \
79 : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0); \
80 X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N)); \
81 X##_f1 >>= (N); \
82 }) \
83 : ({ \
84 S = ((((N) == _FP_W_TYPE_SIZE \
85 ? 0 \
86 : (X##_f1 << (2*_FP_W_TYPE_SIZE - (N)))) \
87 | X##_f0) != 0); \
88 X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE)); \
89 X##_f1 = 0; \
90 }))
91
92 #define _FP_FRAC_SRS_2(X,N,sz) \
93 (void)(((N) < _FP_W_TYPE_SIZE) \
94 ? ({ \
95 X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N) | \
96 (__builtin_constant_p(N) && (N) == 1 \
97 ? X##_f0 & 1 \
98 : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0)); \
99 X##_f1 >>= (N); \
100 }) \
101 : ({ \
102 X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE) | \
103 ((((N) == _FP_W_TYPE_SIZE \
104 ? 0 \
105 : (X##_f1 << (2*_FP_W_TYPE_SIZE - (N)))) \
106 | X##_f0) != 0)); \
107 X##_f1 = 0; \
108 }))
109
110 #define _FP_FRAC_ADDI_2(X,I) \
111 __FP_FRAC_ADDI_2(X##_f1, X##_f0, I)
112
113 #define _FP_FRAC_ADD_2(R,X,Y) \
114 __FP_FRAC_ADD_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
115
116 #define _FP_FRAC_SUB_2(R,X,Y) \
117 __FP_FRAC_SUB_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
118
119 #define _FP_FRAC_DEC_2(X,Y) \
120 __FP_FRAC_DEC_2(X##_f1, X##_f0, Y##_f1, Y##_f0)
121
122 #define _FP_FRAC_CLZ_2(R,X) \
123 do { \
124 if (X##_f1) \
125 __FP_CLZ(R,X##_f1); \
126 else \
127 { \
128 __FP_CLZ(R,X##_f0); \
129 R += _FP_W_TYPE_SIZE; \
130 } \
131 } while(0)
132
133 /* Predicates */
134 #define _FP_FRAC_NEGP_2(X) ((_FP_WS_TYPE)X##_f1 < 0)
135 #define _FP_FRAC_ZEROP_2(X) ((X##_f1 | X##_f0) == 0)
136 #define _FP_FRAC_OVERP_2(fs,X) (_FP_FRAC_HIGH_##fs(X) & _FP_OVERFLOW_##fs)
137 #define _FP_FRAC_CLEAR_OVERP_2(fs,X) (_FP_FRAC_HIGH_##fs(X) &= ~_FP_OVERFLOW_##fs)
138 #define _FP_FRAC_EQ_2(X, Y) (X##_f1 == Y##_f1 && X##_f0 == Y##_f0)
139 #define _FP_FRAC_GT_2(X, Y) \
140 (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 > Y##_f0))
141 #define _FP_FRAC_GE_2(X, Y) \
142 (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 >= Y##_f0))
143
144 #define _FP_ZEROFRAC_2 0, 0
145 #define _FP_MINFRAC_2 0, 1
146 #define _FP_MAXFRAC_2 (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0)
147
148 /*
149 * Internals
150 */
151
152 #define __FP_FRAC_SET_2(X,I1,I0) (X##_f0 = I0, X##_f1 = I1)
153
154 #define __FP_CLZ_2(R, xh, xl) \
155 do { \
156 if (xh) \
157 __FP_CLZ(R,xh); \
158 else \
159 { \
160 __FP_CLZ(R,xl); \
161 R += _FP_W_TYPE_SIZE; \
162 } \
163 } while(0)
164
165 #if 0
166
167 #ifndef __FP_FRAC_ADDI_2
168 #define __FP_FRAC_ADDI_2(xh, xl, i) \
169 (xh += ((xl += i) < i))
170 #endif
171 #ifndef __FP_FRAC_ADD_2
172 #define __FP_FRAC_ADD_2(rh, rl, xh, xl, yh, yl) \
173 (rh = xh + yh + ((rl = xl + yl) < xl))
174 #endif
175 #ifndef __FP_FRAC_SUB_2
176 #define __FP_FRAC_SUB_2(rh, rl, xh, xl, yh, yl) \
177 (rh = xh - yh - ((rl = xl - yl) > xl))
178 #endif
179 #ifndef __FP_FRAC_DEC_2
180 #define __FP_FRAC_DEC_2(xh, xl, yh, yl) \
181 do { \
182 UWtype _t = xl; \
183 xh -= yh + ((xl -= yl) > _t); \
184 } while (0)
185 #endif
186
187 #else
188
189 #undef __FP_FRAC_ADDI_2
190 #define __FP_FRAC_ADDI_2(xh, xl, i) add_ssaaaa(xh, xl, xh, xl, 0, i)
191 #undef __FP_FRAC_ADD_2
192 #define __FP_FRAC_ADD_2 add_ssaaaa
193 #undef __FP_FRAC_SUB_2
194 #define __FP_FRAC_SUB_2 sub_ddmmss
195 #undef __FP_FRAC_DEC_2
196 #define __FP_FRAC_DEC_2(xh, xl, yh, yl) sub_ddmmss(xh, xl, xh, xl, yh, yl)
197
198 #endif
199
200 /*
201 * Unpack the raw bits of a native fp value. Do not classify or
202 * normalize the data.
203 */
204
205 #define _FP_UNPACK_RAW_2(fs, X, val) \
206 do { \
207 union _FP_UNION_##fs _flo; _flo.flt = (val); \
208 \
209 X##_f0 = _flo.bits.frac0; \
210 X##_f1 = _flo.bits.frac1; \
211 X##_e = _flo.bits.exp; \
212 X##_s = _flo.bits.sign; \
213 } while (0)
214
215 #define _FP_UNPACK_RAW_2_P(fs, X, val) \
216 do { \
217 union _FP_UNION_##fs *_flo = \
218 (union _FP_UNION_##fs *)(val); \
219 \
220 X##_f0 = _flo->bits.frac0; \
221 X##_f1 = _flo->bits.frac1; \
222 X##_e = _flo->bits.exp; \
223 X##_s = _flo->bits.sign; \
224 } while (0)
225
226
227 /*
228 * Repack the raw bits of a native fp value.
229 */
230
231 #define _FP_PACK_RAW_2(fs, val, X) \
232 do { \
233 union _FP_UNION_##fs _flo; \
234 \
235 _flo.bits.frac0 = X##_f0; \
236 _flo.bits.frac1 = X##_f1; \
237 _flo.bits.exp = X##_e; \
238 _flo.bits.sign = X##_s; \
239 \
240 (val) = _flo.flt; \
241 } while (0)
242
243 #define _FP_PACK_RAW_2_P(fs, val, X) \
244 do { \
245 union _FP_UNION_##fs *_flo = \
246 (union _FP_UNION_##fs *)(val); \
247 \
248 _flo->bits.frac0 = X##_f0; \
249 _flo->bits.frac1 = X##_f1; \
250 _flo->bits.exp = X##_e; \
251 _flo->bits.sign = X##_s; \
252 } while (0)
253
254
255 /*
256 * Multiplication algorithms:
257 */
258
259 /* Given a 1W * 1W => 2W primitive, do the extended multiplication. */
260
261 #define _FP_MUL_MEAT_2_wide(wfracbits, R, X, Y, doit) \
262 do { \
263 _FP_FRAC_DECL_4(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c); \
264 \
265 doit(_FP_FRAC_WORD_4(_z,1), _FP_FRAC_WORD_4(_z,0), X##_f0, Y##_f0); \
266 doit(_b_f1, _b_f0, X##_f0, Y##_f1); \
267 doit(_c_f1, _c_f0, X##_f1, Y##_f0); \
268 doit(_FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2), X##_f1, Y##_f1); \
269 \
270 __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
271 _FP_FRAC_WORD_4(_z,1), 0, _b_f1, _b_f0, \
272 _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
273 _FP_FRAC_WORD_4(_z,1)); \
274 __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
275 _FP_FRAC_WORD_4(_z,1), 0, _c_f1, _c_f0, \
276 _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
277 _FP_FRAC_WORD_4(_z,1)); \
278 \
279 /* Normalize since we know where the msb of the multiplicands \
280 were (bit B), we know that the msb of the of the product is \
281 at either 2B or 2B-1. */ \
282 _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits); \
283 R##_f0 = _FP_FRAC_WORD_4(_z,0); \
284 R##_f1 = _FP_FRAC_WORD_4(_z,1); \
285 } while (0)
286
287 /* Given a 1W * 1W => 2W primitive, do the extended multiplication.
288 Do only 3 multiplications instead of four. This one is for machines
289 where multiplication is much more expensive than subtraction. */
290
291 #define _FP_MUL_MEAT_2_wide_3mul(wfracbits, R, X, Y, doit) \
292 do { \
293 _FP_FRAC_DECL_4(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c); \
294 _FP_W_TYPE _d; \
295 int _c1, _c2; \
296 \
297 _b_f0 = X##_f0 + X##_f1; \
298 _c1 = _b_f0 < X##_f0; \
299 _b_f1 = Y##_f0 + Y##_f1; \
300 _c2 = _b_f1 < Y##_f0; \
301 doit(_d, _FP_FRAC_WORD_4(_z,0), X##_f0, Y##_f0); \
302 doit(_FP_FRAC_WORD_4(_z,2), _FP_FRAC_WORD_4(_z,1), _b_f0, _b_f1); \
303 doit(_c_f1, _c_f0, X##_f1, Y##_f1); \
304 \
305 _b_f0 &= -_c2; \
306 _b_f1 &= -_c1; \
307 __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
308 _FP_FRAC_WORD_4(_z,1), (_c1 & _c2), 0, _d, \
309 0, _FP_FRAC_WORD_4(_z,2), _FP_FRAC_WORD_4(_z,1)); \
310 __FP_FRAC_ADDI_2(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
311 _b_f0); \
312 __FP_FRAC_ADDI_2(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
313 _b_f1); \
314 __FP_FRAC_DEC_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
315 _FP_FRAC_WORD_4(_z,1), \
316 0, _d, _FP_FRAC_WORD_4(_z,0)); \
317 __FP_FRAC_DEC_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2), \
318 _FP_FRAC_WORD_4(_z,1), 0, _c_f1, _c_f0); \
319 __FP_FRAC_ADD_2(_FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2), \
320 _c_f1, _c_f0, \
321 _FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2)); \
322 \
323 /* Normalize since we know where the msb of the multiplicands \
324 were (bit B), we know that the msb of the of the product is \
325 at either 2B or 2B-1. */ \
326 _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits); \
327 R##_f0 = _FP_FRAC_WORD_4(_z,0); \
328 R##_f1 = _FP_FRAC_WORD_4(_z,1); \
329 } while (0)
330
331 #define _FP_MUL_MEAT_2_gmp(wfracbits, R, X, Y) \
332 do { \
333 _FP_FRAC_DECL_4(_z); \
334 _FP_W_TYPE _x[2], _y[2]; \
335 _x[0] = X##_f0; _x[1] = X##_f1; \
336 _y[0] = Y##_f0; _y[1] = Y##_f1; \
337 \
338 mpn_mul_n(_z_f, _x, _y, 2); \
339 \
340 /* Normalize since we know where the msb of the multiplicands \
341 were (bit B), we know that the msb of the of the product is \
342 at either 2B or 2B-1. */ \
343 _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits); \
344 R##_f0 = _z_f[0]; \
345 R##_f1 = _z_f[1]; \
346 } while (0)
347
348 /* Do at most 120x120=240 bits multiplication using double floating
349 point multiplication. This is useful if floating point
350 multiplication has much bigger throughput than integer multiply.
351 It is supposed to work for _FP_W_TYPE_SIZE 64 and wfracbits
352 between 106 and 120 only.
353 Caller guarantees that X and Y has (1LLL << (wfracbits - 1)) set.
354 SETFETZ is a macro which will disable all FPU exceptions and set rounding
355 towards zero, RESETFE should optionally reset it back. */
356
357 #define _FP_MUL_MEAT_2_120_240_double(wfracbits, R, X, Y, setfetz, resetfe) \
358 do { \
359 static const double _const[] = { \
360 /* 2^-24 */ 5.9604644775390625e-08, \
361 /* 2^-48 */ 3.5527136788005009e-15, \
362 /* 2^-72 */ 2.1175823681357508e-22, \
363 /* 2^-96 */ 1.2621774483536189e-29, \
364 /* 2^28 */ 2.68435456e+08, \
365 /* 2^4 */ 1.600000e+01, \
366 /* 2^-20 */ 9.5367431640625e-07, \
367 /* 2^-44 */ 5.6843418860808015e-14, \
368 /* 2^-68 */ 3.3881317890172014e-21, \
369 /* 2^-92 */ 2.0194839173657902e-28, \
370 /* 2^-116 */ 1.2037062152420224e-35}; \
371 double _a240, _b240, _c240, _d240, _e240, _f240, \
372 _g240, _h240, _i240, _j240, _k240; \
373 union { double d; UDItype i; } _l240, _m240, _n240, _o240, \
374 _p240, _q240, _r240, _s240; \
375 UDItype _t240, _u240, _v240, _w240, _x240, _y240 = 0; \
376 \
377 if (wfracbits < 106 || wfracbits > 120) \
378 abort(); \
379 \
380 setfetz; \
381 \
382 _e240 = (double)(long)(X##_f0 & 0xffffff); \
383 _j240 = (double)(long)(Y##_f0 & 0xffffff); \
384 _d240 = (double)(long)((X##_f0 >> 24) & 0xffffff); \
385 _i240 = (double)(long)((Y##_f0 >> 24) & 0xffffff); \
386 _c240 = (double)(long)(((X##_f1 << 16) & 0xffffff) | (X##_f0 >> 48)); \
387 _h240 = (double)(long)(((Y##_f1 << 16) & 0xffffff) | (Y##_f0 >> 48)); \
388 _b240 = (double)(long)((X##_f1 >> 8) & 0xffffff); \
389 _g240 = (double)(long)((Y##_f1 >> 8) & 0xffffff); \
390 _a240 = (double)(long)(X##_f1 >> 32); \
391 _f240 = (double)(long)(Y##_f1 >> 32); \
392 _e240 *= _const[3]; \
393 _j240 *= _const[3]; \
394 _d240 *= _const[2]; \
395 _i240 *= _const[2]; \
396 _c240 *= _const[1]; \
397 _h240 *= _const[1]; \
398 _b240 *= _const[0]; \
399 _g240 *= _const[0]; \
400 _s240.d = _e240*_j240;\
401 _r240.d = _d240*_j240 + _e240*_i240;\
402 _q240.d = _c240*_j240 + _d240*_i240 + _e240*_h240;\
403 _p240.d = _b240*_j240 + _c240*_i240 + _d240*_h240 + _e240*_g240;\
404 _o240.d = _a240*_j240 + _b240*_i240 + _c240*_h240 + _d240*_g240 + _e240*_f240;\
405 _n240.d = _a240*_i240 + _b240*_h240 + _c240*_g240 + _d240*_f240; \
406 _m240.d = _a240*_h240 + _b240*_g240 + _c240*_f240; \
407 _l240.d = _a240*_g240 + _b240*_f240; \
408 _k240 = _a240*_f240; \
409 _r240.d += _s240.d; \
410 _q240.d += _r240.d; \
411 _p240.d += _q240.d; \
412 _o240.d += _p240.d; \
413 _n240.d += _o240.d; \
414 _m240.d += _n240.d; \
415 _l240.d += _m240.d; \
416 _k240 += _l240.d; \
417 _s240.d -= ((_const[10]+_s240.d)-_const[10]); \
418 _r240.d -= ((_const[9]+_r240.d)-_const[9]); \
419 _q240.d -= ((_const[8]+_q240.d)-_const[8]); \
420 _p240.d -= ((_const[7]+_p240.d)-_const[7]); \
421 _o240.d += _const[7]; \
422 _n240.d += _const[6]; \
423 _m240.d += _const[5]; \
424 _l240.d += _const[4]; \
425 if (_s240.d != 0.0) _y240 = 1; \
426 if (_r240.d != 0.0) _y240 = 1; \
427 if (_q240.d != 0.0) _y240 = 1; \
428 if (_p240.d != 0.0) _y240 = 1; \
429 _t240 = (DItype)_k240; \
430 _u240 = _l240.i; \
431 _v240 = _m240.i; \
432 _w240 = _n240.i; \
433 _x240 = _o240.i; \
434 R##_f1 = (_t240 << (128 - (wfracbits - 1))) \
435 | ((_u240 & 0xffffff) >> ((wfracbits - 1) - 104)); \
436 R##_f0 = ((_u240 & 0xffffff) << (168 - (wfracbits - 1))) \
437 | ((_v240 & 0xffffff) << (144 - (wfracbits - 1))) \
438 | ((_w240 & 0xffffff) << (120 - (wfracbits - 1))) \
439 | ((_x240 & 0xffffff) >> ((wfracbits - 1) - 96)) \
440 | _y240; \
441 resetfe; \
442 } while (0)
443
444 /*
445 * Division algorithms:
446 */
447
448 #define _FP_DIV_MEAT_2_udiv(fs, R, X, Y) \
449 do { \
450 _FP_W_TYPE _n_f2, _n_f1, _n_f0, _r_f1, _r_f0, _m_f1, _m_f0; \
451 if (_FP_FRAC_GT_2(X, Y)) \
452 { \
453 _n_f2 = X##_f1 >> 1; \
454 _n_f1 = X##_f1 << (_FP_W_TYPE_SIZE - 1) | X##_f0 >> 1; \
455 _n_f0 = X##_f0 << (_FP_W_TYPE_SIZE - 1); \
456 } \
457 else \
458 { \
459 R##_e--; \
460 _n_f2 = X##_f1; \
461 _n_f1 = X##_f0; \
462 _n_f0 = 0; \
463 } \
464 \
465 /* Normalize, i.e. make the most significant bit of the \
466 denominator set. */ \
467 _FP_FRAC_SLL_2(Y, _FP_WFRACXBITS_##fs); \
468 \
469 udiv_qrnnd(R##_f1, _r_f1, _n_f2, _n_f1, Y##_f1); \
470 umul_ppmm(_m_f1, _m_f0, R##_f1, Y##_f0); \
471 _r_f0 = _n_f0; \
472 if (_FP_FRAC_GT_2(_m, _r)) \
473 { \
474 R##_f1--; \
475 _FP_FRAC_ADD_2(_r, Y, _r); \
476 if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r)) \
477 { \
478 R##_f1--; \
479 _FP_FRAC_ADD_2(_r, Y, _r); \
480 } \
481 } \
482 _FP_FRAC_DEC_2(_r, _m); \
483 \
484 if (_r_f1 == Y##_f1) \
485 { \
486 /* This is a special case, not an optimization \
487 (_r/Y##_f1 would not fit into UWtype). \
488 As _r is guaranteed to be < Y, R##_f0 can be either \
489 (UWtype)-1 or (UWtype)-2. But as we know what kind \
490 of bits it is (sticky, guard, round), we don't care. \
491 We also don't care what the reminder is, because the \
492 guard bit will be set anyway. -jj */ \
493 R##_f0 = -1; \
494 } \
495 else \
496 { \
497 udiv_qrnnd(R##_f0, _r_f1, _r_f1, _r_f0, Y##_f1); \
498 umul_ppmm(_m_f1, _m_f0, R##_f0, Y##_f0); \
499 _r_f0 = 0; \
500 if (_FP_FRAC_GT_2(_m, _r)) \
501 { \
502 R##_f0--; \
503 _FP_FRAC_ADD_2(_r, Y, _r); \
504 if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r)) \
505 { \
506 R##_f0--; \
507 _FP_FRAC_ADD_2(_r, Y, _r); \
508 } \
509 } \
510 if (!_FP_FRAC_EQ_2(_r, _m)) \
511 R##_f0 |= _FP_WORK_STICKY; \
512 } \
513 } while (0)
514
515
516 #define _FP_DIV_MEAT_2_gmp(fs, R, X, Y) \
517 do { \
518 _FP_W_TYPE _x[4], _y[2], _z[4]; \
519 _y[0] = Y##_f0; _y[1] = Y##_f1; \
520 _x[0] = _x[3] = 0; \
521 if (_FP_FRAC_GT_2(X, Y)) \
522 { \
523 R##_e++; \
524 _x[1] = (X##_f0 << (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE) | \
525 X##_f1 >> (_FP_W_TYPE_SIZE - \
526 (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE))); \
527 _x[2] = X##_f1 << (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE); \
528 } \
529 else \
530 { \
531 _x[1] = (X##_f0 << (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE) | \
532 X##_f1 >> (_FP_W_TYPE_SIZE - \
533 (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE))); \
534 _x[2] = X##_f1 << (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE); \
535 } \
536 \
537 (void) mpn_divrem (_z, 0, _x, 4, _y, 2); \
538 R##_f1 = _z[1]; \
539 R##_f0 = _z[0] | ((_x[0] | _x[1]) != 0); \
540 } while (0)
541
542
543 /*
544 * Square root algorithms:
545 * We have just one right now, maybe Newton approximation
546 * should be added for those machines where division is fast.
547 */
548
549 #define _FP_SQRT_MEAT_2(R, S, T, X, q) \
550 do { \
551 while (q) \
552 { \
553 T##_f1 = S##_f1 + q; \
554 if (T##_f1 <= X##_f1) \
555 { \
556 S##_f1 = T##_f1 + q; \
557 X##_f1 -= T##_f1; \
558 R##_f1 += q; \
559 } \
560 _FP_FRAC_SLL_2(X, 1); \
561 q >>= 1; \
562 } \
563 q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1); \
564 while (q != _FP_WORK_ROUND) \
565 { \
566 T##_f0 = S##_f0 + q; \
567 T##_f1 = S##_f1; \
568 if (T##_f1 < X##_f1 || \
569 (T##_f1 == X##_f1 && T##_f0 <= X##_f0)) \
570 { \
571 S##_f0 = T##_f0 + q; \
572 S##_f1 += (T##_f0 > S##_f0); \
573 _FP_FRAC_DEC_2(X, T); \
574 R##_f0 += q; \
575 } \
576 _FP_FRAC_SLL_2(X, 1); \
577 q >>= 1; \
578 } \
579 if (X##_f0 | X##_f1) \
580 { \
581 if (S##_f1 < X##_f1 || \
582 (S##_f1 == X##_f1 && S##_f0 < X##_f0)) \
583 R##_f0 |= _FP_WORK_ROUND; \
584 R##_f0 |= _FP_WORK_STICKY; \
585 } \
586 } while (0)
587
588
589 /*
590 * Assembly/disassembly for converting to/from integral types.
591 * No shifting or overflow handled here.
592 */
593
594 #define _FP_FRAC_ASSEMBLE_2(r, X, rsize) \
595 (void)((rsize <= _FP_W_TYPE_SIZE) \
596 ? ({ r = X##_f0; }) \
597 : ({ \
598 r = X##_f1; \
599 r <<= _FP_W_TYPE_SIZE; \
600 r += X##_f0; \
601 }))
602
603 #define _FP_FRAC_DISASSEMBLE_2(X, r, rsize) \
604 do { \
605 X##_f0 = r; \
606 X##_f1 = (rsize <= _FP_W_TYPE_SIZE ? 0 : r >> _FP_W_TYPE_SIZE); \
607 } while (0)
608
609 /*
610 * Convert FP values between word sizes
611 */
612
613 #define _FP_FRAC_COPY_1_2(D, S) (D##_f = S##_f0)
614
615 #define _FP_FRAC_COPY_2_1(D, S) ((D##_f0 = S##_f), (D##_f1 = 0))
Unaltered CVS history