darwin-ldouble.c revision 132718 
1/* 128-bit long double support routines for Darwin.
2   Copyright (C) 1993, 2003, 2004 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__))
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 _xlqadd (double, double, double, double);
62extern long double _xlqsub (double, double, double, double);
63extern long double _xlqmul (double, double, double, double);
64extern long double _xlqdiv (double, double, double, double);
65
66typedef union
67{
68  long double ldval;
69  double dval[2];
70} longDblUnion;
71
72static const double FPKINF = 1.0/0.0;
73
74/* Add two 'long double' values and return the result.	*/
75long double
76_xlqadd (double a, double b, double c, double d)
77{
78  longDblUnion z;
79  double t, tau, u, FPR_zero, FPR_PosInf;
80
81  FPR_zero = 0.0;
82  FPR_PosInf = FPKINF;
83
84  if (unlikely (a != a) || unlikely (c != c))
85    return a + c;  /* NaN result.  */
86
87  /* Ordered operands are arranged in order of their magnitudes.  */
88
89  /* Switch inputs if |(c,d)| > |(a,b)|. */
90  if (fabs (c) > fabs (a))
91    {
92      t = a;
93      tau = b;
94      a = c;
95      b = d;
96      c = t;
97      d = tau;
98    }
99
100  /* b <- second largest magnitude double. */
101  if (fabs (c) > fabs (b))
102    {
103      t = b;
104      b = c;
105      c = t;
106    }
107
108  /* Thanks to commutivity, sum is invariant w.r.t. the next
109     conditional exchange. */
110  tau = d + c;
111
112  /* Order the smallest magnitude doubles.  */
113  if (fabs (d) > fabs (c))
114    {
115      t = c;
116      c = d;
117      d = t;
118    }
119
120  t = (tau + b) + a;	     /* Sum values in ascending magnitude order.  */
121
122  /* Infinite or zero result.  */
123  if (unlikely (t == FPR_zero) || unlikely (fabs (t) == FPR_PosInf))
124    return t;
125
126  /* Usual case.  */
127  tau = (((a-t) + b) + c) + d;
128  u = t + tau;
129  z.dval[0] = u;	       /* Final fixup for long double result.  */
130  z.dval[1] = (t - u) + tau;
131  return z.ldval;
132}
133
134long double
135_xlqsub (double a, double b, double c, double d)
136{
137  return _xlqadd (a, b, -c, -d);
138}
139
140long double
141_xlqmul (double a, double b, double c, double d)
142{
143  longDblUnion z;
144  double t, tau, u, v, w, FPR_zero, FPR_PosInf;
145
146  FPR_zero = 0.0;
147  FPR_PosInf = FPKINF;
148
149  t = a * c;			/* Highest order double term.  */
150
151  if (unlikely (t != t) || unlikely (t == FPR_zero)
152      || unlikely (fabs (t) == FPR_PosInf))
153    return t;
154
155  /* Finite nonzero result requires summing of terms of two highest
156     orders.	*/
157
158  /* Use fused multiply-add to get low part of a * c.	 */
159  asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
160  v = a*d;
161  w = b*c;
162  tau += v + w;	    /* Add in other second-order terms.	 */
163  u = t + tau;
164
165  /* Construct long double result.  */
166  z.dval[0] = u;
167  z.dval[1] = (t - u) + tau;
168  return z.ldval;
169}
170
171long double
172_xlqdiv (double a, double b, double c, double d)
173{
174  longDblUnion z;
175  double s, sigma, t, tau, u, v, w, FPR_zero, FPR_PosInf;
176
177  FPR_zero = 0.0;
178  FPR_PosInf = FPKINF;
179
180  t = a / c;                    /* highest order double term */
181
182  if (unlikely (t != t) || unlikely (t == FPR_zero)
183      || unlikely (fabs (t) == FPR_PosInf))
184    return t;
185
186  /* Finite nonzero result requires corrections to the highest order term.  */
187
188  s = c * t;                    /* (s,sigma) = c*t exactly. */
189  w = -(-b + d * t);	/* Written to get fnmsub for speed, but not
190			   numerically necessary.  */
191
192  /* Use fused multiply-add to get low part of c * t.	 */
193  asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
194  v = a - s;
195
196  tau = ((v-sigma)+w)/c;   /* Correction to t. */
197  u = t + tau;
198
199  /* Construct long double result. */
200  z.dval[0] = u;
201  z.dval[1] = (t - u) + tau;
202  return z.ldval;
203}
204
205#endif
206