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