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