1/* $OpenBSD: e_sqrtl.c,v 1.3 2016/09/12 19:47:02 guenther Exp $ */ 2/*- 3 * Copyright (c) 2007 Steven G. Kargl 4 * All rights reserved. 5 * 6 * Redistribution and use in source and binary forms, with or without 7 * modification, are permitted provided that the following conditions 8 * are met: 9 * 1. Redistributions of source code must retain the above copyright 10 * notice unmodified, this list of conditions, and the following 11 * disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 16 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 17 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 18 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 19 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 20 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 21 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 22 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 23 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 25 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 */ 27 28#include <sys/types.h> 29#include <machine/ieee.h> /* for struct ieee_ext */ 30#include <fenv.h> 31#include <float.h> 32#include <math.h> 33 34#ifdef EXT_IMPLICIT_NBIT 35#define LDBL_NBIT 0 36#else /* EXT_IMPLICIT_NBIT */ 37#define LDBL_NBIT 0x80000000 38#endif /* EXT_IMPLICIT_NBIT */ 39 40/* Return (x + ulp) for normal positive x. Assumes no overflow. */ 41static inline long double 42inc(long double x) 43{ 44 struct ieee_ext *p = (struct ieee_ext *)&x; 45 46#ifdef EXT_FRACHMBITS 47 uint64_t frach; 48 49 frach = ((uint64_t)p->ext_frach << EXT_FRACHMBITS) | p->ext_frachm; 50 frach++; 51 p->ext_frach = frach >> EXT_FRACHMBITS; 52 p->ext_frachm = frach & 0x00000000ffffffff; 53#else /* EXT_FRACHMBITS */ 54 uint32_t frach; 55 56 p->ext_frach++; 57 frach = p->ext_frach; 58#endif /* EXT_FRACHMBITS */ 59 60 if (frach == 0) { 61 62#ifdef EXT_FRACLMBITS 63 uint64_t fracl; 64 65 fracl = ((uint64_t)p->ext_fraclm << EXT_FRACLBITS) | 66 p->ext_fracl; 67 fracl++; 68 p->ext_fraclm = fracl >> EXT_FRACLBITS; 69 p->ext_fracl = fracl & 0x00000000ffffffff; 70#else /* EXT_FRACLMBITS */ 71 uint32_t fracl; 72 73 p->ext_fracl++; 74 fracl = p->ext_fracl; 75#endif /* EXT_FRACLMBITS */ 76 77 if (fracl == 0) { 78 p->ext_exp++; 79 p->ext_frach |= LDBL_NBIT; 80 } 81 } 82 83 return x; 84} 85 86/* Return (x - ulp) for normal positive x. Assumes no underflow. */ 87static inline long double 88dec(long double x) 89{ 90 struct ieee_ext *p = (struct ieee_ext *)&x; 91 92#ifdef EXT_FRACLMBITS 93 uint64_t fracl; 94 95 fracl = ((uint64_t)p->ext_fraclm << EXT_FRACLBITS) | p->ext_fracl; 96 fracl--; 97 p->ext_fraclm = fracl >> EXT_FRACLBITS; 98 p->ext_fracl = fracl & 0x00000000ffffffff; 99#else /* EXT_FRACLMBITS */ 100 uint32_t fracl; 101 102 p->ext_fracl--; 103 fracl = p->ext_fracl; 104#endif /* EXT_FRACLMBITS */ 105 106 if (fracl == 0) { 107 108#ifdef EXT_FRACHMBITS 109 uint64_t frach; 110 111 frach = ((uint64_t)p->ext_frach << EXT_FRACHMBITS) | 112 p->ext_frachm; 113 frach--; 114 p->ext_frach = frach >> EXT_FRACHMBITS; 115 p->ext_frachm = frach & 0x00000000ffffffff; 116#else /* EXT_FRACHMBITS */ 117 uint32_t frach; 118 119 p->ext_frach--; 120 frach = p->ext_frach; 121#endif /* EXT_FRACHMBITS */ 122 123 if (frach == LDBL_NBIT) { 124 p->ext_exp--; 125 p->ext_frach |= LDBL_NBIT; 126 } 127 } 128 129 return x; 130} 131 132/* 133 * This is slow, but simple and portable. You should use hardware sqrt 134 * if possible. 135 */ 136 137long double 138sqrtl(long double x) 139{ 140 union { 141 long double e; 142 struct ieee_ext bits; 143 } u; 144 int k, r; 145 long double lo, xn; 146 147 u.e = x; 148 149 /* If x = NaN, then sqrt(x) = NaN. */ 150 /* If x = Inf, then sqrt(x) = Inf. */ 151 /* If x = -Inf, then sqrt(x) = NaN. */ 152 if (u.bits.ext_exp == LDBL_MAX_EXP * 2 - 1) 153 return (x * x + x); 154 155 /* If x = +-0, then sqrt(x) = +-0. */ 156 if ((u.bits.ext_frach 157#ifdef EXT_FRACHMBITS 158 | u.bits.ext_frachm 159#endif /* EXT_FRACHMBITS */ 160#ifdef EXT_FRACLMBITS 161 | u.bits.ext_fraclm 162#endif /* EXT_FRACLMBITS */ 163 | u.bits.ext_fracl | u.bits.ext_exp) == 0) 164 return (x); 165 166 /* If x < 0, then raise invalid and return NaN */ 167 if (u.bits.ext_sign) 168 return ((x - x) / (x - x)); 169 170 if (u.bits.ext_exp == 0) { 171 /* Adjust subnormal numbers. */ 172 u.e *= 0x1.0p514; 173 k = -514; 174 } else { 175 k = 0; 176 } 177 /* 178 * u.e is a normal number, so break it into u.e = e*2^n where 179 * u.e = (2*e)*2^2k for odd n and u.e = (4*e)*2^2k for even n. 180 */ 181 if ((u.bits.ext_exp - 0x3ffe) & 1) { /* n is odd. */ 182 k += u.bits.ext_exp - 0x3fff; /* 2k = n - 1. */ 183 u.bits.ext_exp = 0x3fff; /* u.e in [1,2). */ 184 } else { 185 k += u.bits.ext_exp - 0x4000; /* 2k = n - 2. */ 186 u.bits.ext_exp = 0x4000; /* u.e in [2,4). */ 187 } 188 189 /* 190 * Newton's iteration. 191 * Split u.e into a high and low part to achieve additional precision. 192 */ 193 xn = sqrt(u.e); /* 53-bit estimate of sqrtl(x). */ 194#if LDBL_MANT_DIG > 100 195 xn = (xn + (u.e / xn)) * 0.5; /* 106-bit estimate. */ 196#endif 197 lo = u.e; 198 u.bits.ext_fracl = 0; /* Zero out lower bits. */ 199#ifdef EXT_FRACLMBITS 200 u.bits.ext_fraclm = 0; 201#endif /* EXT_FRACLMBITS */ 202 lo = (lo - u.e) / xn; /* Low bits divided by xn. */ 203 xn = xn + (u.e / xn); /* High portion of estimate. */ 204 u.e = xn + lo; /* Combine everything. */ 205 u.bits.ext_exp += (k >> 1) - 1; 206 207 feclearexcept(FE_INEXACT); 208 r = fegetround(); 209 fesetround(FE_TOWARDZERO); /* Set to round-toward-zero. */ 210 xn = x / u.e; /* Chopped quotient (inexact?). */ 211 212 if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */ 213 if (xn == u.e) { 214 fesetround(r); 215 return (u.e); 216 } 217 /* Round correctly for inputs like x = y**2 - ulp. */ 218 xn = dec(xn); /* xn = xn - ulp. */ 219 } 220 221 if (r == FE_TONEAREST) { 222 xn = inc(xn); /* xn = xn + ulp. */ 223 } else if (r == FE_UPWARD) { 224 u.e = inc(u.e); /* u.e = u.e + ulp. */ 225 xn = inc(xn); /* xn = xn + ulp. */ 226 } 227 u.e = u.e + xn; /* Chopped sum. */ 228 fesetround(r); /* Restore env and raise inexact */ 229 u.bits.ext_exp--; 230 return (u.e); 231} 232DEF_STD(sqrtl); 233