xref: /llvm-project/compiler-rt/lib/builtins/ppc/gcc_qdiv.c (revision 0ba22f51d128bee9d69756c56c4678097270e10b)
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