| /* |
| * IBM Accurate Mathematical Library |
| * written by International Business Machines Corp. |
| * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU Lesser General Public License as published by |
| * the Free Software Foundation; either version 2.1 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU Lesser General Public License for more details. |
| * |
| * You should have received a copy of the GNU Lesser General Public License |
| * along with this program; if not, see <http://www.gnu.org/licenses/>. |
| */ |
| /*********************************************************************/ |
| /* MODULE_NAME: uroot.c */ |
| /* */ |
| /* FUNCTION: usqrt */ |
| /* */ |
| /* FILES NEEDED: dla.h endian.h mydefs.h */ |
| /* uroot.tbl */ |
| /* */ |
| /* An ultimate sqrt routine. Given an IEEE double machine number x */ |
| /* it computes the correctly rounded (to nearest) value of square */ |
| /* root of x. */ |
| /* Assumption: Machine arithmetic operations are performed in */ |
| /* round to nearest mode of IEEE 754 standard. */ |
| /* */ |
| /*********************************************************************/ |
| |
| #include "endian.h" |
| #include "mydefs.h" |
| #include <dla.h> |
| #include "MathLib.h" |
| #include "root.tbl" |
| #include <math_private.h> |
| |
| /*********************************************************************/ |
| /* An ultimate sqrt routine. Given an IEEE double machine number x */ |
| /* it computes the correctly rounded (to nearest) value of square */ |
| /* root of x. */ |
| /*********************************************************************/ |
| double |
| __ieee754_sqrt (double x) |
| { |
| static const double |
| rt0 = 9.99999999859990725855365213134618E-01, |
| rt1 = 4.99999999495955425917856814202739E-01, |
| rt2 = 3.75017500867345182581453026130850E-01, |
| rt3 = 3.12523626554518656309172508769531E-01; |
| static const double big = 134217728.0; |
| double y, t, del, res, res1, hy, z, zz, p, hx, tx, ty, s; |
| mynumber a, c = { { 0, 0 } }; |
| int4 k; |
| |
| a.x = x; |
| k = a.i[HIGH_HALF]; |
| a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000; |
| t = inroot[(k & 0x001fffff) >> 14]; |
| s = a.x; |
| /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ |
| if (k > 0x000fffff && k < 0x7ff00000) |
| { |
| int rm = __fegetround (); |
| fenv_t env; |
| libc_feholdexcept_setround (&env, FE_TONEAREST); |
| double ret; |
| y = 1.0 - t * (t * s); |
| t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3))); |
| c.i[HIGH_HALF] = 0x20000000 + ((k & 0x7fe00000) >> 1); |
| y = t * s; |
| hy = (y + big) - big; |
| del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy)); |
| res = y + del; |
| if (res == (res + 1.002 * ((y - res) + del))) |
| ret = res * c.x; |
| else |
| { |
| res1 = res + 1.5 * ((y - res) + del); |
| EMULV (res, res1, z, zz, p, hx, tx, hy, ty); /* (z+zz)=res*res1 */ |
| res = ((((z - s) + zz) < 0) ? max (res, res1) : |
| min (res, res1)); |
| ret = res * c.x; |
| } |
| math_force_eval (ret); |
| libc_fesetenv (&env); |
| double dret = x / ret; |
| if (dret != ret) |
| { |
| double force_inexact = 1.0 / 3.0; |
| math_force_eval (force_inexact); |
| /* The square root is inexact, ret is the round-to-nearest |
| value which may need adjusting for other rounding |
| modes. */ |
| switch (rm) |
| { |
| #ifdef FE_UPWARD |
| case FE_UPWARD: |
| if (dret > ret) |
| ret = (res + 0x1p-1022) * c.x; |
| break; |
| #endif |
| |
| #ifdef FE_DOWNWARD |
| case FE_DOWNWARD: |
| #endif |
| #ifdef FE_TOWARDZERO |
| case FE_TOWARDZERO: |
| #endif |
| #if defined FE_DOWNWARD || defined FE_TOWARDZERO |
| if (dret < ret) |
| ret = (res - 0x1p-1022) * c.x; |
| break; |
| #endif |
| |
| default: |
| break; |
| } |
| } |
| /* Otherwise (x / ret == ret), either the square root was exact or |
| the division was inexact. */ |
| return ret; |
| } |
| else |
| { |
| if ((k & 0x7ff00000) == 0x7ff00000) |
| return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ |
| if (x == 0) |
| return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */ |
| if (k < 0) |
| return (x - x) / (x - x); /* sqrt(-ve)=sNaN */ |
| return 0x1p-256 * __ieee754_sqrt (x * 0x1p512); |
| } |
| } |
| strong_alias (__ieee754_sqrt, __sqrt_finite) |