OSDN Git Service

* config/rs6000/darwin-ldouble.c: Build file for SOFT_FLOAT.
[pf3gnuchains/gcc-fork.git] / gcc / config / rs6000 / darwin-ldouble.c
1 /* 128-bit long double support routines for Darwin.
2    Copyright (C) 1993, 2003, 2004, 2005, 2006
3    Free Software Foundation, Inc.
4
5 This file is part of GCC.
6
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
10 version.
11
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
19 executable.)
20
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
24 for more details.
25
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
29 02110-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
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);
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
96 typedef union
97 {
98   long double ldval;
99   double dval[2];
100 } longDblUnion;
101
102 /* Add two 'long double' values and return the result.  */
103 long 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
142 long 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
149 static double fmsub (double, double, double);
150 #endif
151
152 long 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
185 long 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 #ifdef _SOFT_FLOAT
223
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);
238
239 extern int __eqdf2 (double, double);
240 extern int __ledf2 (double, double);
241 extern int __gedf2 (double, double);
242 extern int __unorddf2 (double, double);
243
244 /* Negate 'long double' value and return the result.    */
245 long 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.  */
256 int
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
264 strong_alias (__gcc_qeq, __gcc_qne);
265
266 /* Compare two 'long double' values for less than or equal.  */
267 int
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
275 strong_alias (__gcc_qle, __gcc_qlt);
276
277 /* Compare two 'long double' values for greater than or equal.  */
278 int
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
286 strong_alias (__gcc_qge, __gcc_qgt);
287
288 /* Compare two 'long double' values for unordered.  */
289 int
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.  */
298 long 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.  */
310 long 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.  */
322 float
323 __gcc_qtos (double a, double aa __attribute__ ((__unused__)))
324 {
325   return (float) a;
326 }
327
328 /* Convert long double to double.  */
329 double
330 __gcc_qtod (double a, double aa __attribute__ ((__unused__)))
331 {
332   return a;
333 }
334
335 /* Convert long double to int.  */
336 int
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.  */
344 unsigned 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.  */
352 long double
353 __gcc_itoq (int a)
354 {
355   return __gcc_dtoq ((double) a);
356 }
357
358 /* Convert unsigned int to long double.  */
359 long 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.  */
370 static double
371 fmsub (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