84770fceb53e053cc8e2647c21710729c11bf051
[safe/jmp/linux-2.6] / include / math-emu / extended.h
1 /* Software floating-point emulation.
2    Definitions for IEEE Extended Precision.
3    Copyright (C) 1999 Free Software Foundation, Inc.
4    This file is part of the GNU C Library.
5    Contributed by Jakub Jelinek (jj@ultra.linux.cz).
6
7    The GNU C Library is free software; you can redistribute it and/or
8    modify it under the terms of the GNU Library General Public License as
9    published by the Free Software Foundation; either version 2 of the
10    License, or (at your option) any later version.
11
12    The GNU C Library is distributed in the hope that it will be useful,
13    but WITHOUT ANY WARRANTY; without even the implied warranty of
14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15    Library General Public License for more details.
16
17    You should have received a copy of the GNU Library General Public
18    License along with the GNU C Library; see the file COPYING.LIB.  If
19    not, write to the Free Software Foundation, Inc.,
20    59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
21
22
23 #ifndef    __MATH_EMU_EXTENDED_H__
24 #define    __MATH_EMU_EXTENDED_H__
25
26 #if _FP_W_TYPE_SIZE < 32
27 #error "Here's a nickel, kid. Go buy yourself a real computer."
28 #endif
29
30 #if _FP_W_TYPE_SIZE < 64
31 #define _FP_FRACTBITS_E         (4*_FP_W_TYPE_SIZE)
32 #else
33 #define _FP_FRACTBITS_E         (2*_FP_W_TYPE_SIZE)
34 #endif
35
36 #define _FP_FRACBITS_E          64
37 #define _FP_FRACXBITS_E         (_FP_FRACTBITS_E - _FP_FRACBITS_E)
38 #define _FP_WFRACBITS_E         (_FP_WORKBITS + _FP_FRACBITS_E)
39 #define _FP_WFRACXBITS_E        (_FP_FRACTBITS_E - _FP_WFRACBITS_E)
40 #define _FP_EXPBITS_E           15
41 #define _FP_EXPBIAS_E           16383
42 #define _FP_EXPMAX_E            32767
43
44 #define _FP_QNANBIT_E           \
45         ((_FP_W_TYPE)1 << (_FP_FRACBITS_E-2) % _FP_W_TYPE_SIZE)
46 #define _FP_IMPLBIT_E           \
47         ((_FP_W_TYPE)1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE)
48 #define _FP_OVERFLOW_E          \
49         ((_FP_W_TYPE)1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE))
50
51 #if _FP_W_TYPE_SIZE < 64
52
53 union _FP_UNION_E
54 {
55    long double flt;
56    struct 
57    {
58 #if __BYTE_ORDER == __BIG_ENDIAN
59       unsigned long pad1 : _FP_W_TYPE_SIZE;
60       unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
61       unsigned long sign : 1;
62       unsigned long exp : _FP_EXPBITS_E;
63       unsigned long frac1 : _FP_W_TYPE_SIZE;
64       unsigned long frac0 : _FP_W_TYPE_SIZE;
65 #else
66       unsigned long frac0 : _FP_W_TYPE_SIZE;
67       unsigned long frac1 : _FP_W_TYPE_SIZE;
68       unsigned exp : _FP_EXPBITS_E;
69       unsigned sign : 1;
70 #endif /* not bigendian */
71    } bits __attribute__((packed));
72 };
73
74
75 #define FP_DECL_E(X)            _FP_DECL(4,X)
76
77 #define FP_UNPACK_RAW_E(X, val)                         \
78   do {                                                  \
79     union _FP_UNION_E _flo; _flo.flt = (val);           \
80                                                         \
81     X##_f[2] = 0; X##_f[3] = 0;                         \
82     X##_f[0] = _flo.bits.frac0;                         \
83     X##_f[1] = _flo.bits.frac1;                         \
84     X##_e  = _flo.bits.exp;                             \
85     X##_s  = _flo.bits.sign;                            \
86     if (!X##_e && (X##_f[1] || X##_f[0])                \
87         && !(X##_f[1] & _FP_IMPLBIT_E))                 \
88       {                                                 \
89         X##_e++;                                        \
90         FP_SET_EXCEPTION(FP_EX_DENORM);                 \
91       }                                                 \
92   } while (0)
93
94 #define FP_UNPACK_RAW_EP(X, val)                        \
95   do {                                                  \
96     union _FP_UNION_E *_flo =                           \
97     (union _FP_UNION_E *)(val);                         \
98                                                         \
99     X##_f[2] = 0; X##_f[3] = 0;                         \
100     X##_f[0] = _flo->bits.frac0;                        \
101     X##_f[1] = _flo->bits.frac1;                        \
102     X##_e  = _flo->bits.exp;                            \
103     X##_s  = _flo->bits.sign;                           \
104     if (!X##_e && (X##_f[1] || X##_f[0])                \
105         && !(X##_f[1] & _FP_IMPLBIT_E))                 \
106       {                                                 \
107         X##_e++;                                        \
108         FP_SET_EXCEPTION(FP_EX_DENORM);                 \
109       }                                                 \
110   } while (0)
111
112 #define FP_PACK_RAW_E(val, X)                           \
113   do {                                                  \
114     union _FP_UNION_E _flo;                             \
115                                                         \
116     if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;               \
117     else X##_f[1] &= ~(_FP_IMPLBIT_E);                  \
118     _flo.bits.frac0 = X##_f[0];                         \
119     _flo.bits.frac1 = X##_f[1];                         \
120     _flo.bits.exp   = X##_e;                            \
121     _flo.bits.sign  = X##_s;                            \
122                                                         \
123     (val) = _flo.flt;                                   \
124   } while (0)
125
126 #define FP_PACK_RAW_EP(val, X)                          \
127   do {                                                  \
128     if (!FP_INHIBIT_RESULTS)                            \
129       {                                                 \
130         union _FP_UNION_E *_flo =                       \
131           (union _FP_UNION_E *)(val);                   \
132                                                         \
133         if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;           \
134         else X##_f[1] &= ~(_FP_IMPLBIT_E);              \
135         _flo->bits.frac0 = X##_f[0];                    \
136         _flo->bits.frac1 = X##_f[1];                    \
137         _flo->bits.exp   = X##_e;                       \
138         _flo->bits.sign  = X##_s;                       \
139       }                                                 \
140   } while (0)
141
142 #define FP_UNPACK_E(X,val)              \
143   do {                                  \
144     FP_UNPACK_RAW_E(X,val);             \
145     _FP_UNPACK_CANONICAL(E,4,X);        \
146   } while (0)
147
148 #define FP_UNPACK_EP(X,val)             \
149   do {                                  \
150     FP_UNPACK_RAW_2_P(X,val);           \
151     _FP_UNPACK_CANONICAL(E,4,X);        \
152   } while (0)
153
154 #define FP_PACK_E(val,X)                \
155   do {                                  \
156     _FP_PACK_CANONICAL(E,4,X);          \
157     FP_PACK_RAW_E(val,X);               \
158   } while (0)
159
160 #define FP_PACK_EP(val,X)               \
161   do {                                  \
162     _FP_PACK_CANONICAL(E,4,X);          \
163     FP_PACK_RAW_EP(val,X);              \
164   } while (0)
165
166 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,4,X)
167 #define FP_NEG_E(R,X)           _FP_NEG(E,4,R,X)
168 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,4,R,X,Y)
169 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,4,R,X,Y)
170 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,4,R,X,Y)
171 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,4,R,X,Y)
172 #define FP_SQRT_E(R,X)          _FP_SQRT(E,4,R,X)
173
174 /*
175  * Square root algorithms:
176  * We have just one right now, maybe Newton approximation
177  * should be added for those machines where division is fast.
178  * This has special _E version because standard _4 square
179  * root would not work (it has to start normally with the
180  * second word and not the first), but as we have to do it
181  * anyway, we optimize it by doing most of the calculations
182  * in two UWtype registers instead of four.
183  */
184  
185 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
186   do {                                                  \
187     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
188     _FP_FRAC_SRL_4(X, (_FP_WORKBITS));                  \
189     while (q)                                           \
190       {                                                 \
191         T##_f[1] = S##_f[1] + q;                        \
192         if (T##_f[1] <= X##_f[1])                       \
193           {                                             \
194             S##_f[1] = T##_f[1] + q;                    \
195             X##_f[1] -= T##_f[1];                       \
196             R##_f[1] += q;                              \
197           }                                             \
198         _FP_FRAC_SLL_2(X, 1);                           \
199         q >>= 1;                                        \
200       }                                                 \
201     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
202     while (q)                                           \
203       {                                                 \
204         T##_f[0] = S##_f[0] + q;                        \
205         T##_f[1] = S##_f[1];                            \
206         if (T##_f[1] < X##_f[1] ||                      \
207             (T##_f[1] == X##_f[1] &&                    \
208              T##_f[0] <= X##_f[0]))                     \
209           {                                             \
210             S##_f[0] = T##_f[0] + q;                    \
211             S##_f[1] += (T##_f[0] > S##_f[0]);          \
212             _FP_FRAC_DEC_2(X, T);                       \
213             R##_f[0] += q;                              \
214           }                                             \
215         _FP_FRAC_SLL_2(X, 1);                           \
216         q >>= 1;                                        \
217       }                                                 \
218     _FP_FRAC_SLL_4(R, (_FP_WORKBITS));                  \
219     if (X##_f[0] | X##_f[1])                            \
220       {                                                 \
221         if (S##_f[1] < X##_f[1] ||                      \
222             (S##_f[1] == X##_f[1] &&                    \
223              S##_f[0] < X##_f[0]))                      \
224           R##_f[0] |= _FP_WORK_ROUND;                   \
225         R##_f[0] |= _FP_WORK_STICKY;                    \
226       }                                                 \
227   } while (0)
228
229 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,4,r,X,Y,un)
230 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,4,r,X,Y)
231
232 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,4,r,X,rsz,rsg)
233 #define FP_TO_INT_ROUND_E(r,X,rsz,rsg)  _FP_TO_INT_ROUND(E,4,r,X,rsz,rsg)
234 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,4,X,r,rs,rt)
235
236 #define _FP_FRAC_HIGH_E(X)      (X##_f[2])
237 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f[1])
238
239 #else   /* not _FP_W_TYPE_SIZE < 64 */
240 union _FP_UNION_E
241 {
242   long double flt /* __attribute__((mode(TF))) */ ;
243   struct {
244 #if __BYTE_ORDER == __BIG_ENDIAN
245     unsigned long pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
246     unsigned sign  : 1;
247     unsigned exp   : _FP_EXPBITS_E;
248     unsigned long frac : _FP_W_TYPE_SIZE;
249 #else
250     unsigned long frac : _FP_W_TYPE_SIZE;
251     unsigned exp   : _FP_EXPBITS_E;
252     unsigned sign  : 1;
253 #endif
254   } bits;
255 };
256
257 #define FP_DECL_E(X)            _FP_DECL(2,X)
258
259 #define FP_UNPACK_RAW_E(X, val)                                 \
260   do {                                                          \
261     union _FP_UNION_E _flo; _flo.flt = (val);                   \
262                                                                 \
263     X##_f0 = _flo.bits.frac;                                    \
264     X##_f1 = 0;                                                 \
265     X##_e = _flo.bits.exp;                                      \
266     X##_s = _flo.bits.sign;                                     \
267     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
268       {                                                         \
269         X##_e++;                                                \
270         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
271       }                                                         \
272   } while (0)
273
274 #define FP_UNPACK_RAW_EP(X, val)                                \
275   do {                                                          \
276     union _FP_UNION_E *_flo =                                   \
277       (union _FP_UNION_E *)(val);                               \
278                                                                 \
279     X##_f0 = _flo->bits.frac;                                   \
280     X##_f1 = 0;                                                 \
281     X##_e = _flo->bits.exp;                                     \
282     X##_s = _flo->bits.sign;                                    \
283     if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))          \
284       {                                                         \
285         X##_e++;                                                \
286         FP_SET_EXCEPTION(FP_EX_DENORM);                         \
287       }                                                         \
288   } while (0)
289
290 #define FP_PACK_RAW_E(val, X)                                   \
291   do {                                                          \
292     union _FP_UNION_E _flo;                                     \
293                                                                 \
294     if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                         \
295     else X##_f0 &= ~(_FP_IMPLBIT_E);                            \
296     _flo.bits.frac = X##_f0;                                    \
297     _flo.bits.exp  = X##_e;                                     \
298     _flo.bits.sign = X##_s;                                     \
299                                                                 \
300     (val) = _flo.flt;                                           \
301   } while (0)
302
303 #define FP_PACK_RAW_EP(fs, val, X)                              \
304   do {                                                          \
305     if (!FP_INHIBIT_RESULTS)                                    \
306       {                                                         \
307         union _FP_UNION_E *_flo =                               \
308           (union _FP_UNION_E *)(val);                           \
309                                                                 \
310         if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                     \
311         else X##_f0 &= ~(_FP_IMPLBIT_E);                        \
312         _flo->bits.frac = X##_f0;                               \
313         _flo->bits.exp  = X##_e;                                \
314         _flo->bits.sign = X##_s;                                \
315       }                                                         \
316   } while (0)
317
318
319 #define FP_UNPACK_E(X,val)              \
320   do {                                  \
321     FP_UNPACK_RAW_E(X,val);             \
322     _FP_UNPACK_CANONICAL(E,2,X);        \
323   } while (0)
324
325 #define FP_UNPACK_EP(X,val)             \
326   do {                                  \
327     FP_UNPACK_RAW_EP(X,val);            \
328     _FP_UNPACK_CANONICAL(E,2,X);        \
329   } while (0)
330
331 #define FP_PACK_E(val,X)                \
332   do {                                  \
333     _FP_PACK_CANONICAL(E,2,X);          \
334     FP_PACK_RAW_E(val,X);               \
335   } while (0)
336
337 #define FP_PACK_EP(val,X)               \
338   do {                                  \
339     _FP_PACK_CANONICAL(E,2,X);          \
340     FP_PACK_RAW_EP(val,X);              \
341   } while (0)
342
343 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,2,X)
344 #define FP_NEG_E(R,X)           _FP_NEG(E,2,R,X)
345 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,2,R,X,Y)
346 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,2,R,X,Y)
347 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,2,R,X,Y)
348 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,2,R,X,Y)
349 #define FP_SQRT_E(R,X)          _FP_SQRT(E,2,R,X)
350
351 /*
352  * Square root algorithms:
353  * We have just one right now, maybe Newton approximation
354  * should be added for those machines where division is fast.
355  * We optimize it by doing most of the calculations
356  * in one UWtype registers instead of two, although we don't
357  * have to.
358  */
359 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
360   do {                                                  \
361     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
362     _FP_FRAC_SRL_2(X, (_FP_WORKBITS));                  \
363     while (q)                                           \
364       {                                                 \
365         T##_f0 = S##_f0 + q;                            \
366         if (T##_f0 <= X##_f0)                           \
367           {                                             \
368             S##_f0 = T##_f0 + q;                        \
369             X##_f0 -= T##_f0;                           \
370             R##_f0 += q;                                \
371           }                                             \
372         _FP_FRAC_SLL_1(X, 1);                           \
373         q >>= 1;                                        \
374       }                                                 \
375     _FP_FRAC_SLL_2(R, (_FP_WORKBITS));                  \
376     if (X##_f0)                                         \
377       {                                                 \
378         if (S##_f0 < X##_f0)                            \
379           R##_f0 |= _FP_WORK_ROUND;                     \
380         R##_f0 |= _FP_WORK_STICKY;                      \
381       }                                                 \
382   } while (0)
383  
384 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,2,r,X,Y,un)
385 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,2,r,X,Y)
386
387 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,2,r,X,rsz,rsg)
388 #define FP_TO_INT_ROUND_E(r,X,rsz,rsg)  _FP_TO_INT_ROUND(E,2,r,X,rsz,rsg)
389 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,2,X,r,rs,rt)
390
391 #define _FP_FRAC_HIGH_E(X)      (X##_f1)
392 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f0)
393
394 #endif /* not _FP_W_TYPE_SIZE < 64 */
395
396 #endif /* __MATH_EMU_EXTENDED_H__ */