1/* Software floating-point emulation.
2   Definitions for IEEE Extended Precision.
3   Copyright (C) 1999-2014 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   In addition to the permissions in the GNU Lesser General Public
13   License, the Free Software Foundation gives you unlimited
14   permission to link the compiled version of this file into
15   combinations with other programs, and to distribute those
16   combinations without any restriction coming from the use of this
17   file.  (The Lesser General Public License restrictions do apply in
18   other respects; for example, they cover modification of the file,
19   and distribution when not linked into a combine executable.)
20
21   The GNU C Library is distributed in the hope that it will be useful,
22   but WITHOUT ANY WARRANTY; without even the implied warranty of
23   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
24   Lesser General Public License for more details.
25
26   You should have received a copy of the GNU Lesser General Public
27   License along with the GNU C Library; if not, see
28   <http://www.gnu.org/licenses/>.  */
29
30#if _FP_W_TYPE_SIZE < 32
31# error "Here's a nickel, kid. Go buy yourself a real computer."
32#endif
33
34#if _FP_W_TYPE_SIZE < 64
35# define _FP_FRACTBITS_E	(4*_FP_W_TYPE_SIZE)
36# define _FP_FRACTBITS_DW_E	(8*_FP_W_TYPE_SIZE)
37#else
38# define _FP_FRACTBITS_E	(2*_FP_W_TYPE_SIZE)
39# define _FP_FRACTBITS_DW_E	(4*_FP_W_TYPE_SIZE)
40#endif
41
42#define _FP_FRACBITS_E		64
43#define _FP_FRACXBITS_E		(_FP_FRACTBITS_E - _FP_FRACBITS_E)
44#define _FP_WFRACBITS_E		(_FP_WORKBITS + _FP_FRACBITS_E)
45#define _FP_WFRACXBITS_E	(_FP_FRACTBITS_E - _FP_WFRACBITS_E)
46#define _FP_EXPBITS_E		15
47#define _FP_EXPBIAS_E		16383
48#define _FP_EXPMAX_E		32767
49
50#define _FP_QNANBIT_E		\
51	((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-2) % _FP_W_TYPE_SIZE)
52#define _FP_QNANBIT_SH_E		\
53	((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-2+_FP_WORKBITS) % _FP_W_TYPE_SIZE)
54#define _FP_IMPLBIT_E		\
55	((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE)
56#define _FP_IMPLBIT_SH_E		\
57	((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-1+_FP_WORKBITS) % _FP_W_TYPE_SIZE)
58#define _FP_OVERFLOW_E		\
59	((_FP_W_TYPE) 1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE))
60
61#define _FP_WFRACBITS_DW_E	(2 * _FP_WFRACBITS_E)
62#define _FP_WFRACXBITS_DW_E	(_FP_FRACTBITS_DW_E - _FP_WFRACBITS_DW_E)
63#define _FP_HIGHBIT_DW_E	\
64  ((_FP_W_TYPE) 1 << (_FP_WFRACBITS_DW_E - 1) % _FP_W_TYPE_SIZE)
65
66typedef float XFtype __attribute__ ((mode (XF)));
67
68#if _FP_W_TYPE_SIZE < 64
69
70union _FP_UNION_E
71{
72  XFtype flt;
73  struct _FP_STRUCT_LAYOUT
74  {
75# if __BYTE_ORDER == __BIG_ENDIAN
76    unsigned long pad1 : _FP_W_TYPE_SIZE;
77    unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
78    unsigned long sign : 1;
79    unsigned long exp : _FP_EXPBITS_E;
80    unsigned long frac1 : _FP_W_TYPE_SIZE;
81    unsigned long frac0 : _FP_W_TYPE_SIZE;
82# else
83    unsigned long frac0 : _FP_W_TYPE_SIZE;
84    unsigned long frac1 : _FP_W_TYPE_SIZE;
85    unsigned exp : _FP_EXPBITS_E;
86    unsigned sign : 1;
87# endif /* not bigendian */
88  } bits __attribute__ ((packed));
89};
90
91
92# define FP_DECL_E(X)		_FP_DECL (4, X)
93
94# define FP_UNPACK_RAW_E(X, val)			\
95  do							\
96    {							\
97      union _FP_UNION_E FP_UNPACK_RAW_E_flo;		\
98      FP_UNPACK_RAW_E_flo.flt = (val);			\
99							\
100      X##_f[2] = 0;					\
101      X##_f[3] = 0;					\
102      X##_f[0] = FP_UNPACK_RAW_E_flo.bits.frac0;	\
103      X##_f[1] = FP_UNPACK_RAW_E_flo.bits.frac1;	\
104      X##_e  = FP_UNPACK_RAW_E_flo.bits.exp;		\
105      X##_s  = FP_UNPACK_RAW_E_flo.bits.sign;		\
106    }							\
107  while (0)
108
109# define FP_UNPACK_RAW_EP(X, val)			\
110  do							\
111    {							\
112      union _FP_UNION_E *FP_UNPACK_RAW_EP_flo		\
113	= (union _FP_UNION_E *) (val);			\
114							\
115      X##_f[2] = 0;					\
116      X##_f[3] = 0;					\
117      X##_f[0] = FP_UNPACK_RAW_EP_flo->bits.frac0;	\
118      X##_f[1] = FP_UNPACK_RAW_EP_flo->bits.frac1;	\
119      X##_e  = FP_UNPACK_RAW_EP_flo->bits.exp;		\
120      X##_s  = FP_UNPACK_RAW_EP_flo->bits.sign;		\
121    }							\
122  while (0)
123
124# define FP_PACK_RAW_E(val, X)			\
125  do						\
126    {						\
127      union _FP_UNION_E FP_PACK_RAW_E_flo;	\
128						\
129      if (X##_e)				\
130	X##_f[1] |= _FP_IMPLBIT_E;		\
131      else					\
132	X##_f[1] &= ~(_FP_IMPLBIT_E);		\
133      FP_PACK_RAW_E_flo.bits.frac0 = X##_f[0];	\
134      FP_PACK_RAW_E_flo.bits.frac1 = X##_f[1];	\
135      FP_PACK_RAW_E_flo.bits.exp   = X##_e;	\
136      FP_PACK_RAW_E_flo.bits.sign  = X##_s;	\
137						\
138      (val) = FP_PACK_RAW_E_flo.flt;		\
139    }						\
140  while (0)
141
142# define FP_PACK_RAW_EP(val, X)				\
143  do							\
144    {							\
145      if (!FP_INHIBIT_RESULTS)				\
146	{						\
147	  union _FP_UNION_E *FP_PACK_RAW_EP_flo		\
148	    = (union _FP_UNION_E *) (val);		\
149							\
150	  if (X##_e)					\
151	    X##_f[1] |= _FP_IMPLBIT_E;			\
152	  else						\
153	    X##_f[1] &= ~(_FP_IMPLBIT_E);		\
154	  FP_PACK_RAW_EP_flo->bits.frac0 = X##_f[0];	\
155	  FP_PACK_RAW_EP_flo->bits.frac1 = X##_f[1];	\
156	  FP_PACK_RAW_EP_flo->bits.exp   = X##_e;	\
157	  FP_PACK_RAW_EP_flo->bits.sign  = X##_s;	\
158	}						\
159    }							\
160  while (0)
161
162# define FP_UNPACK_E(X, val)			\
163  do						\
164    {						\
165      FP_UNPACK_RAW_E (X, (val));		\
166      _FP_UNPACK_CANONICAL (E, 4, X);		\
167    }						\
168  while (0)
169
170# define FP_UNPACK_EP(X, val)			\
171  do						\
172    {						\
173      FP_UNPACK_RAW_EP (X, (val));		\
174      _FP_UNPACK_CANONICAL (E, 4, X);		\
175    }						\
176  while (0)
177
178# define FP_UNPACK_SEMIRAW_E(X, val)		\
179  do						\
180    {						\
181      FP_UNPACK_RAW_E (X, (val));		\
182      _FP_UNPACK_SEMIRAW (E, 4, X);		\
183    }						\
184  while (0)
185
186# define FP_UNPACK_SEMIRAW_EP(X, val)		\
187  do						\
188    {						\
189      FP_UNPACK_RAW_EP (X, (val));		\
190      _FP_UNPACK_SEMIRAW (E, 4, X);		\
191    }						\
192  while (0)
193
194# define FP_PACK_E(val, X)			\
195  do						\
196    {						\
197      _FP_PACK_CANONICAL (E, 4, X);		\
198      FP_PACK_RAW_E ((val), X);			\
199    }						\
200  while (0)
201
202# define FP_PACK_EP(val, X)			\
203  do						\
204    {						\
205      _FP_PACK_CANONICAL (E, 4, X);		\
206      FP_PACK_RAW_EP ((val), X);		\
207    }						\
208  while (0)
209
210# define FP_PACK_SEMIRAW_E(val, X)		\
211  do						\
212    {						\
213      _FP_PACK_SEMIRAW (E, 4, X);		\
214      FP_PACK_RAW_E ((val), X);			\
215    }						\
216  while (0)
217
218# define FP_PACK_SEMIRAW_EP(val, X)		\
219  do						\
220    {						\
221      _FP_PACK_SEMIRAW (E, 4, X);		\
222      FP_PACK_RAW_EP ((val), X);		\
223    }						\
224  while (0)
225
226# define FP_ISSIGNAN_E(X)	_FP_ISSIGNAN (E, 4, X)
227# define FP_NEG_E(R, X)		_FP_NEG (E, 4, R, X)
228# define FP_ADD_E(R, X, Y)	_FP_ADD (E, 4, R, X, Y)
229# define FP_SUB_E(R, X, Y)	_FP_SUB (E, 4, R, X, Y)
230# define FP_MUL_E(R, X, Y)	_FP_MUL (E, 4, R, X, Y)
231# define FP_DIV_E(R, X, Y)	_FP_DIV (E, 4, R, X, Y)
232# define FP_SQRT_E(R, X)	_FP_SQRT (E, 4, R, X)
233# define FP_FMA_E(R, X, Y, Z)	_FP_FMA (E, 4, 8, R, X, Y, Z)
234
235/* Square root algorithms:
236   We have just one right now, maybe Newton approximation
237   should be added for those machines where division is fast.
238   This has special _E version because standard _4 square
239   root would not work (it has to start normally with the
240   second word and not the first), but as we have to do it
241   anyway, we optimize it by doing most of the calculations
242   in two UWtype registers instead of four.  */
243
244# define _FP_SQRT_MEAT_E(R, S, T, X, q)			\
245  do							\
246    {							\
247      (q) = (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE - 1);	\
248      _FP_FRAC_SRL_4 (X, (_FP_WORKBITS));		\
249      while (q)						\
250	{						\
251	  T##_f[1] = S##_f[1] + (q);			\
252	  if (T##_f[1] <= X##_f[1])			\
253	    {						\
254	      S##_f[1] = T##_f[1] + (q);		\
255	      X##_f[1] -= T##_f[1];			\
256	      R##_f[1] += (q);				\
257	    }						\
258	  _FP_FRAC_SLL_2 (X, 1);			\
259	  (q) >>= 1;					\
260	}						\
261      (q) = (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE - 1);	\
262      while (q)						\
263	{						\
264	  T##_f[0] = S##_f[0] + (q);			\
265	  T##_f[1] = S##_f[1];				\
266	  if (T##_f[1] < X##_f[1]			\
267	      || (T##_f[1] == X##_f[1]			\
268		  && T##_f[0] <= X##_f[0]))		\
269	    {						\
270	      S##_f[0] = T##_f[0] + (q);		\
271	      S##_f[1] += (T##_f[0] > S##_f[0]);	\
272	      _FP_FRAC_DEC_2 (X, T);			\
273	      R##_f[0] += (q);				\
274	    }						\
275	  _FP_FRAC_SLL_2 (X, 1);			\
276	  (q) >>= 1;					\
277	}						\
278      _FP_FRAC_SLL_4 (R, (_FP_WORKBITS));		\
279      if (X##_f[0] | X##_f[1])				\
280	{						\
281	  if (S##_f[1] < X##_f[1]			\
282	      || (S##_f[1] == X##_f[1]			\
283		  && S##_f[0] < X##_f[0]))		\
284	    R##_f[0] |= _FP_WORK_ROUND;			\
285	  R##_f[0] |= _FP_WORK_STICKY;			\
286	}						\
287    }							\
288  while (0)
289
290# define FP_CMP_E(r, X, Y, un, ex)	_FP_CMP (E, 4, (r), X, Y, (un), (ex))
291# define FP_CMP_EQ_E(r, X, Y, ex)	_FP_CMP_EQ (E, 4, (r), X, Y, (ex))
292# define FP_CMP_UNORD_E(r, X, Y, ex)	_FP_CMP_UNORD (E, 4, (r), X, Y, (ex))
293
294# define FP_TO_INT_E(r, X, rsz, rsg)	_FP_TO_INT (E, 4, (r), X, (rsz), (rsg))
295# define FP_FROM_INT_E(X, r, rs, rt)	_FP_FROM_INT (E, 4, X, (r), (rs), rt)
296
297# define _FP_FRAC_HIGH_E(X)	(X##_f[2])
298# define _FP_FRAC_HIGH_RAW_E(X)	(X##_f[1])
299
300# define _FP_FRAC_HIGH_DW_E(X)	(X##_f[4])
301
302#else   /* not _FP_W_TYPE_SIZE < 64 */
303union _FP_UNION_E
304{
305  XFtype flt;
306  struct _FP_STRUCT_LAYOUT
307  {
308# if __BYTE_ORDER == __BIG_ENDIAN
309    _FP_W_TYPE pad  : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
310    unsigned sign   : 1;
311    unsigned exp    : _FP_EXPBITS_E;
312    _FP_W_TYPE frac : _FP_W_TYPE_SIZE;
313# else
314    _FP_W_TYPE frac : _FP_W_TYPE_SIZE;
315    unsigned exp    : _FP_EXPBITS_E;
316    unsigned sign   : 1;
317# endif
318  } bits;
319};
320
321# define FP_DECL_E(X)		_FP_DECL (2, X)
322
323# define FP_UNPACK_RAW_E(X, val)		\
324  do						\
325    {						\
326      union _FP_UNION_E FP_UNPACK_RAW_E_flo;	\
327      FP_UNPACK_RAW_E_flo.flt = (val);		\
328						\
329      X##_f0 = FP_UNPACK_RAW_E_flo.bits.frac;	\
330      X##_f1 = 0;				\
331      X##_e = FP_UNPACK_RAW_E_flo.bits.exp;	\
332      X##_s = FP_UNPACK_RAW_E_flo.bits.sign;	\
333    }						\
334  while (0)
335
336# define FP_UNPACK_RAW_EP(X, val)		\
337  do						\
338    {						\
339      union _FP_UNION_E *FP_UNPACK_RAW_EP_flo	\
340	= (union _FP_UNION_E *) (val);		\
341						\
342      X##_f0 = FP_UNPACK_RAW_EP_flo->bits.frac;	\
343      X##_f1 = 0;				\
344      X##_e = FP_UNPACK_RAW_EP_flo->bits.exp;	\
345      X##_s = FP_UNPACK_RAW_EP_flo->bits.sign;	\
346    }						\
347  while (0)
348
349# define FP_PACK_RAW_E(val, X)			\
350  do						\
351    {						\
352      union _FP_UNION_E FP_PACK_RAW_E_flo;	\
353						\
354      if (X##_e)				\
355	X##_f0 |= _FP_IMPLBIT_E;		\
356      else					\
357	X##_f0 &= ~(_FP_IMPLBIT_E);		\
358      FP_PACK_RAW_E_flo.bits.frac = X##_f0;	\
359      FP_PACK_RAW_E_flo.bits.exp  = X##_e;	\
360      FP_PACK_RAW_E_flo.bits.sign = X##_s;	\
361						\
362      (val) = FP_PACK_RAW_E_flo.flt;		\
363    }						\
364  while (0)
365
366# define FP_PACK_RAW_EP(fs, val, X)			\
367  do							\
368    {							\
369      if (!FP_INHIBIT_RESULTS)				\
370	{						\
371	  union _FP_UNION_E *FP_PACK_RAW_EP_flo		\
372	    = (union _FP_UNION_E *) (val);		\
373							\
374	  if (X##_e)					\
375	    X##_f0 |= _FP_IMPLBIT_E;			\
376	  else						\
377	    X##_f0 &= ~(_FP_IMPLBIT_E);			\
378	  FP_PACK_RAW_EP_flo->bits.frac = X##_f0;	\
379	  FP_PACK_RAW_EP_flo->bits.exp  = X##_e;	\
380	  FP_PACK_RAW_EP_flo->bits.sign = X##_s;	\
381	}						\
382    }							\
383  while (0)
384
385
386# define FP_UNPACK_E(X, val)			\
387  do						\
388    {						\
389      FP_UNPACK_RAW_E (X, (val));		\
390      _FP_UNPACK_CANONICAL (E, 2, X);		\
391    }						\
392  while (0)
393
394# define FP_UNPACK_EP(X, val)			\
395  do						\
396    {						\
397      FP_UNPACK_RAW_EP (X, (val));		\
398      _FP_UNPACK_CANONICAL (E, 2, X);		\
399    }						\
400  while (0)
401
402# define FP_UNPACK_SEMIRAW_E(X, val)		\
403  do						\
404    {						\
405      FP_UNPACK_RAW_E (X, (val));		\
406      _FP_UNPACK_SEMIRAW (E, 2, X);		\
407    }						\
408  while (0)
409
410# define FP_UNPACK_SEMIRAW_EP(X, val)		\
411  do						\
412    {						\
413      FP_UNPACK_RAW_EP (X, (val));		\
414      _FP_UNPACK_SEMIRAW (E, 2, X);		\
415    }						\
416  while (0)
417
418# define FP_PACK_E(val, X)			\
419  do						\
420    {						\
421      _FP_PACK_CANONICAL (E, 2, X);		\
422      FP_PACK_RAW_E ((val), X);			\
423    }						\
424  while (0)
425
426# define FP_PACK_EP(val, X)			\
427  do						\
428    {						\
429      _FP_PACK_CANONICAL (E, 2, X);		\
430      FP_PACK_RAW_EP ((val), X);		\
431    }						\
432  while (0)
433
434# define FP_PACK_SEMIRAW_E(val, X)		\
435  do						\
436    {						\
437      _FP_PACK_SEMIRAW (E, 2, X);		\
438      FP_PACK_RAW_E ((val), X);			\
439    }						\
440  while (0)
441
442# define FP_PACK_SEMIRAW_EP(val, X)		\
443  do						\
444    {						\
445      _FP_PACK_SEMIRAW (E, 2, X);		\
446      FP_PACK_RAW_EP ((val), X);		\
447    }						\
448  while (0)
449
450# define FP_ISSIGNAN_E(X)	_FP_ISSIGNAN (E, 2, X)
451# define FP_NEG_E(R, X)		_FP_NEG (E, 2, R, X)
452# define FP_ADD_E(R, X, Y)	_FP_ADD (E, 2, R, X, Y)
453# define FP_SUB_E(R, X, Y)	_FP_SUB (E, 2, R, X, Y)
454# define FP_MUL_E(R, X, Y)	_FP_MUL (E, 2, R, X, Y)
455# define FP_DIV_E(R, X, Y)	_FP_DIV (E, 2, R, X, Y)
456# define FP_SQRT_E(R, X)	_FP_SQRT (E, 2, R, X)
457# define FP_FMA_E(R, X, Y, Z)	_FP_FMA (E, 2, 4, R, X, Y, Z)
458
459/* Square root algorithms:
460   We have just one right now, maybe Newton approximation
461   should be added for those machines where division is fast.
462   We optimize it by doing most of the calculations
463   in one UWtype registers instead of two, although we don't
464   have to.  */
465# define _FP_SQRT_MEAT_E(R, S, T, X, q)			\
466  do							\
467    {							\
468      (q) = (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE - 1);	\
469      _FP_FRAC_SRL_2 (X, (_FP_WORKBITS));		\
470      while (q)						\
471	{						\
472	  T##_f0 = S##_f0 + (q);			\
473	  if (T##_f0 <= X##_f0)				\
474	    {						\
475	      S##_f0 = T##_f0 + (q);			\
476	      X##_f0 -= T##_f0;				\
477	      R##_f0 += (q);				\
478	    }						\
479	  _FP_FRAC_SLL_1 (X, 1);			\
480	  (q) >>= 1;					\
481	}						\
482      _FP_FRAC_SLL_2 (R, (_FP_WORKBITS));		\
483      if (X##_f0)					\
484	{						\
485	  if (S##_f0 < X##_f0)				\
486	    R##_f0 |= _FP_WORK_ROUND;			\
487	  R##_f0 |= _FP_WORK_STICKY;			\
488	}						\
489    }							\
490  while (0)
491
492# define FP_CMP_E(r, X, Y, un, ex)	_FP_CMP (E, 2, (r), X, Y, (un), (ex))
493# define FP_CMP_EQ_E(r, X, Y, ex)	_FP_CMP_EQ (E, 2, (r), X, Y, (ex))
494# define FP_CMP_UNORD_E(r, X, Y, ex)	_FP_CMP_UNORD (E, 2, (r), X, Y, (ex))
495
496# define FP_TO_INT_E(r, X, rsz, rsg)	_FP_TO_INT (E, 2, (r), X, (rsz), (rsg))
497# define FP_FROM_INT_E(X, r, rs, rt)	_FP_FROM_INT (E, 2, X, (r), (rs), rt)
498
499# define _FP_FRAC_HIGH_E(X)	(X##_f1)
500# define _FP_FRAC_HIGH_RAW_E(X)	(X##_f0)
501
502# define _FP_FRAC_HIGH_DW_E(X)	(X##_f[2])
503
504#endif /* not _FP_W_TYPE_SIZE < 64 */
505