53 lines · c
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-exception4 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 11long 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 cases23 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