1 /* 128-bit long double support routines for Darwin.
2 Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
3 Free Software Foundation, Inc.
5 This file is part of GCC.
7 GCC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU General Public License as published by the Free
9 Software Foundation; either version 2, or (at your option) any later
12 In addition to the permissions in the GNU General Public License, the
13 Free Software Foundation gives you unlimited permission to link the
14 compiled version of this file into combinations with other programs,
15 and to distribute those combinations without any restriction coming
16 from the use of this file. (The General Public License restrictions
17 do apply in other respects; for example, they cover modification of
18 the file, and distribution when not linked into a combine
21 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
22 WARRANTY; without even the implied warranty of MERCHANTABILITY or
23 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
26 You should have received a copy of the GNU General Public License
27 along with GCC; see the file COPYING. If not, write to the Free
28 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
31 /* Implementations of floating-point long double basic arithmetic
32 functions called by the IBM C compiler when generating code for
33 PowerPC platforms. In particular, the following functions are
34 implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
35 Double-double algorithms are based on the paper "Doubled-Precision
36 IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
37 1987. An alternative published reference is "Software for
38 Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
39 ACM TOMS vol 7 no 3, September 1981, pages 272-283. */
41 /* Each long double is made up of two IEEE doubles. The value of the
42 long double is the sum of the values of the two parts. The most
43 significant part is required to be the value of the long double
44 rounded to the nearest double, as specified by IEEE. For Inf
45 values, the least significant part is required to be one of +0.0 or
46 -0.0. No other requirements are made; so, for example, 1.0 may be
47 represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
50 This code currently assumes big-endian. */
52 #if ((!defined (__NO_FPRS__) || defined (_SOFT_FLOAT)) \
53 && !defined (__LITTLE_ENDIAN__) \
54 && (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))
56 #define fabs(x) __builtin_fabs(x)
57 #define isless(x, y) __builtin_isless (x, y)
58 #define inf() __builtin_inf()
60 #define unlikely(x) __builtin_expect ((x), 0)
62 #define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
64 /* Define ALIASNAME as a strong alias for NAME. */
65 # define strong_alias(name, aliasname) _strong_alias(name, aliasname)
66 # define _strong_alias(name, aliasname) \
67 extern __typeof (name) aliasname __attribute__ ((alias (#name)));
69 /* All these routines actually take two long doubles as parameters,
70 but GCC currently generates poor code when a union is used to turn
71 a long double into a pair of doubles. */
73 long double __gcc_qadd (double, double, double, double);
74 long double __gcc_qsub (double, double, double, double);
75 long double __gcc_qmul (double, double, double, double);
76 long double __gcc_qdiv (double, double, double, double);
78 #if defined __ELF__ && defined SHARED \
79 && (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))
80 /* Provide definitions of the old symbol names to satisfy apps and
81 shared libs built against an older libgcc. To access the _xlq
82 symbols an explicit version reference is needed, so these won't
83 satisfy an unadorned reference like _xlqadd. If dot symbols are
84 not needed, the assembler will remove the aliases from the symbol
86 __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
87 ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
88 ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
89 ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
90 ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
91 ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
92 ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
93 ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
102 /* Add two 'long double' values and return the result. */
104 __gcc_qadd (double a, double aa, double c, double cc)
116 x.dval[0] = z; /* Will always be DBL_MAX. */
118 if (fabs(a) > fabs(c))
119 x.dval[1] = a - z + c + zz;
121 x.dval[1] = c - z + a + zz;
126 zz = q + c + (a - (q + z)) + aa + cc;
128 /* Keep -0 result. */
137 x.dval[1] = z - xh + zz;
143 __gcc_qsub (double a, double b, double c, double d)
145 return __gcc_qadd (a, b, -c, -d);
149 static double fmsub (double, double, double);
153 __gcc_qmul (double a, double b, double c, double d)
156 double t, tau, u, v, w;
158 t = a * c; /* Highest order double term. */
160 if (unlikely (t == 0) /* Preserve -0. */
164 /* Sum terms of two highest orders. */
166 /* Use fused multiply-add to get low part of a * c. */
168 asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
170 tau = fmsub (a, c, t);
174 tau += v + w; /* Add in other second-order terms. */
177 /* Construct long double result. */
181 z.dval[1] = (t - u) + tau;
186 __gcc_qdiv (double a, double b, double c, double d)
189 double s, sigma, t, tau, u, v, w;
191 t = a / c; /* highest order double term */
193 if (unlikely (t == 0) /* Preserve -0. */
197 /* Finite nonzero result requires corrections to the highest order term. */
199 s = c * t; /* (s,sigma) = c*t exactly. */
200 w = -(-b + d * t); /* Written to get fnmsub for speed, but not
201 numerically necessary. */
203 /* Use fused multiply-add to get low part of c * t. */
205 asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
207 sigma = fmsub (c, t, s);
211 tau = ((v-sigma)+w)/c; /* Correction to t. */
214 /* Construct long double result. */
218 z.dval[1] = (t - u) + tau;
222 #if defined (_SOFT_FLOAT) && defined (__LONG_DOUBLE_128__)
224 long double __gcc_qneg (double, double);
225 int __gcc_qeq (double, double, double, double);
226 int __gcc_qne (double, double, double, double);
227 int __gcc_qge (double, double, double, double);
228 int __gcc_qle (double, double, double, double);
229 int __gcc_qunord (double, double, double, double);
230 long double __gcc_stoq (float);
231 long double __gcc_dtoq (double);
232 float __gcc_qtos (double, double);
233 double __gcc_qtod (double, double);
234 int __gcc_qtoi (double, double);
235 unsigned int __gcc_qtou (double, double);
236 long double __gcc_itoq (int);
237 long double __gcc_utoq (unsigned int);
239 extern int __eqdf2 (double, double);
240 extern int __ledf2 (double, double);
241 extern int __gedf2 (double, double);
242 extern int __unorddf2 (double, double);
244 /* Negate 'long double' value and return the result. */
246 __gcc_qneg (double a, double aa)
255 /* Compare two 'long double' values for equality. */
257 __gcc_qeq (double a, double aa, double c, double cc)
259 if (__eqdf2 (a, c) == 0)
260 return __eqdf2 (aa, cc);
264 strong_alias (__gcc_qeq, __gcc_qne);
266 /* Compare two 'long double' values for less than or equal. */
268 __gcc_qle (double a, double aa, double c, double cc)
270 if (__eqdf2 (a, c) == 0)
271 return __ledf2 (aa, cc);
272 return __ledf2 (a, c);
275 strong_alias (__gcc_qle, __gcc_qlt);
277 /* Compare two 'long double' values for greater than or equal. */
279 __gcc_qge (double a, double aa, double c, double cc)
281 if (__eqdf2 (a, c) == 0)
282 return __gedf2 (aa, cc);
283 return __gedf2 (a, c);
286 strong_alias (__gcc_qge, __gcc_qgt);
288 /* Compare two 'long double' values for unordered. */
290 __gcc_qunord (double a, double aa, double c, double cc)
292 if (__eqdf2 (a, c) == 0)
293 return __unorddf2 (aa, cc);
294 return __unorddf2 (a, c);
297 /* Convert single to long double. */
303 x.dval[0] = (double) a;
309 /* Convert double to long double. */
311 __gcc_dtoq (double a)
321 /* Convert long double to single. */
323 __gcc_qtos (double a, double aa __attribute__ ((__unused__)))
328 /* Convert long double to double. */
330 __gcc_qtod (double a, double aa __attribute__ ((__unused__)))
335 /* Convert long double to int. */
337 __gcc_qtoi (double a, double aa)
343 /* Convert long double to unsigned int. */
345 __gcc_qtou (double a, double aa)
348 return (unsigned int) z;
351 /* Convert int to long double. */
355 return __gcc_dtoq ((double) a);
358 /* Convert unsigned int to long double. */
360 __gcc_utoq (unsigned int a)
362 return __gcc_dtoq ((double) a);
365 #include "config/soft-fp/soft-fp.h"
366 #include "config/soft-fp/double.h"
367 #include "config/soft-fp/quad.h"
369 /* Compute floating point multiply-subtract with higher (quad) precision. */
371 fmsub (double a, double b, double c)
384 long double u, v, x, y, z;
387 FP_UNPACK_RAW_D (A, a);
388 FP_UNPACK_RAW_D (B, b);
389 FP_UNPACK_RAW_D (C, c);
391 /* Extend double to quad. */
392 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
393 FP_EXTEND(Q,D,4,2,X,A);
394 FP_EXTEND(Q,D,4,2,Y,B);
395 FP_EXTEND(Q,D,4,2,Z,C);
397 FP_EXTEND(Q,D,2,1,X,A);
398 FP_EXTEND(Q,D,2,1,Y,B);
399 FP_EXTEND(Q,D,2,1,Z,C);
404 FP_HANDLE_EXCEPTIONS;
412 FP_HANDLE_EXCEPTIONS;
416 FP_UNPACK_SEMIRAW_Q(U,u);
417 FP_UNPACK_SEMIRAW_Q(Z,z);
419 FP_PACK_SEMIRAW_Q(v,V);
420 FP_HANDLE_EXCEPTIONS;
422 /* Truncate quad to double. */
424 FP_UNPACK_SEMIRAW_Q(V,v);
425 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
426 FP_TRUNC(D,Q,2,4,R,V);
428 FP_TRUNC(D,Q,1,2,R,V);
430 FP_PACK_SEMIRAW_D(r,R);
431 FP_HANDLE_EXCEPTIONS;