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