1/* 128-bit long double support routines for Darwin. 2 Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007 3 Free Software Foundation, Inc. 4 5This file is part of GCC. 6 7GCC is free software; you can redistribute it and/or modify it under 8the terms of the GNU General Public License as published by the Free 9Software Foundation; either version 2, or (at your option) any later 10version. 11 12In addition to the permissions in the GNU General Public License, the 13Free Software Foundation gives you unlimited permission to link the 14compiled version of this file into combinations with other programs, 15and to distribute those combinations without any restriction coming 16from the use of this file. (The General Public License restrictions 17do apply in other respects; for example, they cover modification of 18the file, and distribution when not linked into a combine 19executable.) 20 21GCC is distributed in the hope that it will be useful, but WITHOUT ANY 22WARRANTY; without even the implied warranty of MERCHANTABILITY or 23FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 24for more details. 25 26You should have received a copy of the GNU General Public License 27along with GCC; see the file COPYING. If not, write to the Free 28Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 2902110-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 73long double __gcc_qadd (double, double, double, double); 74long double __gcc_qsub (double, double, double, double); 75long double __gcc_qmul (double, double, double, double); 76long 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 96typedef union 97{ 98 long double ldval; 99 double dval[2]; 100} longDblUnion; 101 102/* Add two 'long double' values and return the result. */ 103long 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 142long 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 149static double fmsub (double, double, double); 150#endif 151 152long 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 185long 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#if defined (_SOFT_FLOAT) && defined (__LONG_DOUBLE_128__) 223 224long double __gcc_qneg (double, double); 225int __gcc_qeq (double, double, double, double); 226int __gcc_qne (double, double, double, double); 227int __gcc_qge (double, double, double, double); 228int __gcc_qle (double, double, double, double); 229int __gcc_qunord (double, double, double, double); 230long double __gcc_stoq (float); 231long double __gcc_dtoq (double); 232float __gcc_qtos (double, double); 233double __gcc_qtod (double, double); 234int __gcc_qtoi (double, double); 235unsigned int __gcc_qtou (double, double); 236long double __gcc_itoq (int); 237long double __gcc_utoq (unsigned int); 238 239extern int __eqdf2 (double, double); 240extern int __ledf2 (double, double); 241extern int __gedf2 (double, double); 242extern int __unorddf2 (double, double); 243 244/* Negate 'long double' value and return the result. */ 245long 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. */ 256int 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 264strong_alias (__gcc_qeq, __gcc_qne); 265 266/* Compare two 'long double' values for less than or equal. */ 267int 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 275strong_alias (__gcc_qle, __gcc_qlt); 276 277/* Compare two 'long double' values for greater than or equal. */ 278int 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 286strong_alias (__gcc_qge, __gcc_qgt); 287 288/* Compare two 'long double' values for unordered. */ 289int 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. */ 298long 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. */ 310long 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. */ 322float 323__gcc_qtos (double a, double aa __attribute__ ((__unused__))) 324{ 325 return (float) a; 326} 327 328/* Convert long double to double. */ 329double 330__gcc_qtod (double a, double aa __attribute__ ((__unused__))) 331{ 332 return a; 333} 334 335/* Convert long double to int. */ 336int 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. */ 344unsigned 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. */ 352long double 353__gcc_itoq (int a) 354{ 355 return __gcc_dtoq ((double) a); 356} 357 358/* Convert unsigned int to long double. */ 359long 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. */ 370static double 371fmsub (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 439