1276789Sdim/* This file is distributed under the University of Illinois Open Source 2276789Sdim * License. See LICENSE.TXT for details. 3276789Sdim */ 4276789Sdim 5276789Sdim/* long double __gcc_qmul(long double x, long double y); 6276789Sdim * This file implements the PowerPC 128-bit double-double multiply operation. 7276789Sdim * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) 8276789Sdim */ 9276789Sdim 10276789Sdim#include "DD.h" 11276789Sdim 12276789Sdimlong double __gcc_qmul(long double x, long double y) 13276789Sdim{ 14276789Sdim static const uint32_t infinityHi = UINT32_C(0x7ff00000); 15276789Sdim DD dst = { .ld = x }, src = { .ld = y }; 16276789Sdim 17276789Sdim register double A = dst.s.hi, a = dst.s.lo, 18276789Sdim B = src.s.hi, b = src.s.lo; 19276789Sdim 20276789Sdim double aHi, aLo, bHi, bLo; 21276789Sdim double ab, tmp, tau; 22276789Sdim 23276789Sdim ab = A * B; 24276789Sdim 25276789Sdim /* Detect special cases */ 26276789Sdim if (ab == 0.0) { 27276789Sdim dst.s.hi = ab; 28276789Sdim dst.s.lo = 0.0; 29276789Sdim return dst.ld; 30276789Sdim } 31276789Sdim 32276789Sdim const doublebits abBits = { .d = ab }; 33276789Sdim if (((uint32_t)(abBits.x >> 32) & infinityHi) == infinityHi) { 34276789Sdim dst.s.hi = ab; 35276789Sdim dst.s.lo = 0.0; 36276789Sdim return dst.ld; 37276789Sdim } 38276789Sdim 39276789Sdim /* Generic cases handled here. */ 40276789Sdim aHi = high26bits(A); 41276789Sdim bHi = high26bits(B); 42276789Sdim aLo = A - aHi; 43276789Sdim bLo = B - bHi; 44276789Sdim 45276789Sdim tmp = LOWORDER(ab, aHi, aLo, bHi, bLo); 46276789Sdim tmp += (A * b + a * B); 47276789Sdim tau = ab + tmp; 48276789Sdim 49276789Sdim dst.s.lo = (ab - tau) + tmp; 50276789Sdim dst.s.hi = tau; 51276789Sdim 52276789Sdim return dst.ld; 53276789Sdim} 54