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_qdiv(long double x, long double y);
6*0b57cec5SDimitry Andric // This file implements the PowerPC 128-bit double-double division 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_qdiv(long double a,long double b)11*0b57cec5SDimitry Andric long double __gcc_qdiv(long double a, long double b) {
12*0b57cec5SDimitry Andric static const uint32_t infinityHi = UINT32_C(0x7ff00000);
13*0b57cec5SDimitry Andric DD dst = {.ld = a}, src = {.ld = b};
14*0b57cec5SDimitry Andric
15*0b57cec5SDimitry Andric register double x = dst.s.hi, x1 = dst.s.lo, y = src.s.hi, y1 = src.s.lo;
16*0b57cec5SDimitry Andric
17*0b57cec5SDimitry Andric double yHi, yLo, qHi, qLo;
18*0b57cec5SDimitry Andric double yq, tmp, q;
19*0b57cec5SDimitry Andric
20*0b57cec5SDimitry Andric q = x / y;
21*0b57cec5SDimitry Andric
22*0b57cec5SDimitry Andric // Detect special cases
23*0b57cec5SDimitry Andric if (q == 0.0) {
24*0b57cec5SDimitry Andric dst.s.hi = q;
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 qBits = {.d = q};
30*0b57cec5SDimitry Andric if (((uint32_t)(qBits.x >> 32) & infinityHi) == infinityHi) {
31*0b57cec5SDimitry Andric dst.s.hi = q;
32*0b57cec5SDimitry Andric dst.s.lo = 0.0;
33*0b57cec5SDimitry Andric return dst.ld;
34*0b57cec5SDimitry Andric }
35*0b57cec5SDimitry Andric
36*0b57cec5SDimitry Andric yHi = high26bits(y);
37*0b57cec5SDimitry Andric qHi = high26bits(q);
38*0b57cec5SDimitry Andric
39*0b57cec5SDimitry Andric yq = y * q;
40*0b57cec5SDimitry Andric yLo = y - yHi;
41*0b57cec5SDimitry Andric qLo = q - qHi;
42*0b57cec5SDimitry Andric
43*0b57cec5SDimitry Andric tmp = LOWORDER(yq, yHi, yLo, qHi, qLo);
44*0b57cec5SDimitry Andric tmp = (x - yq) - tmp;
45*0b57cec5SDimitry Andric tmp = ((tmp + x1) - y1 * q) / y;
46*0b57cec5SDimitry Andric x = q + tmp;
47*0b57cec5SDimitry Andric
48*0b57cec5SDimitry Andric dst.s.lo = (q - x) + tmp;
49*0b57cec5SDimitry Andric dst.s.hi = x;
50*0b57cec5SDimitry Andric
51*0b57cec5SDimitry Andric return dst.ld;
52*0b57cec5SDimitry Andric }
53