darwin-ldouble.c revision 132718
175295Sdes/* 128-bit long double support routines for Darwin. 275295Sdes Copyright (C) 1993, 2003, 2004 Free Software Foundation, Inc. 375295Sdes 475295SdesThis file is part of GCC. 575295Sdes 675295SdesGCC is free software; you can redistribute it and/or modify it under 775295Sdesthe terms of the GNU General Public License as published by the Free 875295SdesSoftware Foundation; either version 2, or (at your option) any later 975295Sdesversion. 1075295Sdes 1175295SdesIn addition to the permissions in the GNU General Public License, the 1275295SdesFree Software Foundation gives you unlimited permission to link the 1375295Sdescompiled version of this file into combinations with other programs, 1475295Sdesand to distribute those combinations without any restriction coming 1575295Sdesfrom the use of this file. (The General Public License restrictions 1675295Sdesdo apply in other respects; for example, they cover modification of 1775295Sdesthe file, and distribution when not linked into a combine 1875295Sdesexecutable.) 1975295Sdes 2075295SdesGCC is distributed in the hope that it will be useful, but WITHOUT ANY 2175295SdesWARRANTY; without even the implied warranty of MERCHANTABILITY or 2275295SdesFITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 2375295Sdesfor more details. 2475295Sdes 2575295SdesYou should have received a copy of the GNU General Public License 2675295Sdesalong with GCC; see the file COPYING. If not, write to the Free 2775295SdesSoftware Foundation, 59 Temple Place - Suite 330, Boston, MA 2875295Sdes02111-1307, USA. */ 2975295Sdes 3075295Sdes/* Implementations of floating-point long double basic arithmetic 3175295Sdes functions called by the IBM C compiler when generating code for 3275295Sdes PowerPC platforms. In particular, the following functions are 3375295Sdes implemented: _xlqadd, _xlqsub, _xlqmul, and _xlqdiv. Double-double 3478073Sdes algorithms are based on the paper "Doubled-Precision IEEE Standard 3575295Sdes 754 Floating-Point Arithmetic" by W. Kahan, February 26, 1987. An 3677965Sdes alternative published reference is "Software for Doubled-Precision 3784246Sdes Floating-Point Computations", by Seppo Linnainmaa, ACM TOMS vol 7 3875295Sdes no 3, September 1961, pages 272-283. */ 3975295Sdes 4075295Sdes/* Each long double is made up of two IEEE doubles. The value of the 4175295Sdes long double is the sum of the values of the two parts. The most 4275295Sdes significant part is required to be the value of the long double 4375295Sdes rounded to the nearest double, as specified by IEEE. For Inf 4477998Sdes values, the least significant part is required to be one of +0.0 or 4575295Sdes -0.0. No other requirements are made; so, for example, 1.0 may be 4675295Sdes represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a 4784246Sdes NaN is don't-care. 4884246Sdes 4975295Sdes This code currently assumes big-endian. */ 5075295Sdes 5175295Sdes#if !_SOFT_FLOAT && (defined (__MACH__) || defined (__powerpc64__)) 5275295Sdes 5384246Sdes#define fabs(x) __builtin_fabs(x) 5484246Sdes 5584246Sdes#define unlikely(x) __builtin_expect ((x), 0) 5684246Sdes 5784246Sdes/* All these routines actually take two long doubles as parameters, 5884246Sdes but GCC currently generates poor code when a union is used to turn 5984246Sdes a long double into a pair of doubles. */ 6084246Sdes 6184246Sdesextern long double _xlqadd (double, double, double, double); 6284246Sdesextern long double _xlqsub (double, double, double, double); 6375295Sdesextern long double _xlqmul (double, double, double, double); 6484246Sdesextern long double _xlqdiv (double, double, double, double); 6584246Sdes 6675295Sdestypedef union 6775295Sdes{ 6875295Sdes long double ldval; 6984246Sdes double dval[2]; 7084246Sdes} longDblUnion; 7175295Sdes 7275295Sdesstatic const double FPKINF = 1.0/0.0; 7375295Sdes 7475295Sdes/* Add two 'long double' values and return the result. */ 7575295Sdeslong double 7675295Sdes_xlqadd (double a, double b, double c, double d) 7775295Sdes{ 7875295Sdes longDblUnion z; 7975295Sdes double t, tau, u, FPR_zero, FPR_PosInf; 8075295Sdes 8184246Sdes FPR_zero = 0.0; 8284246Sdes FPR_PosInf = FPKINF; 8384246Sdes 8475295Sdes if (unlikely (a != a) || unlikely (c != c)) 8575295Sdes return a + c; /* NaN result. */ 8675295Sdes 8775295Sdes /* Ordered operands are arranged in order of their magnitudes. */ 8875295Sdes 8975295Sdes /* Switch inputs if |(c,d)| > |(a,b)|. */ 9075295Sdes if (fabs (c) > fabs (a)) 9175295Sdes { 9284246Sdes t = a; 9384386Sdes tau = b; 9484386Sdes a = c; 9584386Sdes b = d; 9675295Sdes c = t; 9775295Sdes d = tau; 9875295Sdes } 9975295Sdes 10075295Sdes /* b <- second largest magnitude double. */ 10175295Sdes if (fabs (c) > fabs (b)) 10275295Sdes { 10377998Sdes t = b; 10477998Sdes b = c; 10575295Sdes c = t; 10677998Sdes } 10775295Sdes 10888234Sdillon /* Thanks to commutivity, sum is invariant w.r.t. the next 10988234Sdillon conditional exchange. */ 11088234Sdillon tau = d + c; 11188234Sdillon 11288234Sdillon /* Order the smallest magnitude doubles. */ 11375295Sdes if (fabs (d) > fabs (c)) 11484246Sdes { 11577998Sdes t = c; 11684246Sdes c = d; 11775295Sdes d = t; 11884246Sdes } 11975295Sdes 12088234Sdillon t = (tau + b) + a; /* Sum values in ascending magnitude order. */ 12188234Sdillon 12275295Sdes /* Infinite or zero result. */ 12375295Sdes if (unlikely (t == FPR_zero) || unlikely (fabs (t) == FPR_PosInf)) 12477998Sdes return t; 12577998Sdes 12677998Sdes /* Usual case. */ 12777998Sdes tau = (((a-t) + b) + c) + d; 12877998Sdes u = t + tau; 12975295Sdes z.dval[0] = u; /* Final fixup for long double result. */ 13075295Sdes z.dval[1] = (t - u) + tau; 13175295Sdes return z.ldval; 13277998Sdes} 13384246Sdes 13484246Sdeslong double 13575295Sdes_xlqsub (double a, double b, double c, double d) 13677998Sdes{ 13775295Sdes return _xlqadd (a, b, -c, -d); 13877998Sdes} 13977998Sdes 14077998Sdeslong double 14175295Sdes_xlqmul (double a, double b, double c, double d) 14275295Sdes{ 14375295Sdes longDblUnion z; 14475295Sdes double t, tau, u, v, w, FPR_zero, FPR_PosInf; 14575295Sdes 14675295Sdes FPR_zero = 0.0; 14784246Sdes FPR_PosInf = FPKINF; 14875295Sdes 14975295Sdes t = a * c; /* Highest order double term. */ 15075295Sdes 15177998Sdes if (unlikely (t != t) || unlikely (t == FPR_zero) 15275295Sdes || unlikely (fabs (t) == FPR_PosInf)) 15375295Sdes return t; 15475295Sdes 15575295Sdes /* Finite nonzero result requires summing of terms of two highest 15675295Sdes orders. */ 15775295Sdes 15875295Sdes /* Use fused multiply-add to get low part of a * c. */ 15975295Sdes asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t)); 16077998Sdes v = a*d; 16177998Sdes w = b*c; 16275295Sdes tau += v + w; /* Add in other second-order terms. */ 16375295Sdes u = t + tau; 16475295Sdes 16584246Sdes /* Construct long double result. */ 16677998Sdes z.dval[0] = u; 16784246Sdes z.dval[1] = (t - u) + tau; 16884246Sdes return z.ldval; 16984246Sdes} 17084246Sdes 17184246Sdeslong double 17275295Sdes_xlqdiv (double a, double b, double c, double d) 17375295Sdes{ 17475295Sdes longDblUnion z; 17575295Sdes double s, sigma, t, tau, u, v, w, FPR_zero, FPR_PosInf; 17675295Sdes 17775295Sdes FPR_zero = 0.0; 17875295Sdes FPR_PosInf = FPKINF; 17975295Sdes 18075295Sdes t = a / c; /* highest order double term */ 18175295Sdes 18277998Sdes if (unlikely (t != t) || unlikely (t == FPR_zero) 18388234Sdillon || unlikely (fabs (t) == FPR_PosInf)) 18488234Sdillon return t; 18575295Sdes 18675295Sdes /* Finite nonzero result requires corrections to the highest order term. */ 18784246Sdes 18884246Sdes s = c * t; /* (s,sigma) = c*t exactly. */ 18984246Sdes w = -(-b + d * t); /* Written to get fnmsub for speed, but not 19084246Sdes numerically necessary. */ 19184246Sdes 19284246Sdes /* Use fused multiply-add to get low part of c * t. */ 19377998Sdes asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s)); 19484246Sdes v = a - s; 19577998Sdes 19684246Sdes tau = ((v-sigma)+w)/c; /* Correction to t. */ 19784246Sdes u = t + tau; 19877998Sdes 19975295Sdes /* Construct long double result. */ 20075295Sdes z.dval[0] = u; 20175295Sdes z.dval[1] = (t - u) + tau; 20284246Sdes return z.ldval; 20384246Sdes} 20484246Sdes 20584246Sdes#endif 20684246Sdes