1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 /* 22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #include "libm.h" /* __k_clog_r */ 30 #include "complex_wrapper.h" 31 32 /* INDENT OFF */ 33 /* 34 * double __k_clog_r(double x, double y, double *e); 35 * 36 * Compute real part of complex natural logarithm of x+iy in extra precision 37 * 38 * __k_clog_r returns log(hypot(x, y)) with a correction term e. 39 * 40 * Accuracy: 70 bits 41 * 42 * Method. 43 * Let Z = x*x + y*y. Z can be normalized as Z = 2^N * z, 1 <= z < 2. 44 * We further break down z into 1 + zk + zh + zt, where 45 * zk = K*(2^-7) matches z to 7.5 significant bits, 0 <= K <= 2^(-7)-1 46 * zh = (z-zk) rounded to 24 bits 47 * zt = (z-zk-zh) rounded. 48 * 49 * z - (1+zk) (zh+zt) 50 * Let s = ------------ = ---------------, then 51 * z + (1+zk) 2(1+zk)+zh+zt 52 * z 53 * log(Z) = N*log2 + log(z) = N*log2 + log(1+zk) + log(------) 54 * 1+zk 55 * 1+s 56 * = N*log2 + log(1+zk) + log(---) 57 * 1-s 58 * 59 * 1 3 1 5 60 * = N*log2 + log(1+zk) + 2s + -- (2s) + -- (2s) + ... 61 * 12 80 62 * 63 * Note 1. For IEEE double precision, a seven degree odd polynomial 64 * 2s + P1*(2s)^3 + P2*(2s)^5 + P3*(2s)^7 65 * is generated by a special remez algorithm to 66 * approx log((1+s)/(1-s)) accurte to 72 bits. 67 * Note 2. 2s can be computed accurately as s2h+s2t by 68 * r = 2/((zh+zt)+2(1+zk)) 69 * s2 = r*(zh+zt) 70 * s2h = s2 rounded to float; v = 0.5*s2h; 71 * s2t = r*((((zh-s2h*(1+zk))-v*zh)+zt)-v*zt) 72 */ 73 /* INDENT ON */ 74 75 static const double 76 zero = 0.0, 77 half = 0.5, 78 two = 2.0, 79 two120 = 1.32922799578491587290e+36, /* 2^120 */ 80 ln2_h = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 81 ln2_t = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 82 P1 = .083333333333333351554108717377986202224765262191125, 83 P2 = .01249999999819227552330700574633767185896464873834375, 84 P3 = .0022321938458645656605471559987512516234702284287265625; 85 86 /* 87 * T[2k, 2k+1] = log(1+k*2^-7) for k = 0, ..., 2^7 - 1, 88 * with T[2k] * 2^40 is an int 89 */ 90 91 static const double TBL_log1k[] = { 92 0.00000000000000000000e+00, 0.00000000000000000000e+00, 93 7.78214044203195953742e-03, 2.29894100462035112076e-14, 94 1.55041865355087793432e-02, 4.56474807636434698847e-13, 95 2.31670592811497044750e-02, 3.84673753843363762372e-13, 96 3.07716586667083902285e-02, 4.52981425779092882775e-14, 97 3.83188643018002039753e-02, 3.36395218465265063278e-13, 98 4.58095360309016541578e-02, 3.92549008891706208826e-13, 99 5.32445145181554835290e-02, 6.56799336898521766515e-13, 100 6.06246218158048577607e-02, 6.29984819938331143924e-13, 101 6.79506619080711971037e-02, 4.36552290856295281946e-13, 102 7.52234212368421140127e-02, 7.45411685916941618656e-13, 103 8.24436692109884461388e-02, 8.61451293608781447223e-14, 104 8.96121586893059429713e-02, 3.81189648692113819551e-13, 105 9.67296264579999842681e-02, 5.51128027471986918274e-13, 106 1.03796793680885457434e-01, 7.58107392301637643358e-13, 107 1.10814366339582193177e-01, 7.07921017612766061755e-13, 108 1.17783035655520507134e-01, 8.62947404296943765415e-13, 109 1.24703478500123310369e-01, 8.33925494898414856118e-13, 110 1.31576357788617315236e-01, 1.01957352237084734958e-13, 111 1.38402322858382831328e-01, 7.36304357708705134617e-13, 112 1.45182009843665582594e-01, 8.32314688404647202319e-13, 113 1.51916042025732167531e-01, 1.09807540998552379211e-13, 114 1.58605030175749561749e-01, 8.89022343972466269900e-13, 115 1.65249572894936136436e-01, 3.71026439894104998399e-13, 116 1.71850256926518341061e-01, 1.40881279371111350341e-13, 117 1.78407657472234859597e-01, 5.83437522462346671423e-13, 118 1.84922338493379356805e-01, 6.32635858668445232946e-13, 119 1.91394852999110298697e-01, 5.19155912393432989209e-13, 120 1.97825743329303804785e-01, 6.16075577558872326221e-13, 121 2.04215541428311553318e-01, 3.79338185766902218086e-13, 122 2.10564769106895255391e-01, 4.54382278998146218219e-13, 123 2.16873938300523150247e-01, 9.12093724991498410553e-14, 124 2.23143551314024080057e-01, 1.85675709597960106615e-13, 125 2.29374101064422575291e-01, 4.23254700234549300166e-13, 126 2.35566071311950508971e-01, 8.16400106820959292914e-13, 127 2.41719936886511277407e-01, 6.33890736899755317832e-13, 128 2.47836163904139539227e-01, 4.41717553713155466566e-13, 129 2.53915209980732470285e-01, 2.30973852175869394892e-13, 130 2.59957524436686071567e-01, 2.39995404842117353465e-13, 131 2.65963548496984003577e-01, 1.53937761744554075681e-13, 132 2.71933715483100968413e-01, 5.40790418614551497411e-13, 133 2.77868451003087102436e-01, 3.69203750820800887027e-13, 134 2.83768173129828937817e-01, 8.15660529536291275782e-13, 135 2.89633292582948342897e-01, 9.43339818951269030846e-14, 136 2.95464212893421063200e-01, 4.14813187042585679830e-13, 137 3.01261330577290209476e-01, 8.71571536970835103739e-13, 138 3.07025035294827830512e-01, 8.40315630479242455758e-14, 139 3.12755710003330023028e-01, 5.66865358290073900922e-13, 140 3.18453731118097493891e-01, 4.37121919574291444278e-13, 141 3.24119468653407238889e-01, 8.04737201185162774515e-13, 142 3.29753286371669673827e-01, 7.98307987877335024112e-13, 143 3.35355541920762334485e-01, 3.75495772572598557174e-13, 144 3.40926586970454081893e-01, 1.39128412121975659358e-13, 145 3.46466767346100823488e-01, 1.07757430375726404546e-13, 146 3.51976423156884266064e-01, 2.93918591876480007730e-13, 147 3.57455888921322184615e-01, 4.81589611172320539489e-13, 148 3.62905493689140712377e-01, 2.27740761140395561986e-13, 149 3.68325561158599157352e-01, 1.08495696229679121506e-13, 150 3.73716409792905324139e-01, 6.78756682315870616582e-13, 151 3.79078352934811846353e-01, 1.57612037739694350287e-13, 152 3.84411698910298582632e-01, 3.34571026954408237380e-14, 153 3.89716751139530970249e-01, 4.94243121138567024911e-13, 154 3.94993808240542421117e-01, 3.26556988969071456956e-13, 155 4.00243164126550254878e-01, 4.62452051668403792833e-13, 156 4.05465108107819105498e-01, 3.45276479520397708744e-13, 157 4.10659924984429380856e-01, 8.39005077851830734139e-13, 158 4.15827895143593195826e-01, 1.17769787513692141889e-13, 159 4.20969294643327884842e-01, 8.01751287156832458079e-13, 160 4.26084395310681429692e-01, 2.18633432932159103190e-13, 161 4.31173464818130014464e-01, 2.41326394913331314894e-13, 162 4.36236766774527495727e-01, 3.90574622098307022265e-13, 163 4.41274560804231441580e-01, 6.43787909737320689684e-13, 164 4.46287102628048160113e-01, 3.71351419195920213229e-13, 165 4.51274644138720759656e-01, 7.37825488412103968058e-13, 166 4.56237433480964682531e-01, 6.22911850193784704748e-13, 167 4.61175715121498797089e-01, 6.71369279138460114513e-13, 168 4.66089729924533457961e-01, 6.57665976858006147528e-14, 169 4.70979715218163619284e-01, 6.27393263311115598424e-13, 170 4.75845904869856894948e-01, 1.07019317621142549209e-13, 171 4.80688529345570714213e-01, 1.81193463664411114729e-13, 172 4.85507815781602403149e-01, 9.84046527823262695501e-14, 173 4.90303988044615834951e-01, 5.78003198945402769376e-13, 174 4.95077266797125048470e-01, 7.26466128212511528295e-13, 175 4.99827869555701909121e-01, 7.47420700205478712293e-13, 176 5.04556010751912253909e-01, 4.83033149495532022300e-13, 177 5.09261901789614057634e-01, 1.93889170049107088943e-13, 178 5.13945751101346104406e-01, 8.88212395185718544720e-13, 179 5.18607764207445143256e-01, 6.00488896640545761201e-13, 180 5.23248143764249107335e-01, 2.98729182044413286731e-13, 181 5.27867089620485785417e-01, 3.56599696633478298092e-13, 182 5.32464798869114019908e-01, 3.57823965912763837621e-13, 183 5.37041465896436420735e-01, 4.47233831757482468946e-13, 184 5.41597282432121573947e-01, 6.22797629172251525649e-13, 185 5.46132437597407260910e-01, 7.28389472720657362987e-13, 186 5.50647117952394182794e-01, 2.68096466152116723636e-13, 187 5.55141507539701706264e-01, 7.99886451312335479470e-13, 188 5.59615787935399566777e-01, 2.31194938380053776320e-14, 189 5.64070138284478161950e-01, 3.24804121719935740729e-13, 190 5.68504735351780254859e-01, 8.88457219261483317716e-13, 191 5.72919753561109246220e-01, 6.76262872317054154667e-13, 192 5.77315365034337446559e-01, 4.86157758891509033842e-13, 193 5.81691739634152327199e-01, 4.70155322075549811780e-13, 194 5.86049045003164792433e-01, 4.13416470738355643357e-13, 195 5.90387446602107957006e-01, 6.84176364159146659095e-14, 196 5.94707107746216934174e-01, 4.75855340044306376333e-13, 197 5.99008189645246602595e-01, 8.36796786747576938145e-13, 198 6.03290851438032404985e-01, 5.18573553063418286042e-14, 199 6.07555250224322662689e-01, 2.19132812293400917731e-13, 200 6.11801541105705837253e-01, 2.87066276408616768331e-13, 201 6.16029877214714360889e-01, 7.99658758518543977451e-13, 202 6.20240409751204424538e-01, 6.53104313776336534177e-13, 203 6.24433288011459808331e-01, 4.33692711555820529733e-13, 204 6.28608659421843185555e-01, 5.30952189118357790115e-13, 205 6.32766669570628437214e-01, 4.09392332186786656392e-13, 206 6.36907462236194987781e-01, 8.74243839148582888557e-13, 207 6.41031179420679109171e-01, 2.52181884568428814231e-13, 208 6.45137961372711288277e-01, 8.73413388168702670246e-13, 209 6.49227946624705509748e-01, 4.04309142530119209805e-13, 210 6.53301272011958644725e-01, 7.86994033233553225797e-13, 211 6.57358072708120744210e-01, 2.39285932153437645135e-13, 212 6.61398482245203922503e-01, 1.61085757539324585156e-13, 213 6.65422632544505177066e-01, 5.85271884362515112697e-13, 214 6.69430653942072240170e-01, 5.57027128793880294600e-13, 215 6.73422675211440946441e-01, 7.25773856816637653180e-13, 216 6.77398823590920073912e-01, 8.86066898134949155668e-13, 217 6.81359224807238206267e-01, 6.64862680714687006264e-13, 218 6.85304003098281100392e-01, 6.38316151706465171657e-13, 219 6.89233281238557538018e-01, 2.51442307283760746611e-13, 220 }; 221 222 /* 223 * Compute N*log2 + log(1+zk+zh+zt) in extra precision 224 */ 225 static double k_log_NKz(int N, int K, double zh, double *zt) 226 { 227 double y, r, w, s2, s2h, s2t, t, zk, v, P; 228 229 ((int *)&zk)[HIWORD] = 0x3ff00000 + (K << 13); 230 ((int *)&zk)[LOWORD] = 0; 231 t = zh + (*zt); 232 r = two / (t + two * zk); 233 s2h = s2 = r * t; 234 ((int *)&s2h)[LOWORD] &= 0xe0000000; 235 v = half * s2h; 236 w = s2 * s2; 237 s2t = r * ((((zh - s2h * zk) - v * zh) + (*zt)) - v * (*zt)); 238 P = s2t + (w * s2) * ((P1 + w * P2) + (w * w) * P3); 239 P += N * ln2_t + TBL_log1k[K + K + 1]; 240 t = N*ln2_h + TBL_log1k[K+K]; 241 y = t + (P + s2h); 242 P -= ((y - t) - s2h); 243 *zt = P; 244 return (y); 245 } 246 247 double 248 __k_clog_r(double x, double y, double *er) 249 { 250 double t1, t2, t3, t4, tk, z, wh, w, zh, zk; 251 int n, k, ix, iy, iz, nx, ny, nz, i, j; 252 unsigned lx, ly; 253 254 ix = (((int *)&x)[HIWORD]) & ~0x80000000; 255 lx = ((unsigned *)&x)[LOWORD]; 256 iy = (((int *)&y)[HIWORD]) & ~0x80000000; 257 ly = ((unsigned *)&y)[LOWORD]; 258 y = fabs(y); x = fabs(x); 259 if (ix < iy || (ix == iy && lx < ly)) { /* force x >= y */ 260 tk = x; x = y; y = tk; 261 n = ix, ix = iy; iy = n; 262 n = lx, lx = ly; ly = n; 263 } 264 *er = zero; 265 nx = ix >> 20; ny = iy >> 20; 266 if (nx >= 0x7ff) { /* x or y is Inf or NaN */ 267 if (ISINF(ix, lx)) 268 return (x); 269 else if (ISINF(iy, ly)) 270 return (y); 271 else 272 return (x+y); 273 } 274 /* 275 * for tiny y (double y < 2^-35, extended y < 2^-46, quad y < 2^-70): 276 * log(sqrt(1+y^2)) = (y^2)/2 - (y^4)/8 + ... ~= (y^2)/2 277 */ 278 if ((((ix - 0x3ff00000) | lx) == 0) && ny < (0x3ff - 35)) { 279 t2 = y * y; 280 if (ny >= 565) { /* compute er = tail of t2 */ 281 ((int *)&wh)[HIWORD] = iy; 282 ((unsigned *)&wh)[LOWORD] = ly & 0xf8000000; 283 *er = half * ((y - wh) * (y + wh) - (t2 - wh * wh)); 284 } 285 return (half * t2); 286 } 287 /* 288 * x or y is subnormal or zero 289 */ 290 if (nx == 0) { 291 if ((ix | lx) == 0) 292 return (-1.0 / x); 293 else { 294 x *= two120; 295 y *= two120; 296 ix = ((int *)&x)[HIWORD]; 297 lx = ((unsigned *)&x)[LOWORD]; 298 iy = ((int *)&y)[HIWORD]; 299 ly = ((unsigned *)&y)[LOWORD]; 300 nx = (ix >> 20) - 120; 301 ny = (iy >> 20) - 120; 302 /* guard subnormal flush to 0 */ 303 if ((ix | lx) == 0) 304 return (-1.0 / x); 305 } 306 } else if (ny == 0) { /* y subnormal, scale it */ 307 y *= two120; 308 iy = ((int *)&y)[HIWORD]; 309 ly = ((unsigned *)&y)[LOWORD]; 310 ny = (iy >> 20) - 120; 311 } 312 n = nx - ny; 313 /* 314 * return log(x) when y is zero or x >> y so that 315 * log(x) ~ log(sqrt(x*x+y*y)) to 27 extra bits 316 * (n > 62 for double, 78 for i386 extended, 122 for quad) 317 */ 318 if (n > 62 || (iy | ly) == 0) { 319 i = (0x000fffff & ix) | 0x3ff00000; /* normalize x */ 320 ((int *)&x)[HIWORD] = i; 321 i += 0x1000; 322 ((int *)&zk)[HIWORD] = i & 0xffffe000; 323 ((int *)&zk)[LOWORD] = 0; /* zk matches 7.5 bits of x */ 324 z = x - zk; 325 zh = (double)((float)z); 326 i >>= 13; 327 k = i & 0x7f; /* index of zk */ 328 n = nx - 0x3ff; 329 *er = z - zh; 330 if (i >> 17) { /* if zk = 2.0, adjust scaling */ 331 n += 1; 332 zh *= 0.5; *er *= 0.5; 333 } 334 w = k_log_NKz(n, k, zh, er); 335 } else { 336 /* 337 * compute z = x*x + y*y 338 */ 339 ix = (ix & 0xfffff) | 0x3ff00000; 340 iy = (iy & 0xfffff) | (0x3ff00000 - (n << 20)); 341 ((int *)&x)[HIWORD] = ix; ((int *)&y)[HIWORD] = iy; 342 t1 = x * x; t2 = y * y; 343 j = ((lx >> 26) + 1) >> 1; 344 ((int *)&wh)[HIWORD] = ix + (j >> 5); 345 ((unsigned *)&wh)[LOWORD] = (j << 27); 346 z = t1+t2; 347 /* 348 * higher precision simulation x*x = t1 + t3, y*y = t2 + t4 349 */ 350 tk = wh - x; 351 t3 = tk * tk - (two * wh * tk - (wh * wh - t1)); 352 j = ((ly >> 26) + 1) >> 1; 353 ((int *)&wh)[HIWORD] = iy + (j >> 5); 354 ((unsigned *)&wh)[LOWORD] = (j << 27); 355 tk = wh - y; 356 t4 = tk * tk - (two * wh * tk - (wh * wh - t2)); 357 /* 358 * find zk matches z to 7.5 bits 359 */ 360 nx -= 0x3ff; 361 iz = ((int *)&z)[HIWORD] + 0x1000; 362 k = (iz >> 13) & 0x7f; 363 nz = (iz >> 20) - 0x3ff; 364 ((int *)&zk)[HIWORD] = iz & 0xffffe000; 365 ((int *)&zk)[LOWORD] = 0; 366 /* 367 * order t1,t2,t3,t4 according to their size 368 */ 369 if (t2 >= fabs(t3)) { 370 if (fabs(t3) < fabs(t4)) { 371 wh = t3; t3 = t4; t4 = wh; 372 } 373 } else { 374 wh = t2; t2 = t3; t3 = wh; 375 } 376 /* 377 * higher precision simulation: x * x + y * y = t1 + t2 + t3 + t4 378 * = zk (7 bits) + zh (24 bits) + *er (tail) and call k_log_NKz 379 */ 380 tk = t1 - zk; 381 zh = ((tk + t2) + t3) + t4; 382 ((int *)&zh)[LOWORD] &= 0xe0000000; 383 w = fabs(zh); 384 if (w >= fabs(t2)) 385 *er = (((tk - zh) + t2) + t3) + t4; 386 else { 387 if (n == 0) { 388 wh = half * zk; 389 wh = (t1 - wh) - (wh - t2); 390 } else 391 wh = tk + t2; 392 if (w >= fabs(t3)) 393 *er = ((wh - zh) + t3) + t4; 394 else { 395 z = t3; 396 t3 += t4; 397 t4 -= t3 - z; 398 if (w >= fabs(t3)) 399 *er = ((wh - zh) + t3) + t4; 400 else 401 *er = ((wh + t3) - zh) + t4; 402 } 403 } 404 if (nz == 3) {zh *= 0.125; *er *= 0.125; } 405 if (nz == 2) {zh *= 0.25; *er *= 0.25; } 406 if (nz == 1) {zh *= half; *er *= half; } 407 nz += nx + nx; 408 w = half * k_log_NKz(nz, k, zh, er); 409 *er *= half; 410 } 411 return (w); 412 } 413