1/* Copyright (C) 1989-2022 Free Software Foundation, Inc. 2 3This file is part of GCC. 4 5GCC is free software; you can redistribute it and/or modify it under 6the terms of the GNU General Public License as published by the Free 7Software Foundation; either version 3, or (at your option) any later 8version. 9 10GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11WARRANTY; without even the implied warranty of MERCHANTABILITY or 12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13for more details. 14 15Under Section 7 of GPL version 3, you are granted additional 16permissions described in the GCC Runtime Library Exception, version 173.1, as published by the Free Software Foundation. 18 19You should have received a copy of the GNU General Public License and 20a copy of the GCC Runtime Library Exception along with this program; 21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22<http://www.gnu.org/licenses/>. */ 23 24/* This is a temporary specialization of code from libgcc/libgcc2.c. */ 25 26#include "soft-fp.h" 27#include "quad-float128.h" 28 29#define COPYSIGN(x,y) __builtin_copysignf128 (x, y) 30#define INFINITY __builtin_inff128 () 31#define FABS __builtin_fabsf128 32#define isnan __builtin_isnan 33#define isinf __builtin_isinf 34#define isfinite __builtin_isfinite 35 36#if defined(FLOAT128_HW_INSNS) && !defined(__divkc3) 37#define __divkc3 __divkc3_sw 38#endif 39 40#ifndef __LONG_DOUBLE_IEEE128__ 41#define RBIG (__LIBGCC_KF_MAX__ / 2) 42#define RMIN (__LIBGCC_KF_MIN__) 43#define RMIN2 (__LIBGCC_KF_EPSILON__) 44#define RMINSCAL (1 / __LIBGCC_KF_EPSILON__) 45#define RMAX2 (RBIG * RMIN2) 46#else 47#define RBIG (__LIBGCC_TF_MAX__ / 2) 48#define RMIN (__LIBGCC_TF_MIN__) 49#define RMIN2 (__LIBGCC_TF_EPSILON__) 50#define RMINSCAL (1 / __LIBGCC_TF_EPSILON__) 51#define RMAX2 (RBIG * RMIN2) 52#endif 53 54TCtype 55__divkc3 (TFtype a, TFtype b, TFtype c, TFtype d) 56{ 57 TFtype denom, ratio, x, y; 58 TCtype res; 59 60 /* long double has significant potential underflow/overflow errors that 61 can be greatly reduced with a limited number of tests and adjustments. 62 */ 63 64 /* Scale by max(c,d) to reduce chances of denominator overflowing. */ 65 if (FABS (c) < FABS (d)) 66 { 67 /* Prevent underflow when denominator is near max representable. */ 68 if (FABS (d) >= RBIG) 69 { 70 a = a / 2; 71 b = b / 2; 72 c = c / 2; 73 d = d / 2; 74 } 75 /* Avoid overflow/underflow issues when c and d are small. 76 Scaling up helps avoid some underflows. 77 No new overflow possible since c&d < RMIN2. */ 78 if (FABS (d) < RMIN2) 79 { 80 a = a * RMINSCAL; 81 b = b * RMINSCAL; 82 c = c * RMINSCAL; 83 d = d * RMINSCAL; 84 } 85 else 86 { 87 if (((FABS (a) < RMIN) && (FABS (b) < RMAX2) && (FABS (d) < RMAX2)) 88 || ((FABS (b) < RMIN) && (FABS (a) < RMAX2) 89 && (FABS (d) < RMAX2))) 90 { 91 a = a * RMINSCAL; 92 b = b * RMINSCAL; 93 c = c * RMINSCAL; 94 d = d * RMINSCAL; 95 } 96 } 97 ratio = c / d; 98 denom = (c * ratio) + d; 99 /* Choose alternate order of computation if ratio is subnormal. */ 100 if (FABS (ratio) > RMIN) 101 { 102 x = ((a * ratio) + b) / denom; 103 y = ((b * ratio) - a) / denom; 104 } 105 else 106 { 107 x = ((c * (a / d)) + b) / denom; 108 y = ((c * (b / d)) - a) / denom; 109 } 110 } 111 else 112 { 113 /* Prevent underflow when denominator is near max representable. */ 114 if (FABS (c) >= RBIG) 115 { 116 a = a / 2; 117 b = b / 2; 118 c = c / 2; 119 d = d / 2; 120 } 121 /* Avoid overflow/underflow issues when both c and d are small. 122 Scaling up helps avoid some underflows. 123 No new overflow possible since both c&d are less than RMIN2. */ 124 if (FABS (c) < RMIN2) 125 { 126 a = a * RMINSCAL; 127 b = b * RMINSCAL; 128 c = c * RMINSCAL; 129 d = d * RMINSCAL; 130 } 131 else 132 { 133 if (((FABS (a) < RMIN) && (FABS (b) < RMAX2) && (FABS (c) < RMAX2)) 134 || ((FABS (b) < RMIN) && (FABS (a) < RMAX2) 135 && (FABS (c) < RMAX2))) 136 { 137 a = a * RMINSCAL; 138 b = b * RMINSCAL; 139 c = c * RMINSCAL; 140 d = d * RMINSCAL; 141 } 142 } 143 ratio = d / c; 144 denom = (d * ratio) + c; 145 /* Choose alternate order of computation if ratio is subnormal. */ 146 if (FABS (ratio) > RMIN) 147 { 148 x = ((b * ratio) + a) / denom; 149 y = (b - (a * ratio)) / denom; 150 } 151 else 152 { 153 x = (a + (d * (b / c))) / denom; 154 y = (b - (d * (a / c))) / denom; 155 } 156 } 157 158 /* Recover infinities and zeros that computed as NaN+iNaN; the only cases 159 are nonzero/zero, infinite/finite, and finite/infinite. */ 160 if (isnan (x) && isnan (y)) 161 { 162 if (c == 0.0 && d == 0.0 && (!isnan (a) || !isnan (b))) 163 { 164 x = COPYSIGN (INFINITY, c) * a; 165 y = COPYSIGN (INFINITY, c) * b; 166 } 167 else if ((isinf (a) || isinf (b)) && isfinite (c) && isfinite (d)) 168 { 169 a = COPYSIGN (isinf (a) ? 1 : 0, a); 170 b = COPYSIGN (isinf (b) ? 1 : 0, b); 171 x = INFINITY * (a * c + b * d); 172 y = INFINITY * (b * c - a * d); 173 } 174 else if ((isinf (c) || isinf (d)) && isfinite (a) && isfinite (b)) 175 { 176 c = COPYSIGN (isinf (c) ? 1 : 0, c); 177 d = COPYSIGN (isinf (d) ? 1 : 0, d); 178 x = 0.0 * (a * c + b * d); 179 y = 0.0 * (b * c - a * d); 180 } 181 } 182 183 __real__ res = x; 184 __imag__ res = y; 185 return res; 186} 187