hypotf.c
1 /* 2 * Copyright (C) 2008-2020 Advanced Micro Devices, Inc. All rights reserved. 3 * 4 * Redistribution and use in source and binary forms, with or without modification, 5 * are permitted provided that the following conditions are met: 6 * 1. Redistributions of source code must retain the above copyright notice, 7 * this list of conditions and the following disclaimer. 8 * 2. Redistributions in binary form must reproduce the above copyright notice, 9 * this list of conditions and the following disclaimer in the documentation 10 * and/or other materials provided with the distribution. 11 * 3. Neither the name of the copyright holder nor the names of its contributors 12 * may be used to endorse or promote products derived from this software without 13 * specific prior written permission. 14 * 15 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 17 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, 19 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 20 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, 21 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 22 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 23 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 24 * POSSIBILITY OF SUCH DAMAGE. 25 * 26 */ 27 28 #include "libm_amd.h" 29 #include "libm_util_amd.h" 30 #include "libm_special.h" 31 32 float FN_PROTOTYPE_REF(hypotf)(float x, float y) 33 { 34 /* Returns sqrt(x*x + y*y) with no overflow or underflow unless 35 the result warrants it */ 36 37 /* Do intermediate computations in double precision 38 and use sqrt instruction from chip if available. */ 39 double dx, dy, dr, retval; 40 41 /* The largest finite float, stored as a double */ 42 double large; 43 44 int x_is_nan,y_is_nan; 45 unsigned int ux, uy, avx, avy; 46 UT32 val; 47 GET_BITS_SP32(x, avx); 48 avx &= ~SIGNBIT_SP32; 49 GET_BITS_SP32(y, avy); 50 avy &= ~SIGNBIT_SP32; 51 ux = (avx >> EXPSHIFTBITS_SP32); 52 uy = (avy >> EXPSHIFTBITS_SP32); 53 x_is_nan = (avx > 0x7f800000); 54 y_is_nan = (avy > 0x7f800000); 55 val.u32 = PINFBITPATT_SP32; 56 57 if (ux == BIASEDEMAX_SP32 + 1 || uy == BIASEDEMAX_SP32 + 1) 58 { 59 /* One or both of the arguments are NaN or infinity. The 60 result will also be NaN or infinity. */ 61 if (((ux == BIASEDEMAX_SP32 + 1) && !(avx & MANTBITS_SP32)) || 62 ((uy == BIASEDEMAX_SP32 + 1) && !(avy & MANTBITS_SP32))) 63 /* x or y is infinity. ISO C99 defines that we must 64 return +infinity, even if the other argument is NaN. 65 Note that the computation of x*x + y*y above will already 66 have raised invalid if either x or y is a signalling NaN. */ 67 { 68 if(x_is_nan) 69 { 70 #ifdef WINDOWS 71 return val.f32 ;// __amd_handle_errorf("hypotf", __amd_hypot, val.u32, _DOMAIN, 0, EDOM, x, y, 2); 72 #else 73 if(!(avx & 0x00400000)) //x is snan 74 return __amd_handle_errorf("hypotf", __amd_hypot, val.u32, _DOMAIN, AMD_F_INVALID, EDOM, x, y, 2); 75 #endif 76 } 77 if(y_is_nan) 78 { 79 #ifdef WINDOWS 80 return val.f32 ; //__amd_handle_errorf("hypotf", __amd_hypot, val.u32, _DOMAIN, 0, EDOM, x, y, 2); 81 #else 82 if(!(avy & 0x00400000)) //y is snan 83 return __amd_handle_errorf("hypotf", __amd_hypot, val.u32, _DOMAIN, AMD_F_INVALID, EDOM, x, y, 2); 84 #endif 85 } 86 return val.f32; 87 } 88 89 if(x_is_nan || y_is_nan) 90 { 91 /* One or both of x or y is NaN, and neither is infinity. 92 Raise invalid if it's a signalling NaN */ 93 if(x_is_nan) 94 { 95 val.f32 = x; 96 #ifdef WINDOWS 97 return __amd_handle_errorf("hypotf", __amd_hypot, val.u32|0x00400000, _DOMAIN, 0, EDOM, x, y, 2); 98 #else 99 if(!(avx & 0x00400000)) //x is snan 100 return __amd_handle_errorf("hypotf", __amd_hypot, val.u32|0x00400000, _DOMAIN, AMD_F_INVALID, EDOM, x, y, 2); 101 #endif 102 } 103 if(y_is_nan) 104 { 105 val.f32 = y; 106 #ifdef WINDOWS 107 return __amd_handle_errorf("hypotf", __amd_hypot, val.u32|0x00400000, _DOMAIN, 0, EDOM, x, y, 2); 108 #else 109 if(!(avy & 0x00400000)) //y is snan 110 return __amd_handle_errorf("hypotf", __amd_hypot, val.u32|0x00400000, _DOMAIN, AMD_F_INVALID, EDOM, x, y, 2); 111 #endif 112 } 113 } 114 } 115 large = 3.40282346638528859812e+38; /* 0x47efffffe0000000 */ 116 dx = x, dy = y; 117 if(avx == 0x0) // to prevent underflow exception when x is small 118 { 119 val.u32 = avy; 120 return val.f32; 121 } 122 if(avy == 0x0)// to prevent underflow exception when y is small 123 { 124 val.u32 = avx; 125 return val.f32; 126 } 127 dr = (dx*dx + dy*dy); 128 129 #ifdef WINDOWS 130 /* VC++ intrinsic call */ 131 _mm_store_sd(&retval, _mm_sqrt_sd(_mm_setzero_pd(), _mm_load_sd(&dr))); 132 #else 133 /* Hammer sqrt instruction */ 134 asm volatile ("sqrtsd %1, %0" : "=x" (retval) : "x" (dr)); 135 #endif 136 137 if (retval > large) 138 return __amd_handle_errorf("hypotf", __amd_hypot, PINFBITPATT_SP32, _OVERFLOW, AMD_F_INEXACT|AMD_F_OVERFLOW, ERANGE, x, y, 2); 139 else if((x !=0.0) && (y!=0)) 140 { 141 val.f32 = (float)(retval); 142 val.u32 >>= EXPSHIFTBITS_SP32; 143 if(val.u32 == 0x0) 144 { 145 val.f32 = (float)retval; 146 return __amd_handle_errorf("hypotf", __amd_hypot, val.u32, _UNDERFLOW, AMD_F_INEXACT|AMD_F_UNDERFLOW, ERANGE, x, y, 2); 147 } 148 else 149 return (float)retval; 150 } 151 else 152 return (float)retval; 153 } 154 155