1/* GCC Quad-Precision Math Library
2   Copyright (C) 2010, 2011 Free Software Foundation, Inc.
3   Written by Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
4
5This file is part of the libquadmath library.
6Libquadmath is free software; you can redistribute it and/or
7modify it under the terms of the GNU Library General Public
8License as published by the Free Software Foundation; either
9version 2 of the License, or (at your option) any later version.
10
11Libquadmath is distributed in the hope that it will be useful,
12but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14Library General Public License for more details.
15
16You should have received a copy of the GNU Library General Public
17License along with libquadmath; see the file COPYING.LIB.  If
18not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19Boston, MA 02110-1301, USA.  */
20
21#ifndef QUADMATH_IMP_H
22#define QUADMATH_IMP_H
23
24#include <errno.h>
25#include <limits.h>
26#include <stdbool.h>
27#include <stdint.h>
28#include <stdlib.h>
29#include "quadmath.h"
30#include "config.h"
31#ifdef HAVE_FENV_H
32# include <fenv.h>
33#endif
34
35
36/* Under IEEE 754, an architecture may determine tininess of
37   floating-point results either "before rounding" or "after
38   rounding", but must do so in the same way for all operations
39   returning binary results.  Define TININESS_AFTER_ROUNDING to 1 for
40   "after rounding" architectures, 0 for "before rounding"
41   architectures.  */
42
43#define TININESS_AFTER_ROUNDING   1
44
45#define HIGH_ORDER_BIT_IS_SET_FOR_SNAN 0
46
47#define FIX_FLT128_LONG_CONVERT_OVERFLOW 0
48#define FIX_FLT128_LLONG_CONVERT_OVERFLOW 0
49
50/* Prototypes for internal functions.  */
51extern int32_t __quadmath_rem_pio2q (__float128, __float128 *);
52extern void __quadmath_kernel_sincosq (__float128, __float128, __float128 *,
53				       __float128 *, int);
54extern __float128 __quadmath_kernel_sinq (__float128, __float128, int);
55extern __float128 __quadmath_kernel_cosq (__float128, __float128);
56extern __float128 __quadmath_kernel_tanq (__float128, __float128, int);
57extern __float128 __quadmath_gamma_productq (__float128, __float128, int,
58					     __float128 *);
59extern __float128 __quadmath_gammaq_r (__float128, int *);
60extern __float128 __quadmath_lgamma_negq (__float128, int *);
61extern __float128 __quadmath_lgamma_productq (__float128, __float128,
62					      __float128, int);
63extern __float128 __quadmath_lgammaq_r (__float128, int *);
64extern __float128 __quadmath_x2y2m1q (__float128 x, __float128 y);
65extern __complex128 __quadmath_kernel_casinhq (__complex128, int);
66
67static inline void
68mul_splitq (__float128 *hi, __float128 *lo, __float128 x, __float128 y)
69{
70  /* Fast built-in fused multiply-add.  */
71  *hi = x * y;
72  *lo = fmaq (x, y, -*hi);
73}
74
75
76
77
78/* Frankly, if you have __float128, you have 64-bit integers, right?  */
79#ifndef UINT64_C
80# error "No way!"
81#endif
82
83
84/* Main union type we use to manipulate the floating-point type.  */
85typedef union
86{
87  __float128 value;
88
89  struct
90#ifdef __MINGW32__
91  /* On mingw targets the ms-bitfields option is active by default.
92     Therefore enforce gnu-bitfield style.  */
93  __attribute__ ((gcc_struct))
94#endif
95  {
96#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
97    unsigned negative:1;
98    unsigned exponent:15;
99    unsigned mantissa0:16;
100    unsigned mantissa1:32;
101    unsigned mantissa2:32;
102    unsigned mantissa3:32;
103#else
104    unsigned mantissa3:32;
105    unsigned mantissa2:32;
106    unsigned mantissa1:32;
107    unsigned mantissa0:16;
108    unsigned exponent:15;
109    unsigned negative:1;
110#endif
111  } ieee;
112
113  struct
114  {
115#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
116    uint64_t high;
117    uint64_t low;
118#else
119    uint64_t low;
120    uint64_t high;
121#endif
122  } words64;
123
124  struct
125  {
126#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
127    uint32_t w0;
128    uint32_t w1;
129    uint32_t w2;
130    uint32_t w3;
131#else
132    uint32_t w3;
133    uint32_t w2;
134    uint32_t w1;
135    uint32_t w0;
136#endif
137  } words32;
138
139  struct
140#ifdef __MINGW32__
141  /* Make sure we are using gnu-style bitfield handling.  */
142  __attribute__ ((gcc_struct))
143#endif
144  {
145#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
146    unsigned negative:1;
147    unsigned exponent:15;
148    unsigned quiet_nan:1;
149    unsigned mantissa0:15;
150    unsigned mantissa1:32;
151    unsigned mantissa2:32;
152    unsigned mantissa3:32;
153#else
154    unsigned mantissa3:32;
155    unsigned mantissa2:32;
156    unsigned mantissa1:32;
157    unsigned mantissa0:15;
158    unsigned quiet_nan:1;
159    unsigned exponent:15;
160    unsigned negative:1;
161#endif
162  } ieee_nan;
163
164} ieee854_float128;
165
166
167/* Get two 64 bit ints from a long double.  */
168#define GET_FLT128_WORDS64(ix0,ix1,d)  \
169do {                                   \
170  ieee854_float128 u;                  \
171  u.value = (d);                       \
172  (ix0) = u.words64.high;              \
173  (ix1) = u.words64.low;               \
174} while (0)
175
176/* Set a long double from two 64 bit ints.  */
177#define SET_FLT128_WORDS64(d,ix0,ix1)  \
178do {                                   \
179  ieee854_float128 u;                  \
180  u.words64.high = (ix0);              \
181  u.words64.low = (ix1);               \
182  (d) = u.value;                       \
183} while (0)
184
185/* Get the more significant 64 bits of a long double mantissa.  */
186#define GET_FLT128_MSW64(v,d)          \
187do {                                   \
188  ieee854_float128 u;                  \
189  u.value = (d);                       \
190  (v) = u.words64.high;                \
191} while (0)
192
193/* Set the more significant 64 bits of a long double mantissa from an int.  */
194#define SET_FLT128_MSW64(d,v)          \
195do {                                   \
196  ieee854_float128 u;                  \
197  u.value = (d);                       \
198  u.words64.high = (v);                \
199  (d) = u.value;                       \
200} while (0)
201
202/* Get the least significant 64 bits of a long double mantissa.  */
203#define GET_FLT128_LSW64(v,d)          \
204do {                                   \
205  ieee854_float128 u;                  \
206  u.value = (d);                       \
207  (v) = u.words64.low;                 \
208} while (0)
209
210
211#define IEEE854_FLOAT128_BIAS 0x3fff
212
213#define QUADFP_NAN		0
214#define QUADFP_INFINITE		1
215#define QUADFP_ZERO		2
216#define QUADFP_SUBNORMAL	3
217#define QUADFP_NORMAL		4
218#define fpclassifyq(x) \
219  __builtin_fpclassify (QUADFP_NAN, QUADFP_INFINITE, QUADFP_NORMAL, \
220			QUADFP_SUBNORMAL, QUADFP_ZERO, x)
221
222#ifndef math_opt_barrier
223# define math_opt_barrier(x) \
224({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; })
225# define math_force_eval(x) \
226({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); })
227#endif
228
229/* math_narrow_eval reduces its floating-point argument to the range
230   and precision of its semantic type.  (The original evaluation may
231   still occur with excess range and precision, so the result may be
232   affected by double rounding.)  */
233#define math_narrow_eval(x) (x)
234
235/* If X (which is not a NaN) is subnormal, force an underflow
236   exception.  */
237#define math_check_force_underflow(x)				\
238  do								\
239    {								\
240      __float128 force_underflow_tmp = (x);			\
241      if (fabsq (force_underflow_tmp) < FLT128_MIN)		\
242	{							\
243	  __float128 force_underflow_tmp2			\
244	    = force_underflow_tmp * force_underflow_tmp;	\
245	  math_force_eval (force_underflow_tmp2);		\
246	}							\
247    }								\
248  while (0)
249/* Likewise, but X is also known to be nonnegative.  */
250#define math_check_force_underflow_nonneg(x)			\
251  do								\
252    {								\
253      __float128 force_underflow_tmp = (x);			\
254      if (force_underflow_tmp < FLT128_MIN)			\
255	{							\
256	  __float128 force_underflow_tmp2			\
257	    = force_underflow_tmp * force_underflow_tmp;	\
258	  math_force_eval (force_underflow_tmp2);		\
259	}							\
260    }								\
261  while (0)
262
263/* Likewise, for both real and imaginary parts of a complex
264   result.  */
265#define math_check_force_underflow_complex(x)				\
266  do									\
267    {									\
268      __typeof (x) force_underflow_complex_tmp = (x);			\
269      math_check_force_underflow (__real__ force_underflow_complex_tmp); \
270      math_check_force_underflow (__imag__ force_underflow_complex_tmp); \
271    }									\
272  while (0)
273
274#ifndef HAVE_FENV_H
275# define feraiseexcept(arg) ((void) 0)
276typedef int fenv_t;
277# define feholdexcept(arg) ((void) 0)
278# define fesetround(arg) ((void) 0)
279# define feupdateenv(arg) ((void) (arg))
280# define fesetenv(arg) ((void) (arg))
281# define fetestexcept(arg) 0
282# define feclearexcept(arg) ((void) 0)
283#else
284# ifndef HAVE_FEHOLDEXCEPT
285#  define feholdexcept(arg) ((void) 0)
286# endif
287# ifndef HAVE_FESETROUND
288#  define fesetround(arg) ((void) 0)
289# endif
290# ifndef HAVE_FEUPDATEENV
291#  define feupdateenv(arg) ((void) (arg))
292# endif
293# ifndef HAVE_FESETENV
294#  define fesetenv(arg) ((void) (arg))
295# endif
296# ifndef HAVE_FETESTEXCEPT
297#  define fetestexcept(arg) 0
298# endif
299#endif
300
301#ifndef __glibc_likely
302# define __glibc_likely(cond)	__builtin_expect ((cond), 1)
303#endif
304
305#ifndef __glibc_unlikely
306# define __glibc_unlikely(cond)	__builtin_expect ((cond), 0)
307#endif
308
309#if defined HAVE_FENV_H && defined HAVE_FESETROUND && defined HAVE_FEUPDATEENV
310struct rm_ctx
311{
312  fenv_t env;
313  bool updated_status;
314};
315
316# define SET_RESTORE_ROUNDF128(RM)					\
317  struct rm_ctx ctx __attribute__((cleanup (libc_feresetround_ctx)));	\
318  libc_feholdsetround_ctx (&ctx, (RM))
319
320static inline __attribute__ ((always_inline)) void
321libc_feholdsetround_ctx (struct rm_ctx *ctx, int round)
322{
323  ctx->updated_status = false;
324
325  /* Update rounding mode only if different.  */
326  if (__glibc_unlikely (round != fegetround ()))
327    {
328      ctx->updated_status = true;
329      fegetenv (&ctx->env);
330      fesetround (round);
331    }
332}
333
334static inline __attribute__ ((always_inline)) void
335libc_feresetround_ctx (struct rm_ctx *ctx)
336{
337  /* Restore the rounding mode if updated.  */
338  if (__glibc_unlikely (ctx->updated_status))
339    feupdateenv (&ctx->env);
340}
341#else
342# define SET_RESTORE_ROUNDF128(RM) ((void) 0)
343#endif
344
345#endif
346