143615Sdcs// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
243615Sdcs// See https://llvm.org/LICENSE.txt for license information.
343615Sdcs// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
443615Sdcs
543615Sdcs// long double __gcc_qdiv(long double x, long double y);
643615Sdcs// This file implements the PowerPC 128-bit double-double division operation.
743615Sdcs// This implementation is shamelessly cribbed from Apple's DDRT, circa 1993(!)
843615Sdcs
943615Sdcs#include "DD.h"
1043615Sdcs
1143615Sdcslong double __gcc_qdiv(long double a, long double b) {
1243615Sdcs  static const uint32_t infinityHi = UINT32_C(0x7ff00000);
1343615Sdcs  DD dst = {.ld = a}, src = {.ld = b};
1443615Sdcs
1543615Sdcs  register double x = dst.s.hi, x1 = dst.s.lo, y = src.s.hi, y1 = src.s.lo;
1643615Sdcs
1743615Sdcs  double yHi, yLo, qHi, qLo;
1843615Sdcs  double yq, tmp, q;
1943615Sdcs
2043615Sdcs  q = x / y;
2143615Sdcs
2243615Sdcs  // Detect special cases
2343615Sdcs  if (q == 0.0) {
2443615Sdcs    dst.s.hi = q;
2544774Sdcs    dst.s.lo = 0.0;
2643615Sdcs    return dst.ld;
2743615Sdcs  }
2843615Sdcs
2943615Sdcs  const doublebits qBits = {.d = q};
3043615Sdcs  if (((uint32_t)(qBits.x >> 32) & infinityHi) == infinityHi) {
3143615Sdcs    dst.s.hi = q;
3243615Sdcs    dst.s.lo = 0.0;
3344774Sdcs    return dst.ld;
3443615Sdcs  }
3543615Sdcs
3643615Sdcs  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