|
0
|
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
|
