| 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_qdiv(long double x, long double y); |
|---|
| 6 | // This file implements the PowerPC 128-bit double-double division operation. |
|---|
| 7 | // This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!) |
|---|
| 8 | |
|---|
| 9 | #include "DD.h" |
|---|
| 10 | |
|---|
| 11 | long double __gcc_qdiv(long double a, long double b) { |
|---|
| 12 | static const uint32_t infinityHi = UINT32_C(0x7ff00000); |
|---|
| 13 | DD dst = {.ld = a}, src = {.ld = b}; |
|---|
| 14 | |
|---|
| 15 | register double x = dst.s.hi, x1 = dst.s.lo, y = src.s.hi, y1 = src.s.lo; |
|---|
| 16 | |
|---|
| 17 | double yHi, yLo, qHi, qLo; |
|---|
| 18 | double yq, tmp, q; |
|---|
| 19 | |
|---|
| 20 | q = x / y; |
|---|
| 21 | |
|---|
| 22 | // Detect special cases |
|---|
| 23 | if (q == 0.0) { |
|---|
| 24 | dst.s.hi = q; |
|---|
| 25 | dst.s.lo = 0.0; |
|---|
| 26 | return dst.ld; |
|---|
| 27 | } |
|---|
| 28 | |
|---|
| 29 | const doublebits qBits = {.d = q}; |
|---|
| 30 | if (((uint32_t)(qBits.x >> 32) & infinityHi) == infinityHi) { |
|---|
| 31 | dst.s.hi = q; |
|---|
| 32 | dst.s.lo = 0.0; |
|---|
| 33 | return dst.ld; |
|---|
| 34 | } |
|---|
| 35 | |
|---|
| 36 | yHi = high26bits(x: y); |
|---|
| 37 | qHi = high26bits(x: q); |
|---|
| 38 | |
|---|
| 39 | yq = y * q; |
|---|
| 40 | yLo = y - yHi; |
|---|
| 41 | qLo = q - qHi; |
|---|
| 42 | |
|---|
| 43 | tmp = LOWORDER(yq, yHi, yLo, qHi, qLo); |
|---|
| 44 | tmp = (x - yq) - tmp; |
|---|
| 45 | tmp = ((tmp + x1) - y1 * q) / y; |
|---|
| 46 | x = q + tmp; |
|---|
| 47 | |
|---|
| 48 | dst.s.lo = (q - x) + tmp; |
|---|
| 49 | dst.s.hi = x; |
|---|
| 50 | |
|---|
| 51 | return dst.ld; |
|---|
| 52 | } |
|---|
| 53 | |
|---|
