~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~

TOMOYO Linux Cross Reference
Linux/include/math-emu/op-2.h

Version: ~ [ linux-5.1-rc2 ] ~ [ linux-5.0.3 ] ~ [ linux-4.20.17 ] ~ [ linux-4.19.30 ] ~ [ linux-4.18.20 ] ~ [ linux-4.17.19 ] ~ [ linux-4.16.18 ] ~ [ linux-4.15.18 ] ~ [ linux-4.14.107 ] ~ [ linux-4.13.16 ] ~ [ linux-4.12.14 ] ~ [ linux-4.11.12 ] ~ [ linux-4.10.17 ] ~ [ linux-4.9.164 ] ~ [ linux-4.8.17 ] ~ [ linux-4.7.10 ] ~ [ linux-4.6.7 ] ~ [ linux-4.5.7 ] ~ [ linux-4.4.176 ] ~ [ linux-4.3.6 ] ~ [ linux-4.2.8 ] ~ [ linux-4.1.52 ] ~ [ linux-4.0.9 ] ~ [ linux-3.19.8 ] ~ [ linux-3.18.136 ] ~ [ linux-3.17.8 ] ~ [ linux-3.16.63 ] ~ [ linux-3.15.10 ] ~ [ linux-3.14.79 ] ~ [ linux-3.13.11 ] ~ [ linux-3.12.74 ] ~ [ linux-3.11.10 ] ~ [ linux-3.10.108 ] ~ [ linux-3.9.11 ] ~ [ linux-3.8.13 ] ~ [ linux-3.7.10 ] ~ [ linux-3.6.11 ] ~ [ linux-3.5.7 ] ~ [ linux-3.4.113 ] ~ [ linux-3.3.8 ] ~ [ linux-3.2.102 ] ~ [ linux-3.1.10 ] ~ [ linux-3.0.101 ] ~ [ linux-2.6.39.4 ] ~ [ linux-2.6.38.8 ] ~ [ linux-2.6.37.6 ] ~ [ linux-2.6.36.4 ] ~ [ linux-2.6.35.14 ] ~ [ linux-2.6.34.15 ] ~ [ linux-2.6.33.20 ] ~ [ linux-2.6.32.71 ] ~ [ linux-2.6.0 ] ~ [ linux-2.4.37.11 ] ~ [ unix-v6-master ] ~ [ ccs-tools-1.8.5 ] ~ [ policy-sample ] ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

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

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~

kernel.org | git.kernel.org | LWN.net | Project Home | Wiki (Japanese) | Wiki (English) | SVN repository | Mail admin

Linux® is a registered trademark of Linus Torvalds in the United States and other countries.
TOMOYO® is a registered trademark of NTT DATA CORPORATION.

osdn.jp