1214152Sed/* This file is distributed under the University of Illinois Open Source 2214152Sed * License. See LICENSE.TXT for details. 3214152Sed */ 4214152Sed 5214152Sed/* long double __gcc_qmul(long double x, long double y); 6214152Sed * This file implements the PowerPC 128-bit double-double multiply operation. 7214152Sed * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) 8214152Sed */ 9214152Sed 10214152Sed#include "DD.h" 11214152Sed 12214152Sedlong double __gcc_qmul(long double x, long double y) 13214152Sed{ 14214152Sed static const uint32_t infinityHi = UINT32_C(0x7ff00000); 15214152Sed DD dst = { .ld = x }, src = { .ld = y }; 16214152Sed 17214152Sed register double A = dst.s.hi, a = dst.s.lo, 18214152Sed B = src.s.hi, b = src.s.lo; 19214152Sed 20214152Sed double aHi, aLo, bHi, bLo; 21214152Sed double ab, tmp, tau; 22214152Sed 23214152Sed ab = A * B; 24214152Sed 25214152Sed /* Detect special cases */ 26214152Sed if (ab == 0.0) { 27214152Sed dst.s.hi = ab; 28214152Sed dst.s.lo = 0.0; 29214152Sed return dst.ld; 30214152Sed } 31214152Sed 32214152Sed const doublebits abBits = { .d = ab }; 33214152Sed if (((uint32_t)(abBits.x >> 32) & infinityHi) == infinityHi) { 34214152Sed dst.s.hi = ab; 35214152Sed dst.s.lo = 0.0; 36214152Sed return dst.ld; 37214152Sed } 38214152Sed 39214152Sed /* Generic cases handled here. */ 40214152Sed aHi = high26bits(A); 41214152Sed bHi = high26bits(B); 42214152Sed aLo = A - aHi; 43214152Sed bLo = B - bHi; 44214152Sed 45214152Sed tmp = LOWORDER(ab, aHi, aLo, bHi, bLo); 46214152Sed tmp += (A * b + a * B); 47214152Sed tau = ab + tmp; 48214152Sed 49214152Sed dst.s.lo = (ab - tau) + tmp; 50214152Sed dst.s.hi = tau; 51214152Sed 52214152Sed return dst.ld; 53214152Sed} 54