Mercurial > hg > CbC > CbC_gcc

comparison gcc/config/rs6000/darwin-ldouble.c @ 0:a06113de4d67

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first commit
author kent <kent@cr.ie.u-ryukyu.ac.jp>
date Fri, 17 Jul 2009 14:47:48 +0900
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a06113de4d67
1 /* 128-bit long double support routines for Darwin.
2 Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007, 2008, 2009
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 3, or (at your option) any later
10 version.
11
12 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 Under Section 7 of GPL version 3, you are granted additional
18 permissions described in the GCC Runtime Library Exception, version
19 3.1, as published by the Free Software Foundation.
20
21 You should have received a copy of the GNU General Public License and
22 a copy of the GCC Runtime Library Exception along with this program;
23 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
24 <http://www.gnu.org/licenses/>. */
25
26
27 /* Implementations of floating-point long double basic arithmetic
28 functions called by the IBM C compiler when generating code for
29 PowerPC platforms. In particular, the following functions are
30 implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
31 Double-double algorithms are based on the paper "Doubled-Precision
32 IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
33 1987. An alternative published reference is "Software for
34 Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
35 ACM TOMS vol 7 no 3, September 1981, pages 272-283. */
36
37 /* Each long double is made up of two IEEE doubles. The value of the
38 long double is the sum of the values of the two parts. The most
39 significant part is required to be the value of the long double
40 rounded to the nearest double, as specified by IEEE. For Inf
41 values, the least significant part is required to be one of +0.0 or
42 -0.0. No other requirements are made; so, for example, 1.0 may be
43 represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
44 NaN is don't-care.
45
46 This code currently assumes big-endian. */
47
48 #if (!defined (__LITTLE_ENDIAN__) \
49 && (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))
50
51 #define fabs(x) __builtin_fabs(x)
52 #define isless(x, y) __builtin_isless (x, y)
53 #define inf() __builtin_inf()
54
55 #define unlikely(x) __builtin_expect ((x), 0)
56
57 #define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
58
59 /* Define ALIASNAME as a strong alias for NAME. */
60 # define strong_alias(name, aliasname) _strong_alias(name, aliasname)
61 # define _strong_alias(name, aliasname) \
62 extern __typeof (name) aliasname __attribute__ ((alias (#name)));
63
64 /* All these routines actually take two long doubles as parameters,
65 but GCC currently generates poor code when a union is used to turn
66 a long double into a pair of doubles. */
67
68 long double __gcc_qadd (double, double, double, double);
69 long double __gcc_qsub (double, double, double, double);
70 long double __gcc_qmul (double, double, double, double);
71 long double __gcc_qdiv (double, double, double, double);
72
73 #if defined __ELF__ && defined SHARED \
74 && (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))
75 /* Provide definitions of the old symbol names to satisfy apps and
76 shared libs built against an older libgcc. To access the _xlq
77 symbols an explicit version reference is needed, so these won't
78 satisfy an unadorned reference like _xlqadd. If dot symbols are
79 not needed, the assembler will remove the aliases from the symbol
80 table. */
81 __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
82 ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
83 ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
84 ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
85 ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
86 ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
87 ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
88 ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
89 #endif
90
91 typedef union
92 {
93 long double ldval;
94 double dval[2];
95 } longDblUnion;
96
97 /* Add two 'long double' values and return the result. */
98 long double
99 __gcc_qadd (double a, double aa, double c, double cc)
100 {
101 longDblUnion x;
102 double z, q, zz, xh;
103
104 z = a + c;
105
106 if (nonfinite (z))
107 {
108 z = cc + aa + c + a;
109 if (nonfinite (z))
110 return z;
111 x.dval[0] = z; /* Will always be DBL_MAX. */
112 zz = aa + cc;
113 if (fabs(a) > fabs(c))
114 x.dval[1] = a - z + c + zz;
115 else
116 x.dval[1] = c - z + a + zz;
117 }
118 else
119 {
120 q = a - z;
121 zz = q + c + (a - (q + z)) + aa + cc;
122
123 /* Keep -0 result. */
124 if (zz == 0.0)
125 return z;
126
127 xh = z + zz;
128 if (nonfinite (xh))
129 return xh;
130
131 x.dval[0] = xh;
132 x.dval[1] = z - xh + zz;
133 }
134 return x.ldval;
135 }
136
137 long double
138 __gcc_qsub (double a, double b, double c, double d)
139 {
140 return __gcc_qadd (a, b, -c, -d);
141 }
142
143 #ifdef __NO_FPRS__
144 static double fmsub (double, double, double);
145 #endif
146
147 long double
148 __gcc_qmul (double a, double b, double c, double d)
149 {
150 longDblUnion z;
151 double t, tau, u, v, w;
152
153 t = a * c; /* Highest order double term. */
154
155 if (unlikely (t == 0) /* Preserve -0. */
156 || nonfinite (t))
157 return t;
158
159 /* Sum terms of two highest orders. */
160
161 /* Use fused multiply-add to get low part of a * c. */
162 #ifndef __NO_FPRS__
163 asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
164 #else
165 tau = fmsub (a, c, t);
166 #endif
167 v = a*d;
168 w = b*c;
169 tau += v + w; /* Add in other second-order terms. */
170 u = t + tau;
171
172 /* Construct long double result. */
173 if (nonfinite (u))
174 return u;
175 z.dval[0] = u;
176 z.dval[1] = (t - u) + tau;
177 return z.ldval;
178 }
179
180 long double
181 __gcc_qdiv (double a, double b, double c, double d)
182 {
183 longDblUnion z;
184 double s, sigma, t, tau, u, v, w;
185
186 t = a / c; /* highest order double term */
187
188 if (unlikely (t == 0) /* Preserve -0. */
189 || nonfinite (t))
190 return t;
191
192 /* Finite nonzero result requires corrections to the highest order term. */
193
194 s = c * t; /* (s,sigma) = c*t exactly. */
195 w = -(-b + d * t); /* Written to get fnmsub for speed, but not
196 numerically necessary. */
197
198 /* Use fused multiply-add to get low part of c * t. */
199 #ifndef __NO_FPRS__
200 asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
201 #else
202 sigma = fmsub (c, t, s);
203 #endif
204 v = a - s;
205
206 tau = ((v-sigma)+w)/c; /* Correction to t. */
207 u = t + tau;
208
209 /* Construct long double result. */
210 if (nonfinite (u))
211 return u;
212 z.dval[0] = u;
213 z.dval[1] = (t - u) + tau;
214 return z.ldval;
215 }
216
217 #if defined (_SOFT_DOUBLE) && defined (__LONG_DOUBLE_128__)
218
219 long double __gcc_qneg (double, double);
220 int __gcc_qeq (double, double, double, double);
221 int __gcc_qne (double, double, double, double);
222 int __gcc_qge (double, double, double, double);
223 int __gcc_qle (double, double, double, double);
224 long double __gcc_stoq (float);
225 long double __gcc_dtoq (double);
226 float __gcc_qtos (double, double);
227 double __gcc_qtod (double, double);
228 int __gcc_qtoi (double, double);
229 unsigned int __gcc_qtou (double, double);
230 long double __gcc_itoq (int);
231 long double __gcc_utoq (unsigned int);
232
233 extern int __eqdf2 (double, double);
234 extern int __ledf2 (double, double);
235 extern int __gedf2 (double, double);
236
237 /* Negate 'long double' value and return the result. */
238 long double
239 __gcc_qneg (double a, double aa)
240 {
241 longDblUnion x;
242
243 x.dval[0] = -a;
244 x.dval[1] = -aa;
245 return x.ldval;
246 }
247
248 /* Compare two 'long double' values for equality. */
249 int
250 __gcc_qeq (double a, double aa, double c, double cc)
251 {
252 if (__eqdf2 (a, c) == 0)
253 return __eqdf2 (aa, cc);
254 return 1;
255 }
256
257 strong_alias (__gcc_qeq, __gcc_qne);
258
259 /* Compare two 'long double' values for less than or equal. */
260 int
261 __gcc_qle (double a, double aa, double c, double cc)
262 {
263 if (__eqdf2 (a, c) == 0)
264 return __ledf2 (aa, cc);
265 return __ledf2 (a, c);
266 }
267
268 strong_alias (__gcc_qle, __gcc_qlt);
269
270 /* Compare two 'long double' values for greater than or equal. */
271 int
272 __gcc_qge (double a, double aa, double c, double cc)
273 {
274 if (__eqdf2 (a, c) == 0)
275 return __gedf2 (aa, cc);
276 return __gedf2 (a, c);
277 }
278
279 strong_alias (__gcc_qge, __gcc_qgt);
280
281 /* Convert single to long double. */
282 long double
283 __gcc_stoq (float a)
284 {
285 longDblUnion x;
286
287 x.dval[0] = (double) a;
288 x.dval[1] = 0.0;
289
290 return x.ldval;
291 }
292
293 /* Convert double to long double. */
294 long double
295 __gcc_dtoq (double a)
296 {
297 longDblUnion x;
298
299 x.dval[0] = a;
300 x.dval[1] = 0.0;
301
302 return x.ldval;
303 }
304
305 /* Convert long double to single. */
306 float
307 __gcc_qtos (double a, double aa __attribute__ ((__unused__)))
308 {
309 return (float) a;
310 }
311
312 /* Convert long double to double. */
313 double
314 __gcc_qtod (double a, double aa __attribute__ ((__unused__)))
315 {
316 return a;
317 }
318
319 /* Convert long double to int. */
320 int
321 __gcc_qtoi (double a, double aa)
322 {
323 double z = a + aa;
324 return (int) z;
325 }
326
327 /* Convert long double to unsigned int. */
328 unsigned int
329 __gcc_qtou (double a, double aa)
330 {
331 double z = a + aa;
332 return (unsigned int) z;
333 }
334
335 /* Convert int to long double. */
336 long double
337 __gcc_itoq (int a)
338 {
339 return __gcc_dtoq ((double) a);
340 }
341
342 /* Convert unsigned int to long double. */
343 long double
344 __gcc_utoq (unsigned int a)
345 {
346 return __gcc_dtoq ((double) a);
347 }
348
349 #endif
350
351 #ifdef __NO_FPRS__
352
353 int __gcc_qunord (double, double, double, double);
354
355 extern int __eqdf2 (double, double);
356 extern int __unorddf2 (double, double);
357
358 /* Compare two 'long double' values for unordered. */
359 int
360 __gcc_qunord (double a, double aa, double c, double cc)
361 {
362 if (__eqdf2 (a, c) == 0)
363 return __unorddf2 (aa, cc);
364 return __unorddf2 (a, c);
365 }
366
367 #include "config/soft-fp/soft-fp.h"
368 #include "config/soft-fp/double.h"
369 #include "config/soft-fp/quad.h"
370
371 /* Compute floating point multiply-subtract with higher (quad) precision. */
372 static double
373 fmsub (double a, double b, double c)
374 {
375 FP_DECL_EX;
376 FP_DECL_D(A);
377 FP_DECL_D(B);
378 FP_DECL_D(C);
379 FP_DECL_Q(X);
380 FP_DECL_Q(Y);
381 FP_DECL_Q(Z);
382 FP_DECL_Q(U);
383 FP_DECL_Q(V);
384 FP_DECL_D(R);
385 double r;
386 long double u, x, y, z;
387
388 FP_INIT_ROUNDMODE;
389 FP_UNPACK_RAW_D (A, a);
390 FP_UNPACK_RAW_D (B, b);
391 FP_UNPACK_RAW_D (C, c);
392
393 /* Extend double to quad. */
394 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
395 FP_EXTEND(Q,D,4,2,X,A);
396 FP_EXTEND(Q,D,4,2,Y,B);
397 FP_EXTEND(Q,D,4,2,Z,C);
398 #else
399 FP_EXTEND(Q,D,2,1,X,A);
400 FP_EXTEND(Q,D,2,1,Y,B);
401 FP_EXTEND(Q,D,2,1,Z,C);
402 #endif
403 FP_PACK_RAW_Q(x,X);
404 FP_PACK_RAW_Q(y,Y);
405 FP_PACK_RAW_Q(z,Z);
406 FP_HANDLE_EXCEPTIONS;
407
408 /* Multiply. */
409 FP_INIT_ROUNDMODE;
410 FP_UNPACK_Q(X,x);
411 FP_UNPACK_Q(Y,y);
412 FP_MUL_Q(U,X,Y);
413 FP_PACK_Q(u,U);
414 FP_HANDLE_EXCEPTIONS;
415
416 /* Subtract. */
417 FP_INIT_ROUNDMODE;
418 FP_UNPACK_SEMIRAW_Q(U,u);
419 FP_UNPACK_SEMIRAW_Q(Z,z);
420 FP_SUB_Q(V,U,Z);
421
422 /* Truncate quad to double. */
423 #if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
424 V_f[3] &= 0x0007ffff;
425 FP_TRUNC(D,Q,2,4,R,V);
426 #else
427 V_f1 &= 0x0007ffffffffffffL;
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