1*0b57cec5SDimitry Andric // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
2*0b57cec5SDimitry Andric // See https://llvm.org/LICENSE.txt for license information.
3*0b57cec5SDimitry Andric // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
4*0b57cec5SDimitry Andric
5*0b57cec5SDimitry Andric // long double __gcc_qmul(long double x, long double y);
6*0b57cec5SDimitry Andric // This file implements the PowerPC 128-bit double-double multiply operation.
7*0b57cec5SDimitry Andric // This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!)
8*0b57cec5SDimitry Andric
9*0b57cec5SDimitry Andric #include "DD.h"
10*0b57cec5SDimitry Andric
__gcc_qmul(long double x,long double y)11*0b57cec5SDimitry Andric long double __gcc_qmul(long double x, long double y) {
12*0b57cec5SDimitry Andric static const uint32_t infinityHi = UINT32_C(0x7ff00000);
13*0b57cec5SDimitry Andric DD dst = {.ld = x}, src = {.ld = y};
14*0b57cec5SDimitry Andric
15*0b57cec5SDimitry Andric register double A = dst.s.hi, a = dst.s.lo, B = src.s.hi, b = src.s.lo;
16*0b57cec5SDimitry Andric
17*0b57cec5SDimitry Andric double aHi, aLo, bHi, bLo;
18*0b57cec5SDimitry Andric double ab, tmp, tau;
19*0b57cec5SDimitry Andric
20*0b57cec5SDimitry Andric ab = A * B;
21*0b57cec5SDimitry Andric
22*0b57cec5SDimitry Andric // Detect special cases
23*0b57cec5SDimitry Andric if (ab == 0.0) {
24*0b57cec5SDimitry Andric dst.s.hi = ab;
25*0b57cec5SDimitry Andric dst.s.lo = 0.0;
26*0b57cec5SDimitry Andric return dst.ld;
27*0b57cec5SDimitry Andric }
28*0b57cec5SDimitry Andric
29*0b57cec5SDimitry Andric const doublebits abBits = {.d = ab};
30*0b57cec5SDimitry Andric if (((uint32_t)(abBits.x >> 32) & infinityHi) == infinityHi) {
31*0b57cec5SDimitry Andric dst.s.hi = ab;
32*0b57cec5SDimitry Andric dst.s.lo = 0.0;
33*0b57cec5SDimitry Andric return dst.ld;
34*0b57cec5SDimitry Andric }
35*0b57cec5SDimitry Andric
36*0b57cec5SDimitry Andric // Generic cases handled here.
37*0b57cec5SDimitry Andric aHi = high26bits(A);
38*0b57cec5SDimitry Andric bHi = high26bits(B);
39*0b57cec5SDimitry Andric aLo = A - aHi;
40*0b57cec5SDimitry Andric bLo = B - bHi;
41*0b57cec5SDimitry Andric
42*0b57cec5SDimitry Andric tmp = LOWORDER(ab, aHi, aLo, bHi, bLo);
43*0b57cec5SDimitry Andric tmp += (A * b + a * B);
44*0b57cec5SDimitry Andric tau = ab + tmp;
45*0b57cec5SDimitry Andric
46*0b57cec5SDimitry Andric dst.s.lo = (ab - tau) + tmp;
47*0b57cec5SDimitry Andric dst.s.hi = tau;
48*0b57cec5SDimitry Andric
49*0b57cec5SDimitry Andric return dst.ld;
50*0b57cec5SDimitry Andric }
51