1 /* 128-bit long double support routines for Darwin.
2 Copyright (C) 1993-2023 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 && !defined (__rtems__) \
50 && (defined (__LONG_DOUBLE_128__) || defined (__FLOAT128_TYPE__))
51
52 #define fabs(x) __builtin_fabs(x)
53 #define isless(x, y) __builtin_isless (x, y)
54 #define inf() __builtin_inf()
55
56 #define unlikely(x) __builtin_expect ((x), 0)
57
58 #define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
59
60 /* If we have __float128/_Float128, use __ibm128 instead of long double. On
61 other systems, use long double, because __ibm128 might not have been
62 created. */
63 #ifdef __FLOAT128__
64 #define IBM128_TYPE __ibm128
65 #else
66 #define IBM128_TYPE long double
67 #endif
68
69 /* Define ALIASNAME as a strong alias for NAME. */
70 # define strong_alias(name, aliasname) _strong_alias(name, aliasname)
71 # define _strong_alias(name, aliasname) \
72 extern __typeof (name) aliasname __attribute__ ((alias (#name)));
73
74 /* All these routines actually take two long doubles as parameters,
75 but GCC currently generates poor code when a union is used to turn
76 a long double into a pair of doubles. */
77
78 IBM128_TYPE __gcc_qadd (double, double, double, double);
79 IBM128_TYPE __gcc_qsub (double, double, double, double);
80 IBM128_TYPE __gcc_qmul (double, double, double, double);
81 IBM128_TYPE __gcc_qdiv (double, double, double, double);
82
83 #if defined __ELF__ && defined SHARED \
84 && (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))
85 /* Provide definitions of the old symbol names to satisfy apps and
86 shared libs built against an older libgcc. To access the _xlq
87 symbols an explicit version reference is needed, so these won't
88 satisfy an unadorned reference like _xlqadd. If dot symbols are
89 not needed, the assembler will remove the aliases from the symbol
90 table. */
91 __asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
92 ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
93 ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
94 ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
95 ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
96 ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
97 ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
98 ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
99 #endif
100
101 /* Combine two 'double' values into one 'IBM128_TYPE' and return the result. */
102 static inline IBM128_TYPE
103 pack_ldouble (double dh, double dl)
104 {
105 #if defined (__LONG_DOUBLE_128__) && defined (__LONG_DOUBLE_IBM128__) \
106 && !(defined (_SOFT_FLOAT) || defined (__NO_FPRS__))
107 return __builtin_pack_longdouble (dh, dl);
108 #else
109 union
110 {
111 IBM128_TYPE ldval;
112 double dval[2];
113 } x;
114 x.dval[0] = dh;
115 x.dval[1] = dl;
116 return x.ldval;
117 #endif
118 }
119
120 /* Add two 'IBM128_TYPE' values and return the result. */
121 static inline IBM128_TYPE
122 ldouble_qadd_internal (double a, double aa, double c, double cc)
123 {
124 double xh, xl, z, q, zz;
125
126 z = a + c;
127
128 if (nonfinite (z))
129 {
130 if (fabs (z) != inf())
131 return z;
132 z = cc + aa + c + a;
133 if (nonfinite (z))
134 return z;
135 xh = z; /* Will always be DBL_MAX. */
136 zz = aa + cc;
137 if (fabs(a) > fabs(c))
138 xl = a - z + c + zz;
139 else
140 xl = c - z + a + zz;
141 }
142 else
143 {
144 q = a - z;
145 zz = q + c + (a - (q + z)) + aa + cc;
146
147 /* Keep -0 result. */
148 if (zz == 0.0)
149 return z;
150
151 xh = z + zz;
152 if (nonfinite (xh))
153 return xh;
154
155 xl = z - xh + zz;
156 }
157 return pack_ldouble (xh, xl);
158 }
159
160 IBM128_TYPE
161 __gcc_qadd (double a, double aa, double c, double cc)
162 {
163 return ldouble_qadd_internal (a, aa, c, cc);
164 }
165
166 IBM128_TYPE
167 __gcc_qsub (double a, double aa, double c, double cc)
168 {
169 return ldouble_qadd_internal (a, aa, -c, -cc);
170 }
171
172 #ifdef __NO_FPRS__
173 static double fmsub (double, double, double);
174 #endif
175
176 IBM128_TYPE
177 __gcc_qmul (double a, double b, double c, double d)
178 {
179 double xh, xl, t, tau, u, v, w;
180
181 t = a * c; /* Highest order double term. */
182
183 if (unlikely (t == 0) /* Preserve -0. */
184 || nonfinite (t))
185 return t;
186
187 /* Sum terms of two highest orders. */
188
189 /* Use fused multiply-add to get low part of a * c. */
190 #ifndef __NO_FPRS__
191 asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
192 #else
193 tau = fmsub (a, c, t);
194 #endif
195 v = a*d;
196 w = b*c;
197 tau += v + w; /* Add in other second-order terms. */
198 u = t + tau;
199
200 /* Construct IBM128_TYPE result. */
201 if (nonfinite (u))
202 return u;
203 xh = u;
204 xl = (t - u) + tau;
205 return pack_ldouble (xh, xl);
206 }
207
208 IBM128_TYPE
209 __gcc_qdiv (double a, double b, double c, double d)
210 {
211 double xh, xl, s, sigma, t, tau, u, v, w;
212
213 t = a / c; /* highest order double term */
214
215 if (unlikely (t == 0) /* Preserve -0. */
216 || nonfinite (t))
217 return t;
218
219 /* Finite nonzero result requires corrections to the highest order
220 term. These corrections require the low part of c * t to be
221 exactly represented in double. */
222 if (fabs (a) <= 0x1p-969)
223 {
224 a *= 0x1p106;
225 b *= 0x1p106;
226 c *= 0x1p106;
227 d *= 0x1p106;
228 }
229
230 s = c * t; /* (s,sigma) = c*t exactly. */
231 w = -(-b + d * t); /* Written to get fnmsub for speed, but not
232 numerically necessary. */
233
234 /* Use fused multiply-add to get low part of c * t. */
235 #ifndef __NO_FPRS__
236 asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
237 #else
238 sigma = fmsub (c, t, s);
239 #endif
240 v = a - s;
241
242 tau = ((v-sigma)+w)/c; /* Correction to t. */
243 u = t + tau;
244
245 /* Construct IBM128_TYPE result. */
246 if (nonfinite (u))
247 return u;
248 xh = u;
249 xl = (t - u) + tau;
250 return pack_ldouble (xh, xl);
251 }
252
253 #if defined (_SOFT_DOUBLE) && defined (__LONG_DOUBLE_128__)
254
255 IBM128_TYPE __gcc_qneg (double, double);
256 int __gcc_qeq (double, double, double, double);
257 int __gcc_qne (double, double, double, double);
258 int __gcc_qge (double, double, double, double);
259 int __gcc_qle (double, double, double, double);
260 IBM128_TYPE __gcc_stoq (float);
261 IBM128_TYPE __gcc_dtoq (double);
262 float __gcc_qtos (double, double);
263 double __gcc_qtod (double, double);
264 int __gcc_qtoi (double, double);
265 unsigned int __gcc_qtou (double, double);
266 IBM128_TYPE __gcc_itoq (int);
267 IBM128_TYPE __gcc_utoq (unsigned int);
268
269 extern int __eqdf2 (double, double);
270 extern int __ledf2 (double, double);
271 extern int __gedf2 (double, double);
272
273 /* Negate 'IBM128_TYPE' value and return the result. */
274 IBM128_TYPE
275 __gcc_qneg (double a, double aa)
276 {
277 return pack_ldouble (-a, -aa);
278 }
279
280 /* Compare two 'IBM128_TYPE' values for equality. */
281 int
282 __gcc_qeq (double a, double aa, double c, double cc)
283 {
284 if (__eqdf2 (a, c) == 0)
285 return __eqdf2 (aa, cc);
286 return 1;
287 }
288
289 strong_alias (__gcc_qeq, __gcc_qne);
290
291 /* Compare two 'IBM128_TYPE' values for less than or equal. */
292 int
293 __gcc_qle (double a, double aa, double c, double cc)
294 {
295 if (__eqdf2 (a, c) == 0)
296 return __ledf2 (aa, cc);
297 return __ledf2 (a, c);
298 }
299
300 strong_alias (__gcc_qle, __gcc_qlt);
301
302 /* Compare two 'IBM128_TYPE' values for greater than or equal. */
303 int
304 __gcc_qge (double a, double aa, double c, double cc)
305 {
306 if (__eqdf2 (a, c) == 0)
307 return __gedf2 (aa, cc);
308 return __gedf2 (a, c);
309 }
310
311 strong_alias (__gcc_qge, __gcc_qgt);
312
313 /* Convert single to IBM128_TYPE. */
314 IBM128_TYPE
315 __gcc_stoq (float a)
316 {
317 return pack_ldouble ((double) a, 0.0);
318 }
319
320 /* Convert double to IBM128_TYPE. */
321 IBM128_TYPE
322 __gcc_dtoq (double a)
323 {
324 return pack_ldouble (a, 0.0);
325 }
326
327 /* Convert IBM128_TYPE to single. */
328 float
329 __gcc_qtos (double a, double aa __attribute__ ((__unused__)))
330 {
331 return (float) a;
332 }
333
334 /* Convert IBM128_TYPE to double. */
335 double
336 __gcc_qtod (double a, double aa __attribute__ ((__unused__)))
337 {
338 return a;
339 }
340
341 /* Convert IBM128_TYPE to int. */
342 int
343 __gcc_qtoi (double a, double aa)
344 {
345 double z = a + aa;
346 return (int) z;
347 }
348
349 /* Convert IBM128_TYPE to unsigned int. */
350 unsigned int
351 __gcc_qtou (double a, double aa)
352 {
353 double z = a + aa;
354 return (unsigned int) z;
355 }
356
357 /* Convert int to IBM128_TYPE. */
358 IBM128_TYPE
359 __gcc_itoq (int a)
360 {
361 return __gcc_dtoq ((double) a);
362 }
363
364 /* Convert unsigned int to IBM128_TYPE. */
365 IBM128_TYPE
366 __gcc_utoq (unsigned int a)
367 {
368 return __gcc_dtoq ((double) a);
369 }
370
371 #endif
372
373 #ifdef __NO_FPRS__
374
375 int __gcc_qunord (double, double, double, double);
376
377 extern int __eqdf2 (double, double);
378 extern int __unorddf2 (double, double);
379
380 /* Compare two 'IBM128_TYPE' values for unordered. */
381 int
382 __gcc_qunord (double a, double aa, double c, double cc)
383 {
384 if (__eqdf2 (a, c) == 0)
385 return __unorddf2 (aa, cc);
386 return __unorddf2 (a, c);
387 }
388
389 #include "soft-fp/soft-fp.h"
390 #include "soft-fp/double.h"
391 #include "soft-fp/quad.h"
392
393 /* Compute floating point multiply-subtract with higher (quad) precision. */
394 static double
395 fmsub (double a, double b, double c)
396 {
397 FP_DECL_EX;
398 FP_DECL_D(A);
399 FP_DECL_D(B);
400 FP_DECL_D(C);
401 FP_DECL_Q(X);
402 FP_DECL_Q(Y);
403 FP_DECL_Q(Z);
404 FP_DECL_Q(U);
405 FP_DECL_Q(V);
406 FP_DECL_D(R);
407 double r;
408 IBM128_TYPE u, x, y, z;
409
410 FP_INIT_ROUNDMODE;
411 FP_UNPACK_RAW_D (A, a);
412 FP_UNPACK_RAW_D (B, b);
413 FP_UNPACK_RAW_D (C, c);
414
415 /* Extend double to quad. */
416 #if _FP_W_TYPE_SIZE < 64
417 FP_EXTEND(Q,D,4,2,X,A);
418 FP_EXTEND(Q,D,4,2,Y,B);
419 FP_EXTEND(Q,D,4,2,Z,C);
420 #else
421 FP_EXTEND(Q,D,2,1,X,A);
422 FP_EXTEND(Q,D,2,1,Y,B);
423 FP_EXTEND(Q,D,2,1,Z,C);
424 #endif
425 FP_PACK_RAW_Q(x,X);
426 FP_PACK_RAW_Q(y,Y);
427 FP_PACK_RAW_Q(z,Z);
428 FP_HANDLE_EXCEPTIONS;
429
430 /* Multiply. */
431 FP_INIT_ROUNDMODE;
432 FP_UNPACK_Q(X,x);
433 FP_UNPACK_Q(Y,y);
434 FP_MUL_Q(U,X,Y);
435 FP_PACK_Q(u,U);
436 FP_HANDLE_EXCEPTIONS;
437
438 /* Subtract. */
439 FP_INIT_ROUNDMODE;
440 FP_UNPACK_SEMIRAW_Q(U,u);
441 FP_UNPACK_SEMIRAW_Q(Z,z);
442 FP_SUB_Q(V,U,Z);
443
444 /* Truncate quad to double. */
445 #if _FP_W_TYPE_SIZE < 64
446 V_f[3] &= 0x0007ffff;
447 FP_TRUNC(D,Q,2,4,R,V);
448 #else
449 V_f1 &= 0x0007ffffffffffffL;
450 FP_TRUNC(D,Q,1,2,R,V);
451 #endif
452 FP_PACK_SEMIRAW_D(r,R);
453 FP_HANDLE_EXCEPTIONS;
454
455 return r;
456 }
457
458 #endif
459
460 #endif