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