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