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