xref: /freebsd/contrib/llvm-project/compiler-rt/lib/builtins/ppc/gcc_qdiv.c (revision 95ee2897e98f5d444f26ed2334cc7c439f9c16c6)
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(y);
37   qHi = high26bits(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