darwin-ldouble.c revision 146895
1132718Skan/* 128-bit long double support routines for Darwin. 2146895Skan Copyright (C) 1993, 2003, 2004, 2005 Free Software Foundation, Inc. 3132718Skan 4132718SkanThis file is part of GCC. 5132718Skan 6132718SkanGCC is free software; you can redistribute it and/or modify it under 7132718Skanthe terms of the GNU General Public License as published by the Free 8132718SkanSoftware Foundation; either version 2, or (at your option) any later 9132718Skanversion. 10132718Skan 11132718SkanIn addition to the permissions in the GNU General Public License, the 12132718SkanFree Software Foundation gives you unlimited permission to link the 13132718Skancompiled version of this file into combinations with other programs, 14132718Skanand to distribute those combinations without any restriction coming 15132718Skanfrom the use of this file. (The General Public License restrictions 16132718Skando apply in other respects; for example, they cover modification of 17132718Skanthe file, and distribution when not linked into a combine 18132718Skanexecutable.) 19132718Skan 20132718SkanGCC is distributed in the hope that it will be useful, but WITHOUT ANY 21132718SkanWARRANTY; without even the implied warranty of MERCHANTABILITY or 22132718SkanFITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 23132718Skanfor more details. 24132718Skan 25132718SkanYou should have received a copy of the GNU General Public License 26132718Skanalong with GCC; see the file COPYING. If not, write to the Free 27132718SkanSoftware Foundation, 59 Temple Place - Suite 330, Boston, MA 28132718Skan02111-1307, USA. */ 29132718Skan 30132718Skan/* Implementations of floating-point long double basic arithmetic 31132718Skan functions called by the IBM C compiler when generating code for 32132718Skan PowerPC platforms. In particular, the following functions are 33132718Skan implemented: _xlqadd, _xlqsub, _xlqmul, and _xlqdiv. Double-double 34132718Skan algorithms are based on the paper "Doubled-Precision IEEE Standard 35132718Skan 754 Floating-Point Arithmetic" by W. Kahan, February 26, 1987. An 36132718Skan alternative published reference is "Software for Doubled-Precision 37132718Skan Floating-Point Computations", by Seppo Linnainmaa, ACM TOMS vol 7 38132718Skan no 3, September 1961, pages 272-283. */ 39132718Skan 40132718Skan/* Each long double is made up of two IEEE doubles. The value of the 41132718Skan long double is the sum of the values of the two parts. The most 42132718Skan significant part is required to be the value of the long double 43132718Skan rounded to the nearest double, as specified by IEEE. For Inf 44132718Skan values, the least significant part is required to be one of +0.0 or 45132718Skan -0.0. No other requirements are made; so, for example, 1.0 may be 46132718Skan represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a 47132718Skan NaN is don't-care. 48132718Skan 49132718Skan This code currently assumes big-endian. */ 50132718Skan 51146895Skan#if !_SOFT_FLOAT && (defined (__MACH__) || defined (__powerpc64__) || defined (_AIX)) 52132718Skan 53132718Skan#define fabs(x) __builtin_fabs(x) 54132718Skan 55132718Skan#define unlikely(x) __builtin_expect ((x), 0) 56132718Skan 57132718Skan/* All these routines actually take two long doubles as parameters, 58132718Skan but GCC currently generates poor code when a union is used to turn 59132718Skan a long double into a pair of doubles. */ 60132718Skan 61146895Skanextern long double __gcc_qadd (double, double, double, double); 62146895Skanextern long double __gcc_qsub (double, double, double, double); 63146895Skanextern long double __gcc_qmul (double, double, double, double); 64146895Skanextern long double __gcc_qdiv (double, double, double, double); 65132718Skan 66146895Skan#if defined __ELF__ && defined IN_LIBGCC2_S 67146895Skan/* Provide definitions of the old symbol names to statisfy apps and 68146895Skan shared libs built against an older libgcc. To access the _xlq 69146895Skan symbols an explicit version reference is needed, so these won't 70146895Skan satisfy an unadorned reference like _xlqadd. If dot symbols are 71146895Skan not needed, the assembler will remove the aliases from the symbol 72146895Skan table. */ 73146895Skan__asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t" 74146895Skan ".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t" 75146895Skan ".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t" 76146895Skan ".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t" 77146895Skan ".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t" 78146895Skan ".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t" 79146895Skan ".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t" 80146895Skan ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4"); 81146895Skan#endif 82146895Skan 83132718Skantypedef union 84132718Skan{ 85132718Skan long double ldval; 86132718Skan double dval[2]; 87132718Skan} longDblUnion; 88132718Skan 89132718Skanstatic const double FPKINF = 1.0/0.0; 90132718Skan 91132718Skan/* Add two 'long double' values and return the result. */ 92132718Skanlong double 93146895Skan__gcc_qadd (double a, double b, double c, double d) 94132718Skan{ 95132718Skan longDblUnion z; 96132718Skan double t, tau, u, FPR_zero, FPR_PosInf; 97132718Skan 98132718Skan FPR_zero = 0.0; 99132718Skan FPR_PosInf = FPKINF; 100132718Skan 101132718Skan if (unlikely (a != a) || unlikely (c != c)) 102132718Skan return a + c; /* NaN result. */ 103132718Skan 104132718Skan /* Ordered operands are arranged in order of their magnitudes. */ 105132718Skan 106132718Skan /* Switch inputs if |(c,d)| > |(a,b)|. */ 107132718Skan if (fabs (c) > fabs (a)) 108132718Skan { 109132718Skan t = a; 110132718Skan tau = b; 111132718Skan a = c; 112132718Skan b = d; 113132718Skan c = t; 114132718Skan d = tau; 115132718Skan } 116132718Skan 117132718Skan /* b <- second largest magnitude double. */ 118132718Skan if (fabs (c) > fabs (b)) 119132718Skan { 120132718Skan t = b; 121132718Skan b = c; 122132718Skan c = t; 123132718Skan } 124132718Skan 125132718Skan /* Thanks to commutivity, sum is invariant w.r.t. the next 126132718Skan conditional exchange. */ 127132718Skan tau = d + c; 128132718Skan 129132718Skan /* Order the smallest magnitude doubles. */ 130132718Skan if (fabs (d) > fabs (c)) 131132718Skan { 132132718Skan t = c; 133132718Skan c = d; 134132718Skan d = t; 135132718Skan } 136132718Skan 137132718Skan t = (tau + b) + a; /* Sum values in ascending magnitude order. */ 138132718Skan 139132718Skan /* Infinite or zero result. */ 140132718Skan if (unlikely (t == FPR_zero) || unlikely (fabs (t) == FPR_PosInf)) 141132718Skan return t; 142132718Skan 143132718Skan /* Usual case. */ 144132718Skan tau = (((a-t) + b) + c) + d; 145132718Skan u = t + tau; 146132718Skan z.dval[0] = u; /* Final fixup for long double result. */ 147132718Skan z.dval[1] = (t - u) + tau; 148132718Skan return z.ldval; 149132718Skan} 150132718Skan 151132718Skanlong double 152146895Skan__gcc_qsub (double a, double b, double c, double d) 153132718Skan{ 154146895Skan return __gcc_qadd (a, b, -c, -d); 155132718Skan} 156132718Skan 157132718Skanlong double 158146895Skan__gcc_qmul (double a, double b, double c, double d) 159132718Skan{ 160132718Skan longDblUnion z; 161132718Skan double t, tau, u, v, w, FPR_zero, FPR_PosInf; 162132718Skan 163132718Skan FPR_zero = 0.0; 164132718Skan FPR_PosInf = FPKINF; 165132718Skan 166132718Skan t = a * c; /* Highest order double term. */ 167132718Skan 168132718Skan if (unlikely (t != t) || unlikely (t == FPR_zero) 169132718Skan || unlikely (fabs (t) == FPR_PosInf)) 170132718Skan return t; 171132718Skan 172132718Skan /* Finite nonzero result requires summing of terms of two highest 173132718Skan orders. */ 174132718Skan 175132718Skan /* Use fused multiply-add to get low part of a * c. */ 176132718Skan asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t)); 177132718Skan v = a*d; 178132718Skan w = b*c; 179132718Skan tau += v + w; /* Add in other second-order terms. */ 180132718Skan u = t + tau; 181132718Skan 182132718Skan /* Construct long double result. */ 183132718Skan z.dval[0] = u; 184132718Skan z.dval[1] = (t - u) + tau; 185132718Skan return z.ldval; 186132718Skan} 187132718Skan 188132718Skanlong double 189146895Skan__gcc_qdiv (double a, double b, double c, double d) 190132718Skan{ 191132718Skan longDblUnion z; 192132718Skan double s, sigma, t, tau, u, v, w, FPR_zero, FPR_PosInf; 193132718Skan 194132718Skan FPR_zero = 0.0; 195132718Skan FPR_PosInf = FPKINF; 196132718Skan 197132718Skan t = a / c; /* highest order double term */ 198132718Skan 199132718Skan if (unlikely (t != t) || unlikely (t == FPR_zero) 200132718Skan || unlikely (fabs (t) == FPR_PosInf)) 201132718Skan return t; 202132718Skan 203132718Skan /* Finite nonzero result requires corrections to the highest order term. */ 204132718Skan 205132718Skan s = c * t; /* (s,sigma) = c*t exactly. */ 206132718Skan w = -(-b + d * t); /* Written to get fnmsub for speed, but not 207132718Skan numerically necessary. */ 208132718Skan 209132718Skan /* Use fused multiply-add to get low part of c * t. */ 210132718Skan asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s)); 211132718Skan v = a - s; 212132718Skan 213132718Skan tau = ((v-sigma)+w)/c; /* Correction to t. */ 214132718Skan u = t + tau; 215132718Skan 216132718Skan /* Construct long double result. */ 217132718Skan z.dval[0] = u; 218132718Skan z.dval[1] = (t - u) + tau; 219132718Skan return z.ldval; 220132718Skan} 221132718Skan 222132718Skan#endif 223