1/* 128-bit long double support routines for Darwin.
2   Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
3   Free Software Foundation, Inc.
4
5This file is part of GCC.
6
7GCC is free software; you can redistribute it and/or modify it under
8the terms of the GNU General Public License as published by the Free
9Software Foundation; either version 2, or (at your option) any later
10version.
11
12In addition to the permissions in the GNU General Public License, the
13Free Software Foundation gives you unlimited permission to link the
14compiled version of this file into combinations with other programs,
15and to distribute those combinations without any restriction coming
16from the use of this file.  (The General Public License restrictions
17do apply in other respects; for example, they cover modification of
18the file, and distribution when not linked into a combine
19executable.)
20
21GCC is distributed in the hope that it will be useful, but WITHOUT ANY
22WARRANTY; without even the implied warranty of MERCHANTABILITY or
23FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
24for more details.
25
26You should have received a copy of the GNU General Public License
27along with GCC; see the file COPYING.  If not, write to the Free
28Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
2902110-1301, USA.  */
30
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.  */
40
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
48   NaN is don't-care.
49
50   This code currently assumes big-endian.  */
51
52#if ((!defined (__NO_FPRS__) || defined (_SOFT_FLOAT)) \
53     && !defined (__LITTLE_ENDIAN__) \
54     && (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))
55
56#define fabs(x) __builtin_fabs(x)
57#define isless(x, y) __builtin_isless (x, y)
58#define inf() __builtin_inf()
59
60#define unlikely(x) __builtin_expect ((x), 0)
61
62#define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
63
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)));
68
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.  */
72
73long double __gcc_qadd (double, double, double, double);
74long double __gcc_qsub (double, double, double, double);
75long double __gcc_qmul (double, double, double, double);
76long double __gcc_qdiv (double, double, double, double);
77
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
85   table.  */
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");
94#endif
95
96typedef union
97{
98  long double ldval;
99  double dval[2];
100} longDblUnion;
101
102/* Add two 'long double' values and return the result.	*/
103long double
104__gcc_qadd (double a, double aa, double c, double cc)
105{
106  longDblUnion x;
107  double z, q, zz, xh;
108
109  z = a + c;
110
111  if (nonfinite (z))
112    {
113      z = cc + aa + c + a;
114      if (nonfinite (z))
115	return z;
116      x.dval[0] = z;  /* Will always be DBL_MAX.  */
117      zz = aa + cc;
118      if (fabs(a) > fabs(c))
119	x.dval[1] = a - z + c + zz;
120      else
121	x.dval[1] = c - z + a + zz;
122    }
123  else
124    {
125      q = a - z;
126      zz = q + c + (a - (q + z)) + aa + cc;
127
128      /* Keep -0 result.  */
129      if (zz == 0.0)
130	return z;
131
132      xh = z + zz;
133      if (nonfinite (xh))
134	return xh;
135
136      x.dval[0] = xh;
137      x.dval[1] = z - xh + zz;
138    }
139  return x.ldval;
140}
141
142long double
143__gcc_qsub (double a, double b, double c, double d)
144{
145  return __gcc_qadd (a, b, -c, -d);
146}
147
148#ifdef _SOFT_FLOAT
149static double fmsub (double, double, double);
150#endif
151
152long double
153__gcc_qmul (double a, double b, double c, double d)
154{
155  longDblUnion z;
156  double t, tau, u, v, w;
157
158  t = a * c;			/* Highest order double term.  */
159
160  if (unlikely (t == 0)		/* Preserve -0.  */
161      || nonfinite (t))
162    return t;
163
164  /* Sum terms of two highest orders. */
165
166  /* Use fused multiply-add to get low part of a * c.  */
167#ifndef _SOFT_FLOAT
168  asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
169#else
170  tau = fmsub (a, c, t);
171#endif
172  v = a*d;
173  w = b*c;
174  tau += v + w;	    /* Add in other second-order terms.	 */
175  u = t + tau;
176
177  /* Construct long double result.  */
178  if (nonfinite (u))
179    return u;
180  z.dval[0] = u;
181  z.dval[1] = (t - u) + tau;
182  return z.ldval;
183}
184
185long double
186__gcc_qdiv (double a, double b, double c, double d)
187{
188  longDblUnion z;
189  double s, sigma, t, tau, u, v, w;
190
191  t = a / c;                    /* highest order double term */
192
193  if (unlikely (t == 0)		/* Preserve -0.  */
194      || nonfinite (t))
195    return t;
196
197  /* Finite nonzero result requires corrections to the highest order term.  */
198
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.  */
202
203  /* Use fused multiply-add to get low part of c * t.	 */
204#ifndef _SOFT_FLOAT
205  asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
206#else
207  sigma = fmsub (c, t, s);
208#endif
209  v = a - s;
210
211  tau = ((v-sigma)+w)/c;   /* Correction to t.  */
212  u = t + tau;
213
214  /* Construct long double result.  */
215  if (nonfinite (u))
216    return u;
217  z.dval[0] = u;
218  z.dval[1] = (t - u) + tau;
219  return z.ldval;
220}
221
222#if defined (_SOFT_FLOAT) && defined (__LONG_DOUBLE_128__)
223
224long double __gcc_qneg (double, double);
225int __gcc_qeq (double, double, double, double);
226int __gcc_qne (double, double, double, double);
227int __gcc_qge (double, double, double, double);
228int __gcc_qle (double, double, double, double);
229int __gcc_qunord (double, double, double, double);
230long double __gcc_stoq (float);
231long double __gcc_dtoq (double);
232float __gcc_qtos (double, double);
233double __gcc_qtod (double, double);
234int __gcc_qtoi (double, double);
235unsigned int __gcc_qtou (double, double);
236long double __gcc_itoq (int);
237long double __gcc_utoq (unsigned int);
238
239extern int __eqdf2 (double, double);
240extern int __ledf2 (double, double);
241extern int __gedf2 (double, double);
242extern int __unorddf2 (double, double);
243
244/* Negate 'long double' value and return the result.	*/
245long double
246__gcc_qneg (double a, double aa)
247{
248  longDblUnion x;
249
250  x.dval[0] = -a;
251  x.dval[1] = -aa;
252  return x.ldval;
253}
254
255/* Compare two 'long double' values for equality.  */
256int
257__gcc_qeq (double a, double aa, double c, double cc)
258{
259  if (__eqdf2 (a, c) == 0)
260    return __eqdf2 (aa, cc);
261  return 1;
262}
263
264strong_alias (__gcc_qeq, __gcc_qne);
265
266/* Compare two 'long double' values for less than or equal.  */
267int
268__gcc_qle (double a, double aa, double c, double cc)
269{
270  if (__eqdf2 (a, c) == 0)
271    return __ledf2 (aa, cc);
272  return __ledf2 (a, c);
273}
274
275strong_alias (__gcc_qle, __gcc_qlt);
276
277/* Compare two 'long double' values for greater than or equal.  */
278int
279__gcc_qge (double a, double aa, double c, double cc)
280{
281  if (__eqdf2 (a, c) == 0)
282    return __gedf2 (aa, cc);
283  return __gedf2 (a, c);
284}
285
286strong_alias (__gcc_qge, __gcc_qgt);
287
288/* Compare two 'long double' values for unordered.  */
289int
290__gcc_qunord (double a, double aa, double c, double cc)
291{
292  if (__eqdf2 (a, c) == 0)
293    return __unorddf2 (aa, cc);
294  return __unorddf2 (a, c);
295}
296
297/* Convert single to long double.  */
298long double
299__gcc_stoq (float a)
300{
301  longDblUnion x;
302
303  x.dval[0] = (double) a;
304  x.dval[1] = 0.0;
305
306  return x.ldval;
307}
308
309/* Convert double to long double.  */
310long double
311__gcc_dtoq (double a)
312{
313  longDblUnion x;
314
315  x.dval[0] = a;
316  x.dval[1] = 0.0;
317
318  return x.ldval;
319}
320
321/* Convert long double to single.  */
322float
323__gcc_qtos (double a, double aa __attribute__ ((__unused__)))
324{
325  return (float) a;
326}
327
328/* Convert long double to double.  */
329double
330__gcc_qtod (double a, double aa __attribute__ ((__unused__)))
331{
332  return a;
333}
334
335/* Convert long double to int.  */
336int
337__gcc_qtoi (double a, double aa)
338{
339  double z = a + aa;
340  return (int) z;
341}
342
343/* Convert long double to unsigned int.  */
344unsigned int
345__gcc_qtou (double a, double aa)
346{
347  double z = a + aa;
348  return (unsigned int) z;
349}
350
351/* Convert int to long double.  */
352long double
353__gcc_itoq (int a)
354{
355  return __gcc_dtoq ((double) a);
356}
357
358/* Convert unsigned int to long double.  */
359long double
360__gcc_utoq (unsigned int a)
361{
362  return __gcc_dtoq ((double) a);
363}
364
365#include "config/soft-fp/soft-fp.h"
366#include "config/soft-fp/double.h"
367#include "config/soft-fp/quad.h"
368
369/* Compute floating point multiply-subtract with higher (quad) precision.  */
370static double
371fmsub (double a, double b, double c)
372{
373    FP_DECL_EX;
374    FP_DECL_D(A);
375    FP_DECL_D(B);
376    FP_DECL_D(C);
377    FP_DECL_Q(X);
378    FP_DECL_Q(Y);
379    FP_DECL_Q(Z);
380    FP_DECL_Q(U);
381    FP_DECL_Q(V);
382    FP_DECL_D(R);
383    double r;
384    long double u, v, x, y, z;
385
386    FP_INIT_ROUNDMODE;
387    FP_UNPACK_RAW_D (A, a);
388    FP_UNPACK_RAW_D (B, b);
389    FP_UNPACK_RAW_D (C, c);
390
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);
396#else
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);
400#endif
401    FP_PACK_RAW_Q(x,X);
402    FP_PACK_RAW_Q(y,Y);
403    FP_PACK_RAW_Q(z,Z);
404    FP_HANDLE_EXCEPTIONS;
405
406    /* Multiply.  */
407    FP_INIT_ROUNDMODE;
408    FP_UNPACK_Q(X,x);
409    FP_UNPACK_Q(Y,y);
410    FP_MUL_Q(U,X,Y);
411    FP_PACK_Q(u,U);
412    FP_HANDLE_EXCEPTIONS;
413
414    /* Subtract.  */
415    FP_INIT_ROUNDMODE;
416    FP_UNPACK_SEMIRAW_Q(U,u);
417    FP_UNPACK_SEMIRAW_Q(Z,z);
418    FP_SUB_Q(V,U,Z);
419    FP_PACK_SEMIRAW_Q(v,V);
420    FP_HANDLE_EXCEPTIONS;
421
422    /* Truncate quad to double.  */
423    FP_INIT_ROUNDMODE;
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);
427#else
428    FP_TRUNC(D,Q,1,2,R,V);
429#endif
430    FP_PACK_SEMIRAW_D(r,R);
431    FP_HANDLE_EXCEPTIONS;
432
433    return r;
434}
435
436#endif
437
438#endif
439