177ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao/* This file is distributed under the University of Illinois Open Source 277ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao * License. See LICENSE.TXT for details. 377ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao */ 477ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 577ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao/* long double __gcc_qmul(long double x, long double y); 677ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao * This file implements the PowerPC 128-bit double-double multiply operation. 777ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao * This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) 877ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao */ 977ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 1077ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao#include "DD.h" 1177ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 1277ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liaolong double __gcc_qmul(long double x, long double y) 1377ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao{ 1477ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao static const uint32_t infinityHi = UINT32_C(0x7ff00000); 1577ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao DD dst = { .ld = x }, src = { .ld = y }; 1677ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 1777ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao register double A = dst.s.hi, a = dst.s.lo, 1877ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao B = src.s.hi, b = src.s.lo; 1977ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 2077ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao double aHi, aLo, bHi, bLo; 2177ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao double ab, tmp, tau; 2277ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 2377ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao ab = A * B; 2477ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 2577ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao /* Detect special cases */ 2677ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao if (ab == 0.0) { 2777ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao dst.s.hi = ab; 2877ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao dst.s.lo = 0.0; 2977ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao return dst.ld; 3077ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao } 3177ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 3277ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao const doublebits abBits = { .d = ab }; 3377ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao if (((uint32_t)(abBits.x >> 32) & infinityHi) == infinityHi) { 3477ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao dst.s.hi = ab; 3577ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao dst.s.lo = 0.0; 3677ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao return dst.ld; 3777ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao } 3877ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 3977ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao /* Generic cases handled here. */ 4077ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao aHi = high26bits(A); 4177ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao bHi = high26bits(B); 4277ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao aLo = A - aHi; 4377ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao bLo = B - bHi; 4477ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 4577ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao tmp = LOWORDER(ab, aHi, aLo, bHi, bLo); 4677ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao tmp += (A * b + a * B); 4777ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao tau = ab + tmp; 4877ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 4977ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao dst.s.lo = (ab - tau) + tmp; 5077ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao dst.s.hi = tau; 5177ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao 5277ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao return dst.ld; 5377ed6142daed1e068fbda64405d0de9845e40e1Shih-wei Liao} 54