darwin-ldouble.c revision 146895 
1/* 128-bit long double support routines for Darwin.
2   Copyright (C) 1993, 2003, 2004, 2005 Free Software Foundation, Inc.
3
4This file is part of GCC.
5
6GCC is free software; you can redistribute it and/or modify it under
7the terms of the GNU General Public License as published by the Free
8Software Foundation; either version 2, or (at your option) any later
9version.
10
11In addition to the permissions in the GNU General Public License, the
12Free Software Foundation gives you unlimited permission to link the
13compiled version of this file into combinations with other programs,
14and to distribute those combinations without any restriction coming
15from the use of this file.  (The General Public License restrictions
16do apply in other respects; for example, they cover modification of
17the file, and distribution when not linked into a combine
18executable.)
19
20GCC is distributed in the hope that it will be useful, but WITHOUT ANY
21WARRANTY; without even the implied warranty of MERCHANTABILITY or
22FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
23for more details.
24
25You should have received a copy of the GNU General Public License
26along with GCC; see the file COPYING.  If not, write to the Free
27Software Foundation, 59 Temple Place - Suite 330, Boston, MA
2802111-1307, USA.  */
29
30/* Implementations of floating-point long double basic arithmetic
31   functions called by the IBM C compiler when generating code for
32   PowerPC platforms.  In particular, the following functions are
33   implemented: _xlqadd, _xlqsub, _xlqmul, and _xlqdiv.  Double-double
34   algorithms are based on the paper "Doubled-Precision IEEE Standard
35   754 Floating-Point Arithmetic" by W. Kahan, February 26, 1987.  An
36   alternative published reference is "Software for Doubled-Precision
37   Floating-Point Computations", by Seppo Linnainmaa, ACM TOMS vol 7
38   no 3, September 1961, pages 272-283.  */
39
40/* Each long double is made up of two IEEE doubles.  The value of the
41   long double is the sum of the values of the two parts.  The most
42   significant part is required to be the value of the long double
43   rounded to the nearest double, as specified by IEEE.  For Inf
44   values, the least significant part is required to be one of +0.0 or
45   -0.0.  No other requirements are made; so, for example, 1.0 may be
46   represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
47   NaN is don't-care.
48
49   This code currently assumes big-endian.  */
50
51#if !_SOFT_FLOAT && (defined (__MACH__) || defined (__powerpc64__) || defined (_AIX))
52
53#define fabs(x) __builtin_fabs(x)
54
55#define unlikely(x) __builtin_expect ((x), 0)
56
57/* All these routines actually take two long doubles as parameters,
58   but GCC currently generates poor code when a union is used to turn
59   a long double into a pair of doubles.  */
60
61extern long double __gcc_qadd (double, double, double, double);
62extern long double __gcc_qsub (double, double, double, double);
63extern long double __gcc_qmul (double, double, double, double);
64extern long double __gcc_qdiv (double, double, double, double);
65
66#if defined __ELF__ && defined IN_LIBGCC2_S
67/* Provide definitions of the old symbol names to statisfy apps and
68   shared libs built against an older libgcc.  To access the _xlq
69   symbols an explicit version reference is needed, so these won't
70   satisfy an unadorned reference like _xlqadd.  If dot symbols are
71   not needed, the assembler will remove the aliases from the symbol
72   table.  */
73__asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
74         ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
75         ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
76         ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
77         ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
78         ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
79         ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
80         ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
81#endif
82
83typedef union
84{
85  long double ldval;
86  double dval[2];
87} longDblUnion;
88
89static const double FPKINF = 1.0/0.0;
90
91/* Add two 'long double' values and return the result.	*/
92long double
93__gcc_qadd (double a, double b, double c, double d)
94{
95  longDblUnion z;
96  double t, tau, u, FPR_zero, FPR_PosInf;
97
98  FPR_zero = 0.0;
99  FPR_PosInf = FPKINF;
100
101  if (unlikely (a != a) || unlikely (c != c))
102    return a + c;  /* NaN result.  */
103
104  /* Ordered operands are arranged in order of their magnitudes.  */
105
106  /* Switch inputs if |(c,d)| > |(a,b)|. */
107  if (fabs (c) > fabs (a))
108    {
109      t = a;
110      tau = b;
111      a = c;
112      b = d;
113      c = t;
114      d = tau;
115    }
116
117  /* b <- second largest magnitude double. */
118  if (fabs (c) > fabs (b))
119    {
120      t = b;
121      b = c;
122      c = t;
123    }
124
125  /* Thanks to commutivity, sum is invariant w.r.t. the next
126     conditional exchange. */
127  tau = d + c;
128
129  /* Order the smallest magnitude doubles.  */
130  if (fabs (d) > fabs (c))
131    {
132      t = c;
133      c = d;
134      d = t;
135    }
136
137  t = (tau + b) + a;	     /* Sum values in ascending magnitude order.  */
138
139  /* Infinite or zero result.  */
140  if (unlikely (t == FPR_zero) || unlikely (fabs (t) == FPR_PosInf))
141    return t;
142
143  /* Usual case.  */
144  tau = (((a-t) + b) + c) + d;
145  u = t + tau;
146  z.dval[0] = u;	       /* Final fixup for long double result.  */
147  z.dval[1] = (t - u) + tau;
148  return z.ldval;
149}
150
151long double
152__gcc_qsub (double a, double b, double c, double d)
153{
154  return __gcc_qadd (a, b, -c, -d);
155}
156
157long double
158__gcc_qmul (double a, double b, double c, double d)
159{
160  longDblUnion z;
161  double t, tau, u, v, w, FPR_zero, FPR_PosInf;
162
163  FPR_zero = 0.0;
164  FPR_PosInf = FPKINF;
165
166  t = a * c;			/* Highest order double term.  */
167
168  if (unlikely (t != t) || unlikely (t == FPR_zero)
169      || unlikely (fabs (t) == FPR_PosInf))
170    return t;
171
172  /* Finite nonzero result requires summing of terms of two highest
173     orders.	*/
174
175  /* Use fused multiply-add to get low part of a * c.	 */
176  asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
177  v = a*d;
178  w = b*c;
179  tau += v + w;	    /* Add in other second-order terms.	 */
180  u = t + tau;
181
182  /* Construct long double result.  */
183  z.dval[0] = u;
184  z.dval[1] = (t - u) + tau;
185  return z.ldval;
186}
187
188long double
189__gcc_qdiv (double a, double b, double c, double d)
190{
191  longDblUnion z;
192  double s, sigma, t, tau, u, v, w, FPR_zero, FPR_PosInf;
193
194  FPR_zero = 0.0;
195  FPR_PosInf = FPKINF;
196
197  t = a / c;                    /* highest order double term */
198
199  if (unlikely (t != t) || unlikely (t == FPR_zero)
200      || unlikely (fabs (t) == FPR_PosInf))
201    return t;
202
203  /* Finite nonzero result requires corrections to the highest order term.  */
204
205  s = c * t;                    /* (s,sigma) = c*t exactly. */
206  w = -(-b + d * t);	/* Written to get fnmsub for speed, but not
207			   numerically necessary.  */
208
209  /* Use fused multiply-add to get low part of c * t.	 */
210  asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
211  v = a - s;
212
213  tau = ((v-sigma)+w)/c;   /* Correction to t. */
214  u = t + tau;
215
216  /* Construct long double result. */
217  z.dval[0] = u;
218  z.dval[1] = (t - u) + tau;
219  return z.ldval;
220}
221
222#endif
223