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