1/* Complex square root of a float type. 2 Copyright (C) 1997-2018 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>. 5 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. 6 7 The GNU C Library is free software; you can redistribute it and/or 8 modify it under the terms of the GNU Lesser General Public 9 License as published by the Free Software Foundation; either 10 version 2.1 of the License, or (at your option) any later version. 11 12 The GNU C Library is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15 Lesser General Public License for more details. 16 17 You should have received a copy of the GNU Lesser General Public 18 License along with the GNU C Library; if not, see 19 <http://www.gnu.org/licenses/>. */ 20 21#include "quadmath-imp.h" 22 23__complex128 24csqrtq (__complex128 x) 25{ 26 __complex128 res; 27 int rcls = fpclassifyq (__real__ x); 28 int icls = fpclassifyq (__imag__ x); 29 30 if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE)) 31 { 32 if (icls == QUADFP_INFINITE) 33 { 34 __real__ res = HUGE_VALQ; 35 __imag__ res = __imag__ x; 36 } 37 else if (rcls == QUADFP_INFINITE) 38 { 39 if (__real__ x < 0) 40 { 41 __real__ res = icls == QUADFP_NAN ? nanq ("") : 0; 42 __imag__ res = copysignq (HUGE_VALQ, __imag__ x); 43 } 44 else 45 { 46 __real__ res = __real__ x; 47 __imag__ res = (icls == QUADFP_NAN 48 ? nanq ("") : copysignq (0, __imag__ x)); 49 } 50 } 51 else 52 { 53 __real__ res = nanq (""); 54 __imag__ res = nanq (""); 55 } 56 } 57 else 58 { 59 if (__glibc_unlikely (icls == QUADFP_ZERO)) 60 { 61 if (__real__ x < 0) 62 { 63 __real__ res = 0; 64 __imag__ res = copysignq (sqrtq (-__real__ x), __imag__ x); 65 } 66 else 67 { 68 __real__ res = fabsq (sqrtq (__real__ x)); 69 __imag__ res = copysignq (0, __imag__ x); 70 } 71 } 72 else if (__glibc_unlikely (rcls == QUADFP_ZERO)) 73 { 74 __float128 r; 75 if (fabsq (__imag__ x) >= 2 * FLT128_MIN) 76 r = sqrtq (0.5Q * fabsq (__imag__ x)); 77 else 78 r = 0.5Q * sqrtq (2 * fabsq (__imag__ x)); 79 80 __real__ res = r; 81 __imag__ res = copysignq (r, __imag__ x); 82 } 83 else 84 { 85 __float128 d, r, s; 86 int scale = 0; 87 88 if (fabsq (__real__ x) > FLT128_MAX / 4) 89 { 90 scale = 1; 91 __real__ x = scalbnq (__real__ x, -2 * scale); 92 __imag__ x = scalbnq (__imag__ x, -2 * scale); 93 } 94 else if (fabsq (__imag__ x) > FLT128_MAX / 4) 95 { 96 scale = 1; 97 if (fabsq (__real__ x) >= 4 * FLT128_MIN) 98 __real__ x = scalbnq (__real__ x, -2 * scale); 99 else 100 __real__ x = 0; 101 __imag__ x = scalbnq (__imag__ x, -2 * scale); 102 } 103 else if (fabsq (__real__ x) < 2 * FLT128_MIN 104 && fabsq (__imag__ x) < 2 * FLT128_MIN) 105 { 106 scale = -((FLT128_MANT_DIG + 1) / 2); 107 __real__ x = scalbnq (__real__ x, -2 * scale); 108 __imag__ x = scalbnq (__imag__ x, -2 * scale); 109 } 110 111 d = hypotq (__real__ x, __imag__ x); 112 /* Use the identity 2 Re res Im res = Im x 113 to avoid cancellation error in d +/- Re x. */ 114 if (__real__ x > 0) 115 { 116 r = sqrtq (0.5Q * (d + __real__ x)); 117 if (scale == 1 && fabsq (__imag__ x) < 1) 118 { 119 /* Avoid possible intermediate underflow. */ 120 s = __imag__ x / r; 121 r = scalbnq (r, scale); 122 scale = 0; 123 } 124 else 125 s = 0.5Q * (__imag__ x / r); 126 } 127 else 128 { 129 s = sqrtq (0.5Q * (d - __real__ x)); 130 if (scale == 1 && fabsq (__imag__ x) < 1) 131 { 132 /* Avoid possible intermediate underflow. */ 133 r = fabsq (__imag__ x / s); 134 s = scalbnq (s, scale); 135 scale = 0; 136 } 137 else 138 r = fabsq (0.5Q * (__imag__ x / s)); 139 } 140 141 if (scale) 142 { 143 r = scalbnq (r, scale); 144 s = scalbnq (s, scale); 145 } 146 147 math_check_force_underflow (r); 148 math_check_force_underflow (s); 149 150 __real__ res = r; 151 __imag__ res = copysignq (s, __imag__ x); 152 } 153 } 154 155 return res; 156} 157