98 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#include "../int_math.h"6#include "DD.h"7// Use DOUBLE_PRECISION because the soft-fp method we use is logb (on the upper8// half of the long doubles), even though this file defines complex division for9// 128-bit floats.10#define DOUBLE_PRECISION11#include "../fp_lib.h"12 13#if !defined(CRT_INFINITY) && defined(HUGE_VAL)14#define CRT_INFINITY HUGE_VAL15#endif // CRT_INFINITY16 17#define makeFinite(x) \18 { \19 (x).s.hi = crt_copysign(crt_isinf((x).s.hi) ? 1.0 : 0.0, (x).s.hi); \20 (x).s.lo = 0.0; \21 }22 23long double _Complex __divtc3(long double a, long double b, long double c,24 long double d) {25 DD cDD = {.ld = c};26 DD dDD = {.ld = d};27 28 int ilogbw = 0;29 const double logbw =30 __compiler_rt_logb(__compiler_rt_fmax(crt_fabs(cDD.s.hi),31 crt_fabs(dDD.s.hi)));32 33 if (crt_isfinite(logbw)) {34 ilogbw = (int)logbw;35 36 cDD.s.hi = __compiler_rt_scalbn(cDD.s.hi, -ilogbw);37 cDD.s.lo = __compiler_rt_scalbn(cDD.s.lo, -ilogbw);38 dDD.s.hi = __compiler_rt_scalbn(dDD.s.hi, -ilogbw);39 dDD.s.lo = __compiler_rt_scalbn(dDD.s.lo, -ilogbw);40 }41 42 const long double denom =43 __gcc_qadd(__gcc_qmul(cDD.ld, cDD.ld), __gcc_qmul(dDD.ld, dDD.ld));44 const long double realNumerator =45 __gcc_qadd(__gcc_qmul(a, cDD.ld), __gcc_qmul(b, dDD.ld));46 const long double imagNumerator =47 __gcc_qsub(__gcc_qmul(b, cDD.ld), __gcc_qmul(a, dDD.ld));48 49 DD real = {.ld = __gcc_qdiv(realNumerator, denom)};50 DD imag = {.ld = __gcc_qdiv(imagNumerator, denom)};51 52 real.s.hi = __compiler_rt_scalbn(real.s.hi, -ilogbw);53 real.s.lo = __compiler_rt_scalbn(real.s.lo, -ilogbw);54 imag.s.hi = __compiler_rt_scalbn(imag.s.hi, -ilogbw);55 imag.s.lo = __compiler_rt_scalbn(imag.s.lo, -ilogbw);56 57 if (crt_isnan(real.s.hi) && crt_isnan(imag.s.hi)) {58 DD aDD = {.ld = a};59 DD bDD = {.ld = b};60 DD rDD = {.ld = denom};61 62 if ((rDD.s.hi == 0.0) && (!crt_isnan(aDD.s.hi) || !crt_isnan(bDD.s.hi))) {63 real.s.hi = crt_copysign(CRT_INFINITY, cDD.s.hi) * aDD.s.hi;64 real.s.lo = 0.0;65 imag.s.hi = crt_copysign(CRT_INFINITY, cDD.s.hi) * bDD.s.hi;66 imag.s.lo = 0.0;67 }68 69 else if ((crt_isinf(aDD.s.hi) || crt_isinf(bDD.s.hi)) &&70 crt_isfinite(cDD.s.hi) && crt_isfinite(dDD.s.hi)) {71 makeFinite(aDD);72 makeFinite(bDD);73 real.s.hi = CRT_INFINITY * (aDD.s.hi * cDD.s.hi + bDD.s.hi * dDD.s.hi);74 real.s.lo = 0.0;75 imag.s.hi = CRT_INFINITY * (bDD.s.hi * cDD.s.hi - aDD.s.hi * dDD.s.hi);76 imag.s.lo = 0.0;77 }78 79 else if ((crt_isinf(cDD.s.hi) || crt_isinf(dDD.s.hi)) &&80 crt_isfinite(aDD.s.hi) && crt_isfinite(bDD.s.hi)) {81 makeFinite(cDD);82 makeFinite(dDD);83 real.s.hi =84 crt_copysign(0.0, (aDD.s.hi * cDD.s.hi + bDD.s.hi * dDD.s.hi));85 real.s.lo = 0.0;86 imag.s.hi =87 crt_copysign(0.0, (bDD.s.hi * cDD.s.hi - aDD.s.hi * dDD.s.hi));88 imag.s.lo = 0.0;89 }90 }91 92 long double _Complex z;93 __real__ z = real.ld;94 __imag__ z = imag.ld;95 96 return z;97}98