1*25c28e83SPiotr Jasiukajtis /* 2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START 3*25c28e83SPiotr Jasiukajtis * 4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 7*25c28e83SPiotr Jasiukajtis * 8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 11*25c28e83SPiotr Jasiukajtis * and limitations under the License. 12*25c28e83SPiotr Jasiukajtis * 13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 18*25c28e83SPiotr Jasiukajtis * 19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END 20*25c28e83SPiotr Jasiukajtis */ 21*25c28e83SPiotr Jasiukajtis 22*25c28e83SPiotr Jasiukajtis /* 23*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24*25c28e83SPiotr Jasiukajtis */ 25*25c28e83SPiotr Jasiukajtis /* 26*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27*25c28e83SPiotr Jasiukajtis * Use is subject to license terms. 28*25c28e83SPiotr Jasiukajtis */ 29*25c28e83SPiotr Jasiukajtis 30*25c28e83SPiotr Jasiukajtis #include <sys/isa_defs.h> 31*25c28e83SPiotr Jasiukajtis 32*25c28e83SPiotr Jasiukajtis #ifdef _LITTLE_ENDIAN 33*25c28e83SPiotr Jasiukajtis #define HI(x) *(1+(int*)x) 34*25c28e83SPiotr Jasiukajtis #define LO(x) *(unsigned*)x 35*25c28e83SPiotr Jasiukajtis #else 36*25c28e83SPiotr Jasiukajtis #define HI(x) *(int*)x 37*25c28e83SPiotr Jasiukajtis #define LO(x) *(1+(unsigned*)x) 38*25c28e83SPiotr Jasiukajtis #endif 39*25c28e83SPiotr Jasiukajtis 40*25c28e83SPiotr Jasiukajtis #ifdef __RESTRICT 41*25c28e83SPiotr Jasiukajtis #define restrict _Restrict 42*25c28e83SPiotr Jasiukajtis #else 43*25c28e83SPiotr Jasiukajtis #define restrict 44*25c28e83SPiotr Jasiukajtis #endif 45*25c28e83SPiotr Jasiukajtis 46*25c28e83SPiotr Jasiukajtis /* double pow(double x, double y) 47*25c28e83SPiotr Jasiukajtis * 48*25c28e83SPiotr Jasiukajtis * Method : 49*25c28e83SPiotr Jasiukajtis * 1. Special cases: 50*25c28e83SPiotr Jasiukajtis * for (anything) ** 0 => 1 51*25c28e83SPiotr Jasiukajtis * for (anything) ** NaN => QNaN + invalid 52*25c28e83SPiotr Jasiukajtis * for NaN ** (anything) => QNaN + invalid 53*25c28e83SPiotr Jasiukajtis * for +-1 ** +-Inf => QNaN + invalid 54*25c28e83SPiotr Jasiukajtis * for +-(|x| < 1) ** +Inf => +0 55*25c28e83SPiotr Jasiukajtis * for +-(|x| < 1) ** -Inf => +Inf 56*25c28e83SPiotr Jasiukajtis * for +-(|x| > 1) ** +Inf => +Inf 57*25c28e83SPiotr Jasiukajtis * for +-(|x| > 1) ** -Inf => +0 58*25c28e83SPiotr Jasiukajtis * for +Inf ** (negative) => +0 59*25c28e83SPiotr Jasiukajtis * for +Inf ** (positive) => +Inf 60*25c28e83SPiotr Jasiukajtis * for -Inf ** (negative except odd integer) => +0 61*25c28e83SPiotr Jasiukajtis * for -Inf ** (negative odd integer) => -0 62*25c28e83SPiotr Jasiukajtis * for -Inf ** (positive except odd integer) => +Inf 63*25c28e83SPiotr Jasiukajtis * for -Inf ** (positive odd integer) => -Inf 64*25c28e83SPiotr Jasiukajtis * for (negative) ** (non-integer) => QNaN + invalid 65*25c28e83SPiotr Jasiukajtis * for +0 ** (negative) => +Inf + overflow 66*25c28e83SPiotr Jasiukajtis * for +0 ** (positive) => +0 67*25c28e83SPiotr Jasiukajtis * for -0 ** (negative except odd integer) => +Inf + overflow 68*25c28e83SPiotr Jasiukajtis * for -0 ** (negative odd integer) => -Inf + overflow 69*25c28e83SPiotr Jasiukajtis * for -0 ** (positive except odd integer) => +0 70*25c28e83SPiotr Jasiukajtis * for -0 ** (positive odd integer) => -0 71*25c28e83SPiotr Jasiukajtis * 2. Computes x**y from: 72*25c28e83SPiotr Jasiukajtis * x**y = 2**(y*log2(x)) = 2**(w/256), where w = 256*log2(|x|)*y. 73*25c28e83SPiotr Jasiukajtis * 3. Computes w = 256*log2(|x|)*y from 74*25c28e83SPiotr Jasiukajtis * |x| = m * 2**n => log2(|x|) = n + log2(m). 75*25c28e83SPiotr Jasiukajtis * Let m = m0 + dm, where m0 = 1 + k / 256, 76*25c28e83SPiotr Jasiukajtis * k = [0, 255], 77*25c28e83SPiotr Jasiukajtis * dm = [-1/512, 1/512]. 78*25c28e83SPiotr Jasiukajtis * Then 256*log2(m) = 256*log2(m0 + dm) = 256*log2(m0) + 256*log2((1+z)/(1-z)), 79*25c28e83SPiotr Jasiukajtis * where z = (m-m0)/(m+m0), z = [-1/1025, 1/1025]. 80*25c28e83SPiotr Jasiukajtis * Then 81*25c28e83SPiotr Jasiukajtis * 256*log2(m0) is looked up in a table of 256*log2(1), 256*log2(1+1/128), 82*25c28e83SPiotr Jasiukajtis * ..., 256*log2(1+128/128). 83*25c28e83SPiotr Jasiukajtis * 256*log2((1+z)/(1-z)) is computed using 84*25c28e83SPiotr Jasiukajtis * approximation: 256*log2((1+z)/(1-z)) = a0 * z + a1 * z**3 + a1 * z**5. 85*25c28e83SPiotr Jasiukajtis * Perform w = 256*log2(|x|)*y = w1 + w2 by simulating muti-precision arithmetic. 86*25c28e83SPiotr Jasiukajtis * 3. For w >= 262144 87*25c28e83SPiotr Jasiukajtis * then for (negative) ** (odd integer) => -Inf + overflow 88*25c28e83SPiotr Jasiukajtis * else => +Inf + overflow 89*25c28e83SPiotr Jasiukajtis * For w <= -275200 90*25c28e83SPiotr Jasiukajtis * then for (negative) ** (odd integer) => -0 + underflow 91*25c28e83SPiotr Jasiukajtis * else => +0 + underflow 92*25c28e83SPiotr Jasiukajtis * 4. Computes 2 ** (w/256) from: 93*25c28e83SPiotr Jasiukajtis * 2 ** (w/256) = 2**a * 2**(k/256) * 2**(r/256) 94*25c28e83SPiotr Jasiukajtis * Where: 95*25c28e83SPiotr Jasiukajtis * a = int ( w ) >> 8; 96*25c28e83SPiotr Jasiukajtis * k = int ( w ) & 0xFF; 97*25c28e83SPiotr Jasiukajtis * r = frac ( w ). 98*25c28e83SPiotr Jasiukajtis * Note that: 99*25c28e83SPiotr Jasiukajtis * k = 0, 1, ..., 255; 100*25c28e83SPiotr Jasiukajtis * r = (-1, 1). 101*25c28e83SPiotr Jasiukajtis * Then: 102*25c28e83SPiotr Jasiukajtis * 2**(k/256) is looked up in a table of 2**0, 2**1/256, ... 103*25c28e83SPiotr Jasiukajtis * 2**(r/256) is computed using approximation: 104*25c28e83SPiotr Jasiukajtis * 2**(r/256) = ((((b5 * r + b4) * r + b3) * r + b2) * r + b1) * r + b0 105*25c28e83SPiotr Jasiukajtis * Multiplication by 2**a is done by adding "a" to 106*25c28e83SPiotr Jasiukajtis * the biased exponent. 107*25c28e83SPiotr Jasiukajtis * Perform 2 ** (w/256) by simulating muti-precision arithmetic. 108*25c28e83SPiotr Jasiukajtis * 5. For (negative) ** (odd integer) => -(2**(w/256)) 109*25c28e83SPiotr Jasiukajtis * otherwise => 2**(w/256) 110*25c28e83SPiotr Jasiukajtis * 111*25c28e83SPiotr Jasiukajtis * Accuracy: 112*25c28e83SPiotr Jasiukajtis * Max. relative aproximation error < 2**(-67.94) for 256*log2((1+z)/(1-z)). 113*25c28e83SPiotr Jasiukajtis * Max. relative aproximation error < 2**(-63.15) for 2**(r/256). 114*25c28e83SPiotr Jasiukajtis * Maximum error observed: less than 0.761 ulp after 1.300.000.000 115*25c28e83SPiotr Jasiukajtis * results. 116*25c28e83SPiotr Jasiukajtis */ 117*25c28e83SPiotr Jasiukajtis 118*25c28e83SPiotr Jasiukajtis static void 119*25c28e83SPiotr Jasiukajtis __vpowx(int n, double * restrict px, double * restrict py, 120*25c28e83SPiotr Jasiukajtis int stridey, double * restrict pz, int stridez); 121*25c28e83SPiotr Jasiukajtis 122*25c28e83SPiotr Jasiukajtis static const double __TBL_exp2[] = { 123*25c28e83SPiotr Jasiukajtis /* __TBL_exp2[2*i] = high order bits 2^(i/256), i = [0, 255] */ 124*25c28e83SPiotr Jasiukajtis /* __TBL_exp2[2*i+1] = least bits 2^(i/256), i = [0, 255] */ 125*25c28e83SPiotr Jasiukajtis 1.000000000000000000e+00, 0.000000000000000000e+00, 1.002711275050202522e+00, 126*25c28e83SPiotr Jasiukajtis -3.636615928692263944e-17, 1.005429901112802726e+00, 9.499186535455031757e-17, 127*25c28e83SPiotr Jasiukajtis 1.008155898118417548e+00,-3.252058756084308061e-17, 1.010889286051700475e+00, 128*25c28e83SPiotr Jasiukajtis -1.523477860336857718e-17, 1.013630084951489430e+00, 9.283599768183567587e-18, 129*25c28e83SPiotr Jasiukajtis 1.016378314910953096e+00,-5.772170073199660028e-17, 1.019133996077737914e+00, 130*25c28e83SPiotr Jasiukajtis 3.601904982259661106e-17, 1.021897148654116627e+00, 5.109225028973443894e-17, 131*25c28e83SPiotr Jasiukajtis 1.024667792897135721e+00,-7.561607868487779440e-17, 1.027445949118763746e+00, 132*25c28e83SPiotr Jasiukajtis -4.956074174645370440e-17, 1.030231637686040980e+00, 3.319830041080812944e-17, 133*25c28e83SPiotr Jasiukajtis 1.033024879021228415e+00, 7.600838874027088489e-18, 1.035825693601957198e+00, 134*25c28e83SPiotr Jasiukajtis -7.806782391337636167e-17, 1.038634101961378731e+00, 5.996273788852510618e-17, 135*25c28e83SPiotr Jasiukajtis 1.041450124688316103e+00, 3.784830480287576210e-17, 1.044273782427413755e+00, 136*25c28e83SPiotr Jasiukajtis 8.551889705537964892e-17, 1.047105095879289793e+00, 7.277077243104314749e-17, 137*25c28e83SPiotr Jasiukajtis 1.049944085800687210e+00, 5.592937848127002586e-17, 1.052790773004626423e+00, 138*25c28e83SPiotr Jasiukajtis -9.629482899026935739e-17, 1.055645178360557157e+00, 1.759325738772091599e-18, 139*25c28e83SPiotr Jasiukajtis 1.058507322794512762e+00,-7.152651856637780738e-17, 1.061377227289262093e+00, 140*25c28e83SPiotr Jasiukajtis -1.197353708536565756e-17, 1.064254912884464499e+00, 5.078754198611230394e-17, 141*25c28e83SPiotr Jasiukajtis 1.067140400676823697e+00,-7.899853966841582122e-17, 1.070033711820241873e+00, 142*25c28e83SPiotr Jasiukajtis -9.937162711288919381e-17, 1.072934867525975555e+00,-3.839668843358823807e-18, 143*25c28e83SPiotr Jasiukajtis 1.075843889062791048e+00,-1.000271615114413611e-17, 1.078760797757119860e+00, 144*25c28e83SPiotr Jasiukajtis -6.656660436056592603e-17, 1.081685614993215250e+00,-4.782623902997086266e-17, 145*25c28e83SPiotr Jasiukajtis 1.084618362213309206e+00, 3.166152845816346116e-17, 1.087559060917769660e+00, 146*25c28e83SPiotr Jasiukajtis 5.409349307820290759e-18, 1.090507732665257690e+00,-3.046782079812471147e-17, 147*25c28e83SPiotr Jasiukajtis 1.093464399072885840e+00, 1.441395814726920934e-17, 1.096429081816376883e+00, 148*25c28e83SPiotr Jasiukajtis -5.919933484449315824e-17, 1.099401802630221914e+00, 7.170459599701923225e-17, 149*25c28e83SPiotr Jasiukajtis 1.102382583307840891e+00, 5.266036871570694387e-17, 1.105371445701741173e+00, 150*25c28e83SPiotr Jasiukajtis 8.239288760500213590e-17, 1.108368411723678726e+00,-8.786813845180526616e-17, 151*25c28e83SPiotr Jasiukajtis 1.111373503344817548e+00, 5.563945026669697643e-17, 1.114386742595892432e+00, 152*25c28e83SPiotr Jasiukajtis 1.041027845684557095e-16, 1.117408151567369279e+00,-7.976805902628220456e-17, 153*25c28e83SPiotr Jasiukajtis 1.120437752409606746e+00,-6.201085906554178750e-17, 1.123475567333019898e+00, 154*25c28e83SPiotr Jasiukajtis -9.699737588987042995e-17, 1.126521618608241848e+00, 5.165856758795456737e-17, 155*25c28e83SPiotr Jasiukajtis 1.129575928566288079e+00, 6.712805858726256588e-17, 1.132638519598719196e+00, 156*25c28e83SPiotr Jasiukajtis 3.237356166738000264e-17, 1.135709414157805464e+00, 5.066599926126155859e-17, 157*25c28e83SPiotr Jasiukajtis 1.138788634756691565e+00, 8.912812676025407778e-17, 1.141876203969561576e+00, 158*25c28e83SPiotr Jasiukajtis 4.651091177531412387e-17, 1.144972144431804173e+00, 4.641289892170010657e-17, 159*25c28e83SPiotr Jasiukajtis 1.148076478840178938e+00, 6.897740236627191770e-17, 1.151189229952982673e+00, 160*25c28e83SPiotr Jasiukajtis 3.250710218863827212e-17, 1.154310420590215935e+00, 1.041712894627326619e-16, 161*25c28e83SPiotr Jasiukajtis 1.157440073633751121e+00,-9.123871231134400287e-17, 1.160578212027498779e+00, 162*25c28e83SPiotr Jasiukajtis -3.261040205417393722e-17, 1.163724858777577476e+00, 3.829204836924093499e-17, 163*25c28e83SPiotr Jasiukajtis 1.166880036952481658e+00,-8.791879579999169742e-17, 1.170043769683250190e+00, 164*25c28e83SPiotr Jasiukajtis -1.847744201790004694e-18, 1.173216080163637320e+00,-7.287562586584994479e-17, 165*25c28e83SPiotr Jasiukajtis 1.176396991650281221e+00, 5.554203254218078963e-17, 1.179586527462875845e+00, 166*25c28e83SPiotr Jasiukajtis 1.009231277510039044e-16, 1.182784710984341014e+00, 1.542975430079076058e-17, 167*25c28e83SPiotr Jasiukajtis 1.185991565660993841e+00,-9.209506835293105905e-18, 1.189207115002721027e+00, 168*25c28e83SPiotr Jasiukajtis 3.982015231465646111e-17, 1.192431382583151178e+00, 4.397551415609721443e-17, 169*25c28e83SPiotr Jasiukajtis 1.195664392039827328e+00, 4.616603670481481397e-17, 1.198906167074380580e+00, 170*25c28e83SPiotr Jasiukajtis -9.809193356008423118e-17, 1.202156731452703076e+00, 6.644981499252301245e-17, 171*25c28e83SPiotr Jasiukajtis 1.205416109005123859e+00,-3.357272193267529634e-17, 1.208684323626581625e+00, 172*25c28e83SPiotr Jasiukajtis -4.746725945228984097e-17, 1.211961399276801243e+00,-4.890611077521118357e-17, 173*25c28e83SPiotr Jasiukajtis 1.215247359980468955e+00,-7.712630692681488131e-17, 1.218542229827408452e+00, 174*25c28e83SPiotr Jasiukajtis -9.006726958363837675e-17, 1.221846032972757623e+00,-1.061102121140269116e-16, 175*25c28e83SPiotr Jasiukajtis 1.225158793637145527e+00,-8.903533814269983429e-17, 1.228480536106870025e+00, 176*25c28e83SPiotr Jasiukajtis -1.898781631302529953e-17, 1.231811284734075862e+00, 7.389382471610050247e-17, 177*25c28e83SPiotr Jasiukajtis 1.235151063936933413e+00,-1.075524434430784138e-16, 1.238499898199816540e+00, 178*25c28e83SPiotr Jasiukajtis 2.767702055573967430e-17, 1.241857812073484002e+00, 4.658027591836936791e-17, 179*25c28e83SPiotr Jasiukajtis 1.245224830175257980e+00,-4.677240449846727500e-17, 1.248600977189204819e+00, 180*25c28e83SPiotr Jasiukajtis -8.261810999021963550e-17, 1.251986277866316222e+00, 4.834167152469897600e-17, 181*25c28e83SPiotr Jasiukajtis 1.255380757024691096e+00,-6.711389821296878419e-18, 1.258784439549716527e+00, 182*25c28e83SPiotr Jasiukajtis -8.421782587730599357e-17, 1.262197350394250739e+00,-3.084464887473846465e-17, 183*25c28e83SPiotr Jasiukajtis 1.265619514578806282e+00, 4.250577003450868637e-17, 1.269050957191733220e+00, 184*25c28e83SPiotr Jasiukajtis 2.667932131342186095e-18, 1.272491703389402762e+00,-1.057791626721242103e-17, 185*25c28e83SPiotr Jasiukajtis 1.275941778396392001e+00, 9.915430244214290330e-17, 1.279401207505669325e+00, 186*25c28e83SPiotr Jasiukajtis -9.759095008356062210e-17, 1.282870016078778264e+00, 1.713594918243560968e-17, 187*25c28e83SPiotr Jasiukajtis 1.286348229546025568e+00,-3.416955706936181976e-17, 1.289835873406665723e+00, 188*25c28e83SPiotr Jasiukajtis 8.949257530897591722e-17, 1.293332973229089466e+00,-2.974590443132751646e-17, 189*25c28e83SPiotr Jasiukajtis 1.296839554651009641e+00, 2.538250279488831496e-17, 1.300355643379650594e+00, 190*25c28e83SPiotr Jasiukajtis 5.678728102802217422e-17, 1.303881265191935812e+00, 8.647675598267871179e-17, 191*25c28e83SPiotr Jasiukajtis 1.307416445934677318e+00,-7.336645652878868892e-17, 1.310961211524764414e+00, 192*25c28e83SPiotr Jasiukajtis -7.181536135519453857e-17, 1.314515587949354636e+00, 2.267543315104585645e-17, 193*25c28e83SPiotr Jasiukajtis 1.318079601266064049e+00,-5.457955827149153502e-17, 1.321653277603157539e+00, 194*25c28e83SPiotr Jasiukajtis -2.480638245913021742e-17, 1.325236643159741323e+00,-2.858731210038861373e-17, 195*25c28e83SPiotr Jasiukajtis 1.328829724205954355e+00, 4.089086223910160052e-17, 1.332432547083161500e+00, 196*25c28e83SPiotr Jasiukajtis -5.101586630916743959e-17, 1.336045138204145832e+00,-5.891866356388801353e-17, 197*25c28e83SPiotr Jasiukajtis 1.339667524053302916e+00, 8.927282594831731984e-17, 1.343299731186835322e+00, 198*25c28e83SPiotr Jasiukajtis -5.802580890201437751e-17, 1.346941786232945804e+00, 3.224065101254679169e-17, 199*25c28e83SPiotr Jasiukajtis 1.350593715892034474e+00,-8.287110381462416533e-17, 1.354255546936892651e+00, 200*25c28e83SPiotr Jasiukajtis 7.700948379802989462e-17, 1.357927306212901142e+00,-9.529635744825188867e-17, 201*25c28e83SPiotr Jasiukajtis 1.361609020638224754e+00, 1.533787661270668046e-18, 1.365300717204011915e+00, 202*25c28e83SPiotr Jasiukajtis -1.000536312597476517e-16, 1.369002422974590516e+00, 9.593797919118848773e-17, 203*25c28e83SPiotr Jasiukajtis 1.372714165087668414e+00,-4.495960595234841262e-17, 1.376435970754530169e+00, 204*25c28e83SPiotr Jasiukajtis -6.898588935871801042e-17, 1.380167867260237990e+00, 1.051031457996998395e-16, 205*25c28e83SPiotr Jasiukajtis 1.383909881963832023e+00,-6.770511658794786287e-17, 1.387662042298529075e+00, 206*25c28e83SPiotr Jasiukajtis 8.422984274875415318e-17, 1.391424375771926236e+00,-4.906174865288989325e-17, 207*25c28e83SPiotr Jasiukajtis 1.395196909966200272e+00,-9.329336224225496552e-17, 1.398979672538311236e+00, 208*25c28e83SPiotr Jasiukajtis -9.614213209051323072e-17, 1.402772691220204759e+00,-5.295783249407989223e-17, 209*25c28e83SPiotr Jasiukajtis 1.406575993819015435e+00, 7.034914812136422188e-18, 1.410389608217270663e+00, 210*25c28e83SPiotr Jasiukajtis 4.166548728435062259e-17, 1.414213562373095145e+00,-9.667293313452913451e-17, 211*25c28e83SPiotr Jasiukajtis 1.418047884320415175e+00, 2.274438542185529452e-17, 1.421892602169165576e+00, 212*25c28e83SPiotr Jasiukajtis -1.607782891589024413e-17, 1.425747744105494208e+00, 9.880690758500607284e-17, 213*25c28e83SPiotr Jasiukajtis 1.429613338391970023e+00,-1.203164248905365518e-17, 1.433489413367788901e+00, 214*25c28e83SPiotr Jasiukajtis -5.802454243926826103e-17, 1.437375997448982368e+00,-4.204034016467556612e-17, 215*25c28e83SPiotr Jasiukajtis 1.441273119128625657e+00, 5.602503650878985675e-18, 1.445180806977046650e+00, 216*25c28e83SPiotr Jasiukajtis -3.023758134993987319e-17, 1.449099089642035043e+00,-6.259405000819309254e-17, 217*25c28e83SPiotr Jasiukajtis 1.453027995849052623e+00,-5.779948609396106102e-17, 1.456967554401443765e+00, 218*25c28e83SPiotr Jasiukajtis 5.648679453876998140e-17, 1.460917794180647045e+00,-5.600377186075215800e-17, 219*25c28e83SPiotr Jasiukajtis 1.464878744146405731e+00, 9.530767543587157319e-17, 1.468850433336981842e+00, 220*25c28e83SPiotr Jasiukajtis 8.465882756533627608e-17, 1.472832890869367528e+00, 6.691774081940589372e-17, 221*25c28e83SPiotr Jasiukajtis 1.476826145939499346e+00,-3.483994556892795796e-17, 1.480830227822471867e+00, 222*25c28e83SPiotr Jasiukajtis -9.686952102630618578e-17, 1.484845165872752393e+00, 1.078008676440748076e-16, 223*25c28e83SPiotr Jasiukajtis 1.488870989524397004e+00, 6.155367157742871330e-17, 1.492907728291264835e+00, 224*25c28e83SPiotr Jasiukajtis 1.419292015428403577e-17, 1.496955411767235455e+00,-2.861663253899158211e-17, 225*25c28e83SPiotr Jasiukajtis 1.501014069626425584e+00,-6.413767275790235039e-17, 1.505083731623406473e+00, 226*25c28e83SPiotr Jasiukajtis 7.074710613582846364e-17, 1.509164427593422841e+00,-1.016455327754295039e-16, 227*25c28e83SPiotr Jasiukajtis 1.513256187452609813e+00, 8.884497851338712091e-17, 1.517359041198214742e+00, 228*25c28e83SPiotr Jasiukajtis -4.308699472043340801e-17, 1.521473018908814590e+00,-5.996387675945683420e-18, 229*25c28e83SPiotr Jasiukajtis 1.525598150744538417e+00,-1.102494171234256094e-16, 1.529734466947286986e+00, 230*25c28e83SPiotr Jasiukajtis 3.785792115157219653e-17, 1.533881997840955913e+00, 8.875226844438446141e-17, 231*25c28e83SPiotr Jasiukajtis 1.538040773831656827e+00, 1.017467235116135806e-16, 1.542210825407940744e+00, 232*25c28e83SPiotr Jasiukajtis 7.949834809697620856e-17, 1.546392183141021448e+00, 1.068396000565721980e-16, 233*25c28e83SPiotr Jasiukajtis 1.550584877684999974e+00,-1.460070659068938518e-17, 1.554788939777088652e+00, 234*25c28e83SPiotr Jasiukajtis -8.003161350116035641e-17, 1.559004400237836929e+00, 3.781207053357527502e-17, 235*25c28e83SPiotr Jasiukajtis 1.563231289971357629e+00, 7.484777645590734389e-17, 1.567469639965552997e+00, 236*25c28e83SPiotr Jasiukajtis -1.035206176884972199e-16, 1.571719481292341403e+00,-3.342984004687200069e-17, 237*25c28e83SPiotr Jasiukajtis 1.575980845107886497e+00,-1.013691647127830398e-17, 1.580253762652824578e+00, 238*25c28e83SPiotr Jasiukajtis -5.163402929554468062e-17, 1.584538265252493749e+00,-1.933771703458570293e-17, 239*25c28e83SPiotr Jasiukajtis 1.588834384317163950e+00,-5.994950118824479401e-18, 1.593142151342266999e+00, 240*25c28e83SPiotr Jasiukajtis -1.009440654231196372e-16, 1.597461597908627073e+00, 2.486839279622099613e-17, 241*25c28e83SPiotr Jasiukajtis 1.601792755682693414e+00,-6.054917453527784343e-17, 1.606135656416771029e+00, 242*25c28e83SPiotr Jasiukajtis -1.035454528805999526e-16, 1.610490331949254283e+00, 2.470719256979788785e-17, 243*25c28e83SPiotr Jasiukajtis 1.614856814204860713e+00,-7.316663399125123263e-17, 1.619235135194863728e+00, 244*25c28e83SPiotr Jasiukajtis 2.094133415422909241e-17, 1.623625327017328868e+00,-3.584512851414474710e-17, 245*25c28e83SPiotr Jasiukajtis 1.628027421857347834e+00,-6.712955084707084086e-17, 1.632441451987274972e+00, 246*25c28e83SPiotr Jasiukajtis 9.852819230429992964e-17, 1.636867449766964411e+00, 7.698325071319875575e-17, 247*25c28e83SPiotr Jasiukajtis 1.641305447644006321e+00,-9.247568737640705508e-17, 1.645755478153964946e+00, 248*25c28e83SPiotr Jasiukajtis -1.012567991367477260e-16, 1.650217573920617742e+00, 9.133279588729904190e-18, 249*25c28e83SPiotr Jasiukajtis 1.654691767656194301e+00, 9.643294303196028661e-17, 1.659178092161616158e+00, 250*25c28e83SPiotr Jasiukajtis -7.275545550823050654e-17, 1.663676580326736376e+00, 5.890992696713099670e-17, 251*25c28e83SPiotr Jasiukajtis 1.668187265130582464e+00, 4.269178019570615091e-17, 1.672710179641596628e+00, 252*25c28e83SPiotr Jasiukajtis -5.476715964599563076e-17, 1.677245357017878469e+00, 8.303949509950732785e-17, 253*25c28e83SPiotr Jasiukajtis 1.681792830507429004e+00, 8.199010020581496520e-17, 1.686352633448393368e+00, 254*25c28e83SPiotr Jasiukajtis -7.181463278358010675e-17, 1.690924799269305279e+00,-9.669671474394880166e-17, 255*25c28e83SPiotr Jasiukajtis 1.695509361489332623e+00, 7.238416872845166641e-17, 1.700106353718523478e+00, 256*25c28e83SPiotr Jasiukajtis -8.023719370397700246e-18, 1.704715809658051251e+00,-2.728883284797281563e-17, 257*25c28e83SPiotr Jasiukajtis 1.709337763100462926e+00,-9.868779456632931076e-17, 1.713972247929925974e+00, 258*25c28e83SPiotr Jasiukajtis 6.473975107753367064e-17, 1.718619298122477934e+00,-1.851380418263110988e-17, 259*25c28e83SPiotr Jasiukajtis 1.723278947746273992e+00,-9.522123800393799963e-17, 1.727951230961837670e+00, 260*25c28e83SPiotr Jasiukajtis -1.075098186120464245e-16, 1.732636182022311067e+00,-1.698051074315415494e-18, 261*25c28e83SPiotr Jasiukajtis 1.737333835273706217e+00, 3.164389299292956947e-17, 1.742044225155156445e+00, 262*25c28e83SPiotr Jasiukajtis -1.525959118950788792e-18, 1.746767386199169048e+00,-1.075229048350751450e-16, 263*25c28e83SPiotr Jasiukajtis 1.751503353031878207e+00,-5.124450420596724659e-17, 1.756252160373299454e+00, 264*25c28e83SPiotr Jasiukajtis 2.960140695448873307e-17, 1.761013843037583904e+00,-7.943253125039227711e-17, 265*25c28e83SPiotr Jasiukajtis 1.765788435933272726e+00, 9.461315018083267867e-17, 1.770575974063554714e+00, 266*25c28e83SPiotr Jasiukajtis 5.961794510040555848e-17, 1.775376492526521188e+00, 6.429731796556572034e-17, 267*25c28e83SPiotr Jasiukajtis 1.780190026515424462e+00,-5.284627289091617365e-17, 1.785016611318934965e+00, 268*25c28e83SPiotr Jasiukajtis 1.533040012103131382e-17, 1.789856282321401038e+00,-4.154354660683350387e-17, 269*25c28e83SPiotr Jasiukajtis 1.794709075003107168e+00, 1.822745842791208677e-17, 1.799575024940535117e+00, 270*25c28e83SPiotr Jasiukajtis -2.526889233358897644e-17, 1.804454167806623932e+00,-5.177222408793317883e-17, 271*25c28e83SPiotr Jasiukajtis 1.809346539371031959e+00,-9.032641402450029682e-17, 1.814252175500398856e+00, 272*25c28e83SPiotr Jasiukajtis -9.969531538920348820e-17, 1.819171112158608494e+00, 7.402676901145838890e-17, 273*25c28e83SPiotr Jasiukajtis 1.824103385407053413e+00,-1.015962786227708306e-16, 1.829049031404897274e+00, 274*25c28e83SPiotr Jasiukajtis 6.889192908835695637e-17, 1.834008086409342431e+00, 3.283107224245627204e-17, 275*25c28e83SPiotr Jasiukajtis 1.838980586775893711e+00, 6.918969740272511942e-18, 1.843966568958625984e+00, 276*25c28e83SPiotr Jasiukajtis -5.939742026949964550e-17, 1.848966069510450838e+00, 9.027580446261089288e-17, 277*25c28e83SPiotr Jasiukajtis 1.853979125083385471e+00, 9.761887490727593538e-17, 1.859005772428820480e+00, 278*25c28e83SPiotr Jasiukajtis -9.528705461989940687e-17, 1.864046048397788979e+00, 6.540912680620571711e-17, 279*25c28e83SPiotr Jasiukajtis 1.869099989941238604e+00,-9.938505214255067083e-17, 1.874167634110299963e+00, 280*25c28e83SPiotr Jasiukajtis -6.122763413004142562e-17, 1.879249018056560194e+00,-1.622631555783584478e-17, 281*25c28e83SPiotr Jasiukajtis 1.884344179032334532e+00,-8.226593125533710906e-17, 1.889453154390939194e+00, 282*25c28e83SPiotr Jasiukajtis -9.005168285059126718e-17, 1.894575981586965607e+00, 3.403403535216529671e-17, 283*25c28e83SPiotr Jasiukajtis 1.899712698176555303e+00,-3.859739769378514323e-17, 1.904863341817674138e+00, 284*25c28e83SPiotr Jasiukajtis 6.533857514718278629e-17, 1.910027950270389852e+00,-5.909688006744060237e-17, 285*25c28e83SPiotr Jasiukajtis 1.915206561397147400e+00,-1.061994605619596264e-16, 1.920399213163047403e+00, 286*25c28e83SPiotr Jasiukajtis 7.116681540630314186e-17, 1.925605943636125028e+00,-9.914963769693740927e-17, 287*25c28e83SPiotr Jasiukajtis 1.930826790987627106e+00, 6.167149706169109553e-17, 1.936061793492294347e+00, 288*25c28e83SPiotr Jasiukajtis 1.033238596067632574e-16, 1.941310989528640452e+00,-6.638029891621487990e-17, 289*25c28e83SPiotr Jasiukajtis 1.946574417579233218e+00, 6.811022349533877184e-17, 1.951852116230978318e+00, 290*25c28e83SPiotr Jasiukajtis -2.199016969979351086e-17, 1.957144124175400179e+00, 8.960767791036667768e-17, 291*25c28e83SPiotr Jasiukajtis 1.962450480208927317e+00, 1.097684400091354695e-16, 1.967771223233175881e+00, 292*25c28e83SPiotr Jasiukajtis -1.031492801153113151e-16, 1.973106392255234320e+00,-7.451617863956037486e-18, 293*25c28e83SPiotr Jasiukajtis 1.978456026387950928e+00, 4.038875310927816657e-17, 1.983820164850219392e+00, 294*25c28e83SPiotr Jasiukajtis -2.203454412391062657e-17, 1.989198846967266343e+00, 8.205132638369199416e-18, 295*25c28e83SPiotr Jasiukajtis 1.994592112170940235e+00, 1.790971035200264509e-17 296*25c28e83SPiotr Jasiukajtis }; 297*25c28e83SPiotr Jasiukajtis 298*25c28e83SPiotr Jasiukajtis static const double __TBL_log2[] = { 299*25c28e83SPiotr Jasiukajtis /* __TBL_log2[2*i] = high order rounded 32 bits log2(1+i/256)*256, i = [0, 255] */ 300*25c28e83SPiotr Jasiukajtis /* __TBL_log2[2*i+1] = low order least bits log2(1+i/256)*256, i = [0, 255] */ 301*25c28e83SPiotr Jasiukajtis 0.000000000000000000e+00, 0.000000000000000000e+00, 1.439884185791015625e+00, 302*25c28e83SPiotr Jasiukajtis 4.078417797464839152e-07, 2.874177932739257812e+00,-5.443862030060025621e-07, 303*25c28e83SPiotr Jasiukajtis 4.302921295166015625e+00, 3.525917800357419922e-07, 5.726161956787109375e+00, 304*25c28e83SPiotr Jasiukajtis -1.821502755258614180e-06, 7.143936157226562500e+00,-1.035336134691423741e-06, 305*25c28e83SPiotr Jasiukajtis 8.556289672851562500e+00,-1.279264291071495652e-06, 9.963264465332031250e+00, 306*25c28e83SPiotr Jasiukajtis -3.206502629414843101e-06, 1.136489105224609375e+01, 3.503517986289194222e-06, 307*25c28e83SPiotr Jasiukajtis 1.276123046875000000e+01,-1.809406249049319022e-06, 1.415230560302734375e+01, 308*25c28e83SPiotr Jasiukajtis -2.114722805833714926e-06, 1.553816223144531250e+01,-3.719431504776986979e-06, 309*25c28e83SPiotr Jasiukajtis 1.691883850097656250e+01,-5.743786819670105240e-06, 1.829435729980468750e+01, 310*25c28e83SPiotr Jasiukajtis 7.514691093524705578e-06, 1.966479492187500000e+01,-2.076862291588726520e-06, 311*25c28e83SPiotr Jasiukajtis 2.103015136718750000e+01, 3.219403619538604258e-06, 2.239048767089843750e+01, 312*25c28e83SPiotr Jasiukajtis -3.108115489869591032e-07, 2.374583435058593750e+01,-6.275103710481114264e-06, 313*25c28e83SPiotr Jasiukajtis 2.509620666503906250e+01, 6.572855776743687178e-06, 2.644168090820312500e+01, 314*25c28e83SPiotr Jasiukajtis -1.954725505303359537e-06, 2.778225708007812500e+01, 3.855133152759458770e-06, 315*25c28e83SPiotr Jasiukajtis 2.911799621582031250e+01,-1.707228100041815487e-06, 3.044891357421875000e+01, 316*25c28e83SPiotr Jasiukajtis 1.042999152333371737e-06, 3.177505493164062500e+01, 8.966313933586820042e-07, 317*25c28e83SPiotr Jasiukajtis 3.309646606445312500e+01,-1.372654171244005427e-05, 3.441314697265625000e+01, 318*25c28e83SPiotr Jasiukajtis -8.996099168734074844e-06, 3.572515869140625000e+01,-1.247731510027211536e-05, 319*25c28e83SPiotr Jasiukajtis 3.703250122070312500e+01, 8.944258749129049106e-06, 3.833526611328125000e+01, 320*25c28e83SPiotr Jasiukajtis -3.520082642279872716e-06, 3.963342285156250000e+01, 1.306577612991810031e-05, 321*25c28e83SPiotr Jasiukajtis 4.092706298828125000e+01,-7.730135593513790229e-07, 4.221618652343750000e+01, 322*25c28e83SPiotr Jasiukajtis -1.329446142304436745e-05, 4.350079345703125000e+01, 6.912200714904314733e-06, 323*25c28e83SPiotr Jasiukajtis 4.478097534179687500e+01,-6.216230979739182064e-07, 4.605673217773437500e+01, 324*25c28e83SPiotr Jasiukajtis -5.133911151040936670e-06, 4.732809448242187500e+01,-6.697901206512330627e-06, 325*25c28e83SPiotr Jasiukajtis 4.859509277343750000e+01,-5.700153089154811841e-06, 4.985775756835937500e+01, 326*25c28e83SPiotr Jasiukajtis -2.836263919120346801e-06, 5.111611938476562500e+01, 8.933436604624454391e-07, 327*25c28e83SPiotr Jasiukajtis 5.237020874023437500e+01, 4.187561748309498307e-06, 5.362005615234375000e+01, 328*25c28e83SPiotr Jasiukajtis 5.448667394155597532e-06, 5.486569213867187500e+01, 2.786324169943508531e-06, 329*25c28e83SPiotr Jasiukajtis 5.610714721679687500e+01,-5.978483512667373796e-06, 5.734442138671875000e+01, 330*25c28e83SPiotr Jasiukajtis 7.207996138368885843e-06, 5.857757568359375000e+01, 9.083351754561760127e-06, 331*25c28e83SPiotr Jasiukajtis 5.980664062500000000e+01,-3.374516276140515786e-06, 6.103161621093750000e+01, 332*25c28e83SPiotr Jasiukajtis -2.943717299925017200e-06, 6.225253295898437500e+01, 6.810091060168101732e-06, 333*25c28e83SPiotr Jasiukajtis 6.346945190429687500e+01,-8.462738988588859704e-06, 6.468237304687500000e+01, 334*25c28e83SPiotr Jasiukajtis -2.233961135216831566e-05, 6.589129638671875000e+01,-8.657399896582645111e-06, 335*25c28e83SPiotr Jasiukajtis 6.709625244140625000e+01, 2.797335967336006296e-05, 6.829736328125000000e+01, 336*25c28e83SPiotr Jasiukajtis -8.863355250907819214e-06, 6.949450683593750000e+01, 2.830758238800374038e-05, 337*25c28e83SPiotr Jasiukajtis 7.068786621093750000e+01,-1.846073268549083018e-05, 7.187731933593750000e+01, 338*25c28e83SPiotr Jasiukajtis -2.182503249464459606e-06, 7.306298828125000000e+01,-2.025251442448625989e-05, 339*25c28e83SPiotr Jasiukajtis 7.424481201171875000e+01, 1.280303154355201204e-05, 7.542291259765625000e+01, 340*25c28e83SPiotr Jasiukajtis -8.813997363590295654e-07, 7.659722900390625000e+01, 2.370323712746426047e-05, 341*25c28e83SPiotr Jasiukajtis 7.776788330078125000e+01,-1.176744290134661421e-05, 7.893481445312500000e+01, 342*25c28e83SPiotr Jasiukajtis -2.273743674288609119e-05, 8.009802246093750000e+01, 1.409185747234803696e-05, 343*25c28e83SPiotr Jasiukajtis 8.125762939453125000e+01,-2.707246895087010889e-07, 8.241357421875000000e+01, 344*25c28e83SPiotr Jasiukajtis 1.807241476105480180e-05, 8.356597900390625000e+01,-3.030059664889450720e-05, 345*25c28e83SPiotr Jasiukajtis 8.471472167968750000e+01,-8.823455531875539245e-07, 8.585992431640625000e+01, 346*25c28e83SPiotr Jasiukajtis 6.485238524924182146e-06, 8.700158691406250000e+01, 1.382440142980862947e-05, 347*25c28e83SPiotr Jasiukajtis 8.813977050781250000e+01,-1.808136338482881111e-05, 8.927441406250000000e+01, 348*25c28e83SPiotr Jasiukajtis -6.579344146543672011e-06, 9.040557861328125000e+01, 8.714227880222726313e-06, 349*25c28e83SPiotr Jasiukajtis 9.153332519531250000e+01,-1.201308307454951138e-05, 9.265759277343750000e+01, 350*25c28e83SPiotr Jasiukajtis 1.330278431878087205e-05, 9.377850341796875000e+01,-1.657103990890600482e-05, 351*25c28e83SPiotr Jasiukajtis 9.489599609375000000e+01,-1.995110226941163424e-05, 9.601007080078125000e+01, 352*25c28e83SPiotr Jasiukajtis 2.362403148762806632e-05, 9.712084960937500000e+01, 1.236086810905991142e-05, 353*25c28e83SPiotr Jasiukajtis 9.822827148437500000e+01, 2.738898236946465744e-05, 9.933239746093750000e+01, 354*25c28e83SPiotr Jasiukajtis 2.758741700388469572e-05, 1.004332885742187500e+02,-2.834285611604269955e-05, 355*25c28e83SPiotr Jasiukajtis 1.015308227539062500e+02, 1.228649517068771375e-06, 1.026251220703125000e+02, 356*25c28e83SPiotr Jasiukajtis 1.361792668612316888e-05, 1.037161865234375000e+02, 2.803946653578170389e-05, 357*25c28e83SPiotr Jasiukajtis 1.048040771484375000e+02, 2.502814149567842806e-06, 1.058887329101562500e+02, 358*25c28e83SPiotr Jasiukajtis 1.692003190104140317e-05, 1.069702148437500000e+02, 2.896703985131545672e-05, 359*25c28e83SPiotr Jasiukajtis 1.080485839843750000e+02,-3.844135045484567362e-06, 1.091237792968750000e+02, 360*25c28e83SPiotr Jasiukajtis -2.093137927645659717e-06, 1.101958618164062500e+02,-8.590030211185738579e-06, 361*25c28e83SPiotr Jasiukajtis 1.112648315429687500e+02,-5.267967244023324300e-06, 1.123306884765625000e+02, 362*25c28e83SPiotr Jasiukajtis 2.578347229232600646e-05, 1.133935546875000000e+02,-1.975022555464358195e-05, 363*25c28e83SPiotr Jasiukajtis 1.144533081054687500e+02,-2.195797778964440179e-06, 1.155100708007812500e+02, 364*25c28e83SPiotr Jasiukajtis -2.617170507638525077e-05, 1.165637817382812500e+02,-1.334031370958194516e-05, 365*25c28e83SPiotr Jasiukajtis 1.176145019531250000e+02,-7.581976902412963145e-06, 1.186622314453125000e+02, 366*25c28e83SPiotr Jasiukajtis 8.112109654298731037e-06, 1.197070312500000000e+02,-1.042875265529314613e-05, 367*25c28e83SPiotr Jasiukajtis 1.207488403320312500e+02, 1.455233211877492951e-05, 1.217877807617187500e+02, 368*25c28e83SPiotr Jasiukajtis -2.243432092472914265e-05, 1.228237304687500000e+02, 1.712269952247034061e-05, 369*25c28e83SPiotr Jasiukajtis 1.238568115234375000e+02, 2.745621214456745937e-05, 1.248870239257812500e+02, 370*25c28e83SPiotr Jasiukajtis 2.473291989440979066e-05, 1.259143676757812500e+02, 2.498461547595911484e-05, 371*25c28e83SPiotr Jasiukajtis 1.269389038085937500e+02,-1.692547797717771941e-05, 1.279605712890625000e+02, 372*25c28e83SPiotr Jasiukajtis -2.419576192770340594e-05, 1.289793701171875000e+02, 1.880972467762623192e-05, 373*25c28e83SPiotr Jasiukajtis 1.299954833984375000e+02,-5.550757125543327248e-05, 1.310086669921875000e+02, 374*25c28e83SPiotr Jasiukajtis 1.237226167189998996e-05, 1.320191650390625000e+02,-6.438347630770959254e-06, 375*25c28e83SPiotr Jasiukajtis 1.330268554687500000e+02, 2.525911246920619613e-05, 1.340318603515625000e+02, 376*25c28e83SPiotr Jasiukajtis 3.990327953073019333e-07, 1.350340576171875000e+02, 5.593427389035480335e-05, 377*25c28e83SPiotr Jasiukajtis 1.360336914062500000e+02,-3.751407409478960320e-05, 1.370305175781250000e+02, 378*25c28e83SPiotr Jasiukajtis -2.116319935859897563e-05, 1.380246582031250000e+02,-2.559468964093475045e-06, 379*25c28e83SPiotr Jasiukajtis 1.390161132812500000e+02, 3.270409087092109593e-05, 1.400050048828125000e+02, 380*25c28e83SPiotr Jasiukajtis -2.315157751389992129e-05, 1.409912109375000000e+02,-3.387938973438343638e-05, 381*25c28e83SPiotr Jasiukajtis 1.419747314453125000e+02, 1.458416266727572812e-05, 1.429556884765625000e+02, 382*25c28e83SPiotr Jasiukajtis 1.412021555596584681e-05, 1.439340820312500000e+02,-2.143065540113838312e-05, 383*25c28e83SPiotr Jasiukajtis 1.449097900390625000e+02, 4.373273697503468317e-05, 1.458830566406250000e+02, 384*25c28e83SPiotr Jasiukajtis -2.090790235253405790e-05, 1.468536376953125000e+02, 4.230297794089183646e-05, 385*25c28e83SPiotr Jasiukajtis 1.478217773437500000e+02, 2.633401664450247309e-06, 1.487873535156250000e+02, 386*25c28e83SPiotr Jasiukajtis -4.542835986281740771e-06, 1.497503662109375000e+02, 3.397367848245215483e-05, 387*25c28e83SPiotr Jasiukajtis 1.507109375000000000e+02, 9.209059510146982590e-06, 1.516689453125000000e+02, 388*25c28e83SPiotr Jasiukajtis 5.622812858742714859e-05, 1.526246337890625000e+02,-5.621609346274134244e-05, 389*25c28e83SPiotr Jasiukajtis 1.535776367187500000e+02, 5.088115468603551539e-05, 1.545283203125000000e+02, 390*25c28e83SPiotr Jasiukajtis 2.400396513473623342e-05, 1.554765625000000000e+02,-2.180099663431456814e-06, 391*25c28e83SPiotr Jasiukajtis 1.564223632812500000e+02,-1.517056781617965675e-05, 1.573657226562500000e+02, 392*25c28e83SPiotr Jasiukajtis -2.562756696989711716e-06, 1.583066406250000000e+02, 4.795320325388065854e-05, 393*25c28e83SPiotr Jasiukajtis 1.592452392578125000e+02, 2.652301982429665372e-05, 1.601815185546875000e+02, 394*25c28e83SPiotr Jasiukajtis -5.473018439029181240e-05, 1.611152343750000000e+02, 6.036538006249134820e-05, 395*25c28e83SPiotr Jasiukajtis 1.620467529296875000e+02, 1.753890969321481711e-05, 1.629759521484375000e+02, 396*25c28e83SPiotr Jasiukajtis -4.928926339732922490e-05, 1.639027099609375000e+02,-6.288016979631557560e-06, 397*25c28e83SPiotr Jasiukajtis 1.648271484375000000e+02, 3.614482952210960361e-05, 1.657493896484375000e+02, 398*25c28e83SPiotr Jasiukajtis -3.247597790375142114e-05, 1.666691894531250000e+02, 4.348868072528205213e-05, 399*25c28e83SPiotr Jasiukajtis 1.675867919921875000e+02, 3.131097214651595330e-05, 1.685021972656250000e+02, 400*25c28e83SPiotr Jasiukajtis -5.768116554728405733e-05, 1.694151611328125000e+02, 3.189681619086343127e-05, 401*25c28e83SPiotr Jasiukajtis 1.703260498046875000e+02,-5.500528238559059116e-05, 1.712344970703125000e+02, 402*25c28e83SPiotr Jasiukajtis 5.890184674174263693e-05, 1.721408691406250000e+02, 1.840407787096519837e-05, 403*25c28e83SPiotr Jasiukajtis 1.730450439453125000e+02,-4.351222480150346831e-05, 1.739468994140625000e+02, 404*25c28e83SPiotr Jasiukajtis 6.059331686505290421e-06, 1.748465576171875000e+02, 5.580532332169584454e-05, 405*25c28e83SPiotr Jasiukajtis 1.757441406250000000e+02,-5.666096094448416139e-06, 1.766395263671875000e+02, 406*25c28e83SPiotr Jasiukajtis -4.568380948624016041e-05, 1.775327148437500000e+02,-5.372392273978838048e-05, 407*25c28e83SPiotr Jasiukajtis 1.784237060546875000e+02,-1.933871000131713187e-05, 1.793126220703125000e+02, 408*25c28e83SPiotr Jasiukajtis -5.422619290693841471e-05, 1.801993408203125000e+02,-2.601847861521447132e-05, 409*25c28e83SPiotr Jasiukajtis 1.810839843750000000e+02,-4.656229401600182454e-05, 1.819664306640625000e+02, 410*25c28e83SPiotr Jasiukajtis 1.636297150881445295e-05, 1.828468017578125000e+02, 5.076471489501210225e-05, 411*25c28e83SPiotr Jasiukajtis 1.837252197265625000e+02,-5.542156510357154555e-05, 1.846014404296875000e+02, 412*25c28e83SPiotr Jasiukajtis -4.812064810565531807e-05, 1.854755859375000000e+02,-3.953879286781995545e-05, 413*25c28e83SPiotr Jasiukajtis 1.863476562500000000e+02,-1.988182101010412125e-05, 1.872176513671875000e+02, 414*25c28e83SPiotr Jasiukajtis 2.057522891062264376e-05, 1.880856933593750000e+02,-3.058156040982771239e-05, 415*25c28e83SPiotr Jasiukajtis 1.889516601562500000e+02,-4.169340446171797184e-05, 1.898155517578125000e+02, 416*25c28e83SPiotr Jasiukajtis -3.239118881346662872e-06, 1.906774902343750000e+02,-2.783449132689922134e-05, 417*25c28e83SPiotr Jasiukajtis 1.915373535156250000e+02, 1.597927683340914293e-05, 1.923952636718750000e+02, 418*25c28e83SPiotr Jasiukajtis 1.545493412281261116e-05, 1.932512207031250000e+02,-2.014927705264352875e-05, 419*25c28e83SPiotr Jasiukajtis 1.941051025390625000e+02, 4.043097907577914080e-05, 1.949571533203125000e+02, 420*25c28e83SPiotr Jasiukajtis -3.781452579504048975e-05, 1.958071289062500000e+02,-1.677810793588779092e-06, 421*25c28e83SPiotr Jasiukajtis 1.966551513671875000e+02, 3.577570564777057149e-05, 1.975013427734375000e+02, 422*25c28e83SPiotr Jasiukajtis -3.858128431828155999e-05, 1.983454589843750000e+02, 2.827352539329734468e-05, 423*25c28e83SPiotr Jasiukajtis 1.991877441406250000e+02, 1.020426695132691908e-06, 2.000280761718750000e+02, 424*25c28e83SPiotr Jasiukajtis 1.049043785864183866e-05, 2.008665771484375000e+02,-5.668571223208539910e-05, 425*25c28e83SPiotr Jasiukajtis 2.017030029296875000e+02, 5.227451898157462205e-05, 2.025377197265625000e+02, 426*25c28e83SPiotr Jasiukajtis -2.025647781341857894e-05, 2.033704833984375000e+02,-2.161281037339224341e-05, 427*25c28e83SPiotr Jasiukajtis 2.042012939453125000e+02, 5.667325008632565576e-05, 2.050303955078125000e+02, 428*25c28e83SPiotr Jasiukajtis -2.112821448834358837e-05, 2.058575439453125000e+02,-2.522383155215216853e-06, 429*25c28e83SPiotr Jasiukajtis 2.066828613281250000e+02,-1.281378348494855858e-06, 2.075063476562500000e+02, 430*25c28e83SPiotr Jasiukajtis -9.162516382743561384e-06, 2.083280029296875000e+02,-1.797812601298608335e-05, 431*25c28e83SPiotr Jasiukajtis 2.091478271484375000e+02,-1.959505997696247453e-05, 2.099658203125000000e+02, 432*25c28e83SPiotr Jasiukajtis -5.934211946670452627e-06, 2.107819824218750000e+02, 3.102996118252714271e-05, 433*25c28e83SPiotr Jasiukajtis 2.115964355468750000e+02,-2.280040076415178584e-05, 2.124090576171875000e+02, 434*25c28e83SPiotr Jasiukajtis -3.743515649437846729e-05, 2.132198486328125000e+02,-5.006638631136701490e-06, 435*25c28e83SPiotr Jasiukajtis 2.140289306640625000e+02,-3.976919665668718942e-05, 2.148361816406250000e+02, 436*25c28e83SPiotr Jasiukajtis -1.188780735169185652e-05, 2.156417236328125000e+02,-3.571887766413048520e-05, 437*25c28e83SPiotr Jasiukajtis 2.164454345703125000e+02, 1.847144755636210490e-05, 2.172474365234375000e+02, 438*25c28e83SPiotr Jasiukajtis 3.622647302213163157e-05, 2.180477294921875000e+02, 2.511032323154433900e-05, 439*25c28e83SPiotr Jasiukajtis 2.188463134765625000e+02,-7.361941985081681848e-06, 2.196431884765625000e+02, 440*25c28e83SPiotr Jasiukajtis -5.372390403709574017e-05, 2.204382324218750000e+02, 1.551294579696132803e-05, 441*25c28e83SPiotr Jasiukajtis 2.212316894531250000e+02,-3.642162925932327343e-05, 2.220233154296875000e+02, 442*25c28e83SPiotr Jasiukajtis 4.193598594979618241e-05, 2.228133544921875000e+02, 1.372116405796589833e-05, 443*25c28e83SPiotr Jasiukajtis 2.236016845703125000e+02, 8.233623894335039537e-06, 2.243883056640625000e+02, 444*25c28e83SPiotr Jasiukajtis 3.265657742833052654e-05, 2.251733398437500000e+02,-2.794287750390687326e-05, 445*25c28e83SPiotr Jasiukajtis 2.259566650390625000e+02,-4.440243113774530265e-05, 2.267382812500000000e+02, 446*25c28e83SPiotr Jasiukajtis -9.675114830058622014e-06, 2.275183105468750000e+02,-3.882892066889445600e-05, 447*25c28e83SPiotr Jasiukajtis 2.282966308593750000e+02,-2.835487591479255673e-06, 2.290733642578125000e+02, 448*25c28e83SPiotr Jasiukajtis -1.685097895998181422e-05, 2.298483886718750000e+02, 4.806553595480019518e-05, 449*25c28e83SPiotr Jasiukajtis 2.306219482421875000e+02,-4.539911586906436716e-05, 2.313937988281250000e+02, 450*25c28e83SPiotr Jasiukajtis -4.631966285757620260e-05, 2.321639404296875000e+02, 5.204609324350696002e-05, 451*25c28e83SPiotr Jasiukajtis 2.329326171875000000e+02, 1.225763073721718197e-05, 2.336997070312500000e+02, 452*25c28e83SPiotr Jasiukajtis -3.695637982554016382e-05, 2.344650878906250000e+02, 3.309133292926460016e-05, 453*25c28e83SPiotr Jasiukajtis 2.352290039062500000e+02,-1.516395380482592629e-05, 2.359913330078125000e+02, 454*25c28e83SPiotr Jasiukajtis -5.311674305290968619e-05, 2.367519531250000000e+02, 4.779807991226078768e-05, 455*25c28e83SPiotr Jasiukajtis 2.375111083984375000e+02, 4.989464209345647548e-05, 2.382687988281250000e+02, 456*25c28e83SPiotr Jasiukajtis -4.041202611322311408e-05, 2.390247802734375000e+02, 2.739433433590848536e-05, 457*25c28e83SPiotr Jasiukajtis 2.397792968750000000e+02, 1.550965806406508966e-05, 2.405322265625000000e+02, 458*25c28e83SPiotr Jasiukajtis 5.230206142425020257e-05, 2.412836914062500000e+02, 2.196059540790264514e-05, 459*25c28e83SPiotr Jasiukajtis 2.420335693359375000e+02, 5.277680785141730338e-05, 2.427819824218750000e+02, 460*25c28e83SPiotr Jasiukajtis 2.886380247947272558e-05, 2.435289306640625000e+02,-4.363251767645384661e-05, 461*25c28e83SPiotr Jasiukajtis 2.442742919921875000e+02,-3.653314744654563199e-05, 2.450180664062500000e+02, 462*25c28e83SPiotr Jasiukajtis 5.623369525922526825e-05, 2.457604980468750000e+02,-3.437446279919778004e-06, 463*25c28e83SPiotr Jasiukajtis 2.465013427734375000e+02, 3.459290119679066472e-05, 2.472407226562500000e+02, 464*25c28e83SPiotr Jasiukajtis 5.421724428316440202e-05, 2.479787597656250000e+02,-6.070765164808318435e-05, 465*25c28e83SPiotr Jasiukajtis 2.487152099609375000e+02,-6.014953987030989107e-05, 2.494501953125000000e+02, 466*25c28e83SPiotr Jasiukajtis -6.032228506450037554e-05, 2.501837158203125000e+02,-5.540433388359054134e-05, 467*25c28e83SPiotr Jasiukajtis 2.509157714843750000e+02,-3.960875078622925214e-05, 2.516463623046875000e+02, 468*25c28e83SPiotr Jasiukajtis -7.182944107105660894e-06, 2.523754882812500000e+02, 4.759160516857532540e-05, 469*25c28e83SPiotr Jasiukajtis 2.531032714843750000e+02, 8.329299458439681639e-06, 2.538295898437500000e+02, 470*25c28e83SPiotr Jasiukajtis 2.751627995643241118e-06, 2.545544433593750000e+02, 3.647649263201999678e-05, 471*25c28e83SPiotr Jasiukajtis 2.552779541015625000e+02,-6.981531437649667064e-06 472*25c28e83SPiotr Jasiukajtis }; 473*25c28e83SPiotr Jasiukajtis 474*25c28e83SPiotr Jasiukajtis static const unsigned long long LCONST[] = { 475*25c28e83SPiotr Jasiukajtis 0x3c90000000000000ULL, /* 2**(-54) = 5.551115123125782702e-17 */ 476*25c28e83SPiotr Jasiukajtis 0x3ff0000000000000ULL, /* DONE = 1.0 */ 477*25c28e83SPiotr Jasiukajtis 0x4330000000000000ULL, /* DVAIN52 = 2**52 = 4.503599627370496e15 */ 478*25c28e83SPiotr Jasiukajtis 0xffffffff00000000ULL, /* 0xffffffff00000000 */ 479*25c28e83SPiotr Jasiukajtis 0x000fffffffffffffULL, /* 0x000fffffffffffff */ 480*25c28e83SPiotr Jasiukajtis 0x0000080000000000ULL, /* 0x0000080000000000 */ 481*25c28e83SPiotr Jasiukajtis 0xfffff00000000000ULL, /* 0xfffff00000000000 */ 482*25c28e83SPiotr Jasiukajtis 0x0000000000000000ULL, /* DZERO = 0.0 */ 483*25c28e83SPiotr Jasiukajtis 0x4062776d8ce329bdULL, /* KA5 = 5.77078604860893737986e-01*256 */ 484*25c28e83SPiotr Jasiukajtis 0x406ec709dc39fc99ULL, /* KA3 = 9.61796693925765549423e-01*256 */ 485*25c28e83SPiotr Jasiukajtis 0x3f6d94ae0bf85de6ULL, /* KA1_LO = 1.41052154268147309568e-05*256 */ 486*25c28e83SPiotr Jasiukajtis 0x4087154000000000ULL, /* KA1_HI = 2.8853759765625e+00*256 */ 487*25c28e83SPiotr Jasiukajtis 0x40871547652b82feULL, /* KA1 = 2.885390081777926774e+00*256 */ 488*25c28e83SPiotr Jasiukajtis 0x4110000000000000ULL, /* HTHRESH = 262144.0 */ 489*25c28e83SPiotr Jasiukajtis 0xc110cc0000000000ULL, /* LTHRESH = -275200.0 */ 490*25c28e83SPiotr Jasiukajtis 0x3cd5d52893bc7fecULL, /* KB5 = 1.21195555854068860923e-15 */ 491*25c28e83SPiotr Jasiukajtis 0x3d83b2abc07c93d0ULL, /* KB4 = 2.23939573811855104311e-12 */ 492*25c28e83SPiotr Jasiukajtis 0x3e2c6b08d71f5d1eULL, /* KB3 = 3.30830268126604677436e-09 */ 493*25c28e83SPiotr Jasiukajtis 0x3ecebfbdff82c4edULL, /* KB2 = 3.66556559691003767877e-06 */ 494*25c28e83SPiotr Jasiukajtis 0x3f662e42fefa39efULL, /* KB1 = 2.70760617406228636578e-03 */ 495*25c28e83SPiotr Jasiukajtis 0x01a56e1fc2f8f359ULL, /* _TINY = 1.0e-300 */ 496*25c28e83SPiotr Jasiukajtis 0x7e37e43c8800759cULL /* _HUGE = 1.0e+300 */ 497*25c28e83SPiotr Jasiukajtis }; 498*25c28e83SPiotr Jasiukajtis 499*25c28e83SPiotr Jasiukajtis #define SCALE_ARR ((double*)LCONST + 1) 500*25c28e83SPiotr Jasiukajtis #define _TINY ((double*)LCONST)[20] /* 1.0e-300 */ 501*25c28e83SPiotr Jasiukajtis #define _HUGE ((double*)LCONST)[21] /* 1.0e+300 */ 502*25c28e83SPiotr Jasiukajtis 503*25c28e83SPiotr Jasiukajtis #define RET_SC(I) \ 504*25c28e83SPiotr Jasiukajtis px += stridex; \ 505*25c28e83SPiotr Jasiukajtis py += stridey; \ 506*25c28e83SPiotr Jasiukajtis pz += stridez; \ 507*25c28e83SPiotr Jasiukajtis if (--n <= 0) \ 508*25c28e83SPiotr Jasiukajtis break; \ 509*25c28e83SPiotr Jasiukajtis goto start##I; 510*25c28e83SPiotr Jasiukajtis 511*25c28e83SPiotr Jasiukajtis #define RETURN(I, ret) \ 512*25c28e83SPiotr Jasiukajtis { \ 513*25c28e83SPiotr Jasiukajtis pz[0] = (ret); \ 514*25c28e83SPiotr Jasiukajtis RET_SC(I) \ 515*25c28e83SPiotr Jasiukajtis } 516*25c28e83SPiotr Jasiukajtis 517*25c28e83SPiotr Jasiukajtis #define PREP(I) \ 518*25c28e83SPiotr Jasiukajtis hx = HI(px); \ 519*25c28e83SPiotr Jasiukajtis lx = LO(px); \ 520*25c28e83SPiotr Jasiukajtis hy = HI(py); \ 521*25c28e83SPiotr Jasiukajtis ly = LO(py); \ 522*25c28e83SPiotr Jasiukajtis sx = hx >> 31; \ 523*25c28e83SPiotr Jasiukajtis sy = hy >> 31; \ 524*25c28e83SPiotr Jasiukajtis hx &= 0x7fffffff; \ 525*25c28e83SPiotr Jasiukajtis hy &= 0x7fffffff; \ 526*25c28e83SPiotr Jasiukajtis ull_y0 = *(unsigned long long*)px; \ 527*25c28e83SPiotr Jasiukajtis \ 528*25c28e83SPiotr Jasiukajtis if (hy < 0x3bf00000) /* |Y| < 2^(-64) */ \ 529*25c28e83SPiotr Jasiukajtis { \ 530*25c28e83SPiotr Jasiukajtis y0 = *px; \ 531*25c28e83SPiotr Jasiukajtis if ((hy | ly) == 0) /* pow(X,0) */ \ 532*25c28e83SPiotr Jasiukajtis RETURN (I, DONE) \ 533*25c28e83SPiotr Jasiukajtis if (hx > 0x7ff00000 || (hx == 0x7ff00000 && lx != 0)) /* |X| = Nan */ \ 534*25c28e83SPiotr Jasiukajtis *pz = y0 + y0; \ 535*25c28e83SPiotr Jasiukajtis else if ((hx | lx) == 0 || (hx == 0x7ff00000 && lx == 0)) /* X = 0 or Inf */ \ 536*25c28e83SPiotr Jasiukajtis { \ 537*25c28e83SPiotr Jasiukajtis HI(pz) = hx; \ 538*25c28e83SPiotr Jasiukajtis LO(pz) = lx; \ 539*25c28e83SPiotr Jasiukajtis if (sy) \ 540*25c28e83SPiotr Jasiukajtis *pz = DONE / *pz; \ 541*25c28e83SPiotr Jasiukajtis } \ 542*25c28e83SPiotr Jasiukajtis else \ 543*25c28e83SPiotr Jasiukajtis *pz = (sx) ? DZERO / DZERO : DONE; \ 544*25c28e83SPiotr Jasiukajtis RET_SC(I) \ 545*25c28e83SPiotr Jasiukajtis } \ 546*25c28e83SPiotr Jasiukajtis yisint##I = 0; /* Y - non-integer */ \ 547*25c28e83SPiotr Jasiukajtis exp = hy >> 20; /* Y exponent */ \ 548*25c28e83SPiotr Jasiukajtis ull_y0 &= LMMANT; \ 549*25c28e83SPiotr Jasiukajtis ull_x##I = (ull_y0 | LDONE); \ 550*25c28e83SPiotr Jasiukajtis x##I = *(double*)&ull_x##I; \ 551*25c28e83SPiotr Jasiukajtis ull_ax##I = ((ull_x##I + LMROUND) & LMHI20); \ 552*25c28e83SPiotr Jasiukajtis ax##I = *(double*)&ull_ax##I; \ 553*25c28e83SPiotr Jasiukajtis if (hx >= 0x7ff00000 || exp >= 0x43e) /* X=Inf,Nan or |Y|>2^63,Inf,Nan */ \ 554*25c28e83SPiotr Jasiukajtis { \ 555*25c28e83SPiotr Jasiukajtis y0 = *px; \ 556*25c28e83SPiotr Jasiukajtis if (hx > 0x7ff00000 || (hx == 0x7ff00000 && lx != 0) || \ 557*25c28e83SPiotr Jasiukajtis hy > 0x7ff00000 || (hy == 0x7ff00000 && ly != 0)) /* |X| or |Y| = Nan */ \ 558*25c28e83SPiotr Jasiukajtis RETURN (I, y0 + *py) \ 559*25c28e83SPiotr Jasiukajtis if (hy == 0x7ff00000 && (ly == 0)) /* |Y| = Inf */ \ 560*25c28e83SPiotr Jasiukajtis { \ 561*25c28e83SPiotr Jasiukajtis if (hx == 0x3ff00000 && (lx == 0)) /* +-1 ** +-Inf */ \ 562*25c28e83SPiotr Jasiukajtis *pz = *py - *py; \ 563*25c28e83SPiotr Jasiukajtis else if ((hx < 0x3ff00000) != sy) \ 564*25c28e83SPiotr Jasiukajtis *pz = DZERO; \ 565*25c28e83SPiotr Jasiukajtis else \ 566*25c28e83SPiotr Jasiukajtis { \ 567*25c28e83SPiotr Jasiukajtis HI(pz) = hy; \ 568*25c28e83SPiotr Jasiukajtis LO(pz) = ly; \ 569*25c28e83SPiotr Jasiukajtis } \ 570*25c28e83SPiotr Jasiukajtis RET_SC(I) \ 571*25c28e83SPiotr Jasiukajtis } \ 572*25c28e83SPiotr Jasiukajtis if (exp < 0x43e) /* |Y| < 2^63 */ \ 573*25c28e83SPiotr Jasiukajtis { \ 574*25c28e83SPiotr Jasiukajtis if (sx) /* X = -Inf */ \ 575*25c28e83SPiotr Jasiukajtis { \ 576*25c28e83SPiotr Jasiukajtis if (exp >= 0x434) /* |Y| >= 2^53 */ \ 577*25c28e83SPiotr Jasiukajtis yisint##I = 2; /* Y - even */ \ 578*25c28e83SPiotr Jasiukajtis else \ 579*25c28e83SPiotr Jasiukajtis { \ 580*25c28e83SPiotr Jasiukajtis if (exp >= 0x3ff) /* |Y| >= 1 */ \ 581*25c28e83SPiotr Jasiukajtis { \ 582*25c28e83SPiotr Jasiukajtis if (exp > (20 + 0x3ff)) \ 583*25c28e83SPiotr Jasiukajtis { \ 584*25c28e83SPiotr Jasiukajtis i0 = ly >> (52 - (exp - 0x3ff)); \ 585*25c28e83SPiotr Jasiukajtis if ((i0 << (52 - (exp - 0x3ff))) == ly) \ 586*25c28e83SPiotr Jasiukajtis yisint##I = 2 - (i0 & 1); \ 587*25c28e83SPiotr Jasiukajtis } \ 588*25c28e83SPiotr Jasiukajtis else if (ly == 0) \ 589*25c28e83SPiotr Jasiukajtis { \ 590*25c28e83SPiotr Jasiukajtis i0 = hy >> (20 - (exp - 0x3ff)); \ 591*25c28e83SPiotr Jasiukajtis if ((i0 << (20 - (exp - 0x3ff))) == hy) \ 592*25c28e83SPiotr Jasiukajtis yisint##I = 2 - (i0 & 1); \ 593*25c28e83SPiotr Jasiukajtis } \ 594*25c28e83SPiotr Jasiukajtis } \ 595*25c28e83SPiotr Jasiukajtis } \ 596*25c28e83SPiotr Jasiukajtis } \ 597*25c28e83SPiotr Jasiukajtis if (sy) \ 598*25c28e83SPiotr Jasiukajtis hx = lx = 0; \ 599*25c28e83SPiotr Jasiukajtis hx += yisint##I << 31; \ 600*25c28e83SPiotr Jasiukajtis HI(pz) = hx; \ 601*25c28e83SPiotr Jasiukajtis LO(pz) = lx; \ 602*25c28e83SPiotr Jasiukajtis RET_SC(I) \ 603*25c28e83SPiotr Jasiukajtis } \ 604*25c28e83SPiotr Jasiukajtis else /* |Y| >= 2^63 */ \ 605*25c28e83SPiotr Jasiukajtis { \ 606*25c28e83SPiotr Jasiukajtis /* |X| = 0, 1, Inf */ \ 607*25c28e83SPiotr Jasiukajtis if (lx == 0 && (hx == 0 || hx == 0x3ff00000 || hx == 0x7ff00000)) \ 608*25c28e83SPiotr Jasiukajtis { \ 609*25c28e83SPiotr Jasiukajtis HI(pz) = hx; \ 610*25c28e83SPiotr Jasiukajtis LO(pz) = lx; \ 611*25c28e83SPiotr Jasiukajtis if (sy) \ 612*25c28e83SPiotr Jasiukajtis *pz = DONE / *pz; \ 613*25c28e83SPiotr Jasiukajtis } \ 614*25c28e83SPiotr Jasiukajtis else \ 615*25c28e83SPiotr Jasiukajtis { \ 616*25c28e83SPiotr Jasiukajtis y0 = ((hx < 0x3ff00000) != sy) ? _TINY : _HUGE; \ 617*25c28e83SPiotr Jasiukajtis *pz = y0 * y0; \ 618*25c28e83SPiotr Jasiukajtis } \ 619*25c28e83SPiotr Jasiukajtis RET_SC(I) \ 620*25c28e83SPiotr Jasiukajtis } \ 621*25c28e83SPiotr Jasiukajtis } \ 622*25c28e83SPiotr Jasiukajtis if ((sx || (hx | lx)) == 0) /* X <= 0 */ \ 623*25c28e83SPiotr Jasiukajtis { \ 624*25c28e83SPiotr Jasiukajtis if (exp >= 0x434) /* |Y| >= 2^53 */ \ 625*25c28e83SPiotr Jasiukajtis yisint##I = 2; /* Y - even */ \ 626*25c28e83SPiotr Jasiukajtis else \ 627*25c28e83SPiotr Jasiukajtis { \ 628*25c28e83SPiotr Jasiukajtis if (exp >= 0x3ff) /* |Y| >= 1 */ \ 629*25c28e83SPiotr Jasiukajtis { \ 630*25c28e83SPiotr Jasiukajtis if (exp > (20 + 0x3ff)) \ 631*25c28e83SPiotr Jasiukajtis { \ 632*25c28e83SPiotr Jasiukajtis i0 = ly >> (52 - (exp - 0x3ff)); \ 633*25c28e83SPiotr Jasiukajtis if ((i0 << (52 - (exp - 0x3ff))) == ly) \ 634*25c28e83SPiotr Jasiukajtis yisint##I = 2 - (i0 & 1); \ 635*25c28e83SPiotr Jasiukajtis } \ 636*25c28e83SPiotr Jasiukajtis else if (ly == 0) \ 637*25c28e83SPiotr Jasiukajtis { \ 638*25c28e83SPiotr Jasiukajtis i0 = hy >> (20 - (exp - 0x3ff)); \ 639*25c28e83SPiotr Jasiukajtis if ((i0 << (20 - (exp - 0x3ff))) == hy) \ 640*25c28e83SPiotr Jasiukajtis yisint##I = 2 - (i0 & 1); \ 641*25c28e83SPiotr Jasiukajtis } \ 642*25c28e83SPiotr Jasiukajtis } \ 643*25c28e83SPiotr Jasiukajtis } \ 644*25c28e83SPiotr Jasiukajtis if ((hx | lx) == 0) /* X == 0 */ \ 645*25c28e83SPiotr Jasiukajtis { \ 646*25c28e83SPiotr Jasiukajtis y0 = DZERO; \ 647*25c28e83SPiotr Jasiukajtis if (sy) \ 648*25c28e83SPiotr Jasiukajtis y0 = DONE / y0; \ 649*25c28e83SPiotr Jasiukajtis if (sx & yisint##I) \ 650*25c28e83SPiotr Jasiukajtis y0 = -y0; \ 651*25c28e83SPiotr Jasiukajtis RETURN (I, y0) \ 652*25c28e83SPiotr Jasiukajtis } \ 653*25c28e83SPiotr Jasiukajtis if (yisint##I == 0) /* pow(neg,non-integer) */ \ 654*25c28e83SPiotr Jasiukajtis RETURN (I, DZERO / DZERO) /* NaN */ \ 655*25c28e83SPiotr Jasiukajtis } \ 656*25c28e83SPiotr Jasiukajtis exp = (hx >> 20); \ 657*25c28e83SPiotr Jasiukajtis exp##I = exp - 2046; \ 658*25c28e83SPiotr Jasiukajtis py##I = py; \ 659*25c28e83SPiotr Jasiukajtis pz##I = pz; \ 660*25c28e83SPiotr Jasiukajtis ux##I = x##I + ax##I; \ 661*25c28e83SPiotr Jasiukajtis if (!exp) \ 662*25c28e83SPiotr Jasiukajtis { \ 663*25c28e83SPiotr Jasiukajtis ax##I = (double) ull_y0; \ 664*25c28e83SPiotr Jasiukajtis ull_ax##I = *(unsigned long long*)&ax##I; \ 665*25c28e83SPiotr Jasiukajtis ull_x##I = ((ull_ax##I & LMMANT) | LDONE); \ 666*25c28e83SPiotr Jasiukajtis x##I = *(double*)&ull_x##I; \ 667*25c28e83SPiotr Jasiukajtis exp##I = ((unsigned int*) & ull_ax##I)[0]; \ 668*25c28e83SPiotr Jasiukajtis exp##I = (exp##I >> 20) - (2046 + 1023 + 51); \ 669*25c28e83SPiotr Jasiukajtis ull_ax##I = (ull_x##I + (LMROUND & LMHI20)); \ 670*25c28e83SPiotr Jasiukajtis ax##I = *(double*)&ull_ax##I; \ 671*25c28e83SPiotr Jasiukajtis ux##I = x##I + ax##I; \ 672*25c28e83SPiotr Jasiukajtis } \ 673*25c28e83SPiotr Jasiukajtis ull_x##I = *(unsigned long long *)&ux##I; \ 674*25c28e83SPiotr Jasiukajtis hx##I = HI(&ull_ax##I); \ 675*25c28e83SPiotr Jasiukajtis yd##I = DONE / ux##I; 676*25c28e83SPiotr Jasiukajtis 677*25c28e83SPiotr Jasiukajtis void 678*25c28e83SPiotr Jasiukajtis __vpow(int n, double * restrict px, int stridex, double * restrict py, 679*25c28e83SPiotr Jasiukajtis int stridey, double * restrict pz, int stridez) 680*25c28e83SPiotr Jasiukajtis { 681*25c28e83SPiotr Jasiukajtis double *py0 = 0, *py1 = 0, *py2; 682*25c28e83SPiotr Jasiukajtis double *pz0 = 0, *pz1 = 0, *pz2; 683*25c28e83SPiotr Jasiukajtis double y0, yd0 = 0.0L, u0, s0, s_l0, m_h0; 684*25c28e83SPiotr Jasiukajtis double y1, yd1 = 0.0L, u1, s1, s_l1, m_h1; 685*25c28e83SPiotr Jasiukajtis double y2, yd2, u2, s2, s_l2, m_h2; 686*25c28e83SPiotr Jasiukajtis double ax0 = 0.0L, x0 = 0.0L, s_h0, ux0; 687*25c28e83SPiotr Jasiukajtis double ax1 = 0.0L, x1 = 0.0L, s_h1, ux1; 688*25c28e83SPiotr Jasiukajtis double ax2, x2, s_h2, ux2; 689*25c28e83SPiotr Jasiukajtis int eflag0, gflag0, ind0, i0; 690*25c28e83SPiotr Jasiukajtis int eflag1, gflag1, ind1, i1; 691*25c28e83SPiotr Jasiukajtis int eflag2, gflag2, ind2, i2; 692*25c28e83SPiotr Jasiukajtis int hx0 = 0, yisint0 = 0, exp0 = 0; 693*25c28e83SPiotr Jasiukajtis int hx1 = 0, yisint1 = 0, exp1 = 0; 694*25c28e83SPiotr Jasiukajtis int hx2, yisint2, exp2; 695*25c28e83SPiotr Jasiukajtis int exp, i = 0; 696*25c28e83SPiotr Jasiukajtis unsigned hx, lx, sx, hy, ly, sy; 697*25c28e83SPiotr Jasiukajtis unsigned long long ull_y0, ull_x0, ull_x1, ull_x2, ull_ax0, ull_ax1, ull_ax2; 698*25c28e83SPiotr Jasiukajtis unsigned long long LDONE = ((unsigned long long*)LCONST)[1]; /* 1.0 */ 699*25c28e83SPiotr Jasiukajtis unsigned long long LMMANT = ((unsigned long long*)LCONST)[4]; /* 0x000fffffffffffff */ 700*25c28e83SPiotr Jasiukajtis unsigned long long LMROUND = ((unsigned long long*)LCONST)[5]; /* 0x0000080000000000 */ 701*25c28e83SPiotr Jasiukajtis unsigned long long LMHI20 = ((unsigned long long*)LCONST)[6]; /* 0xfffff00000000000 */ 702*25c28e83SPiotr Jasiukajtis double DONE = ((double*)LCONST)[1]; /* 1.0 */ 703*25c28e83SPiotr Jasiukajtis double DZERO = ((double*)LCONST)[7]; /* 0.0 */ 704*25c28e83SPiotr Jasiukajtis double KA5 = ((double*)LCONST)[8]; /* 5.77078604860893737986e-01*256 */ 705*25c28e83SPiotr Jasiukajtis double KA3 = ((double*)LCONST)[9]; /* 9.61796693925765549423e-01*256 */ 706*25c28e83SPiotr Jasiukajtis double KA1_LO = ((double*)LCONST)[10]; /* 1.41052154268147309568e-05*256 */ 707*25c28e83SPiotr Jasiukajtis double KA1_HI = ((double*)LCONST)[11]; /* 2.8853759765625e+00*256 */ 708*25c28e83SPiotr Jasiukajtis double KA1 = ((double*)LCONST)[12]; /* 2.885390081777926774e+00*256 */ 709*25c28e83SPiotr Jasiukajtis double HTHRESH = ((double*)LCONST)[13]; /* 262144.0 */ 710*25c28e83SPiotr Jasiukajtis double LTHRESH = ((double*)LCONST)[14]; /* -275200.0 */ 711*25c28e83SPiotr Jasiukajtis double KB5 = ((double*)LCONST)[15]; /* 1.21195555854068860923e-15 */ 712*25c28e83SPiotr Jasiukajtis double KB4 = ((double*)LCONST)[16]; /* 2.23939573811855104311e-12 */ 713*25c28e83SPiotr Jasiukajtis double KB3 = ((double*)LCONST)[17]; /* 3.30830268126604677436e-09 */ 714*25c28e83SPiotr Jasiukajtis double KB2 = ((double*)LCONST)[18]; /* 3.66556559691003767877e-06 */ 715*25c28e83SPiotr Jasiukajtis double KB1 = ((double*)LCONST)[19]; /* 2.70760617406228636578e-03 */ 716*25c28e83SPiotr Jasiukajtis 717*25c28e83SPiotr Jasiukajtis if (stridex == 0) 718*25c28e83SPiotr Jasiukajtis { 719*25c28e83SPiotr Jasiukajtis unsigned hx = HI(px); 720*25c28e83SPiotr Jasiukajtis unsigned lx = LO(px); 721*25c28e83SPiotr Jasiukajtis 722*25c28e83SPiotr Jasiukajtis /* if x is a positive normal number not equal to one, 723*25c28e83SPiotr Jasiukajtis call __vpowx */ 724*25c28e83SPiotr Jasiukajtis if (hx >= 0x00100000 && hx < 0x7ff00000 && 725*25c28e83SPiotr Jasiukajtis (hx != 0x3ff00000 || lx != 0)) 726*25c28e83SPiotr Jasiukajtis { 727*25c28e83SPiotr Jasiukajtis __vpowx(n, px, py, stridey, pz, stridez); 728*25c28e83SPiotr Jasiukajtis return; 729*25c28e83SPiotr Jasiukajtis } 730*25c28e83SPiotr Jasiukajtis } 731*25c28e83SPiotr Jasiukajtis 732*25c28e83SPiotr Jasiukajtis do 733*25c28e83SPiotr Jasiukajtis { 734*25c28e83SPiotr Jasiukajtis /* perform si + ydi = 256*log2(xi)*yi */ 735*25c28e83SPiotr Jasiukajtis start0: 736*25c28e83SPiotr Jasiukajtis PREP(0) 737*25c28e83SPiotr Jasiukajtis px += stridex; 738*25c28e83SPiotr Jasiukajtis py += stridey; 739*25c28e83SPiotr Jasiukajtis pz += stridez; 740*25c28e83SPiotr Jasiukajtis i = 1; 741*25c28e83SPiotr Jasiukajtis if (--n <= 0) 742*25c28e83SPiotr Jasiukajtis break; 743*25c28e83SPiotr Jasiukajtis 744*25c28e83SPiotr Jasiukajtis start1: 745*25c28e83SPiotr Jasiukajtis PREP(1) 746*25c28e83SPiotr Jasiukajtis px += stridex; 747*25c28e83SPiotr Jasiukajtis py += stridey; 748*25c28e83SPiotr Jasiukajtis pz += stridez; 749*25c28e83SPiotr Jasiukajtis i = 2; 750*25c28e83SPiotr Jasiukajtis if (--n <= 0) 751*25c28e83SPiotr Jasiukajtis break; 752*25c28e83SPiotr Jasiukajtis 753*25c28e83SPiotr Jasiukajtis start2: 754*25c28e83SPiotr Jasiukajtis PREP(2) 755*25c28e83SPiotr Jasiukajtis 756*25c28e83SPiotr Jasiukajtis u0 = x0 - ax0; 757*25c28e83SPiotr Jasiukajtis u1 = x1 - ax1; 758*25c28e83SPiotr Jasiukajtis u2 = x2 - ax2; 759*25c28e83SPiotr Jasiukajtis 760*25c28e83SPiotr Jasiukajtis s0 = u0 * yd0; 761*25c28e83SPiotr Jasiukajtis LO(&ux0) = 0; 762*25c28e83SPiotr Jasiukajtis s1 = u1 * yd1; 763*25c28e83SPiotr Jasiukajtis LO(&ux1) = 0; 764*25c28e83SPiotr Jasiukajtis s2 = u2 * yd2; 765*25c28e83SPiotr Jasiukajtis LO(&ux2) = 0; 766*25c28e83SPiotr Jasiukajtis 767*25c28e83SPiotr Jasiukajtis y0 = s0 * s0; 768*25c28e83SPiotr Jasiukajtis s_h0 = s0; 769*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 770*25c28e83SPiotr Jasiukajtis y1 = s1 * s1; 771*25c28e83SPiotr Jasiukajtis s_h1 = s1; 772*25c28e83SPiotr Jasiukajtis LO(&s_h1) = 0; 773*25c28e83SPiotr Jasiukajtis y2 = s2 * s2; 774*25c28e83SPiotr Jasiukajtis s_h2 = s2; 775*25c28e83SPiotr Jasiukajtis LO(&s_h2) = 0; 776*25c28e83SPiotr Jasiukajtis 777*25c28e83SPiotr Jasiukajtis s0 = (KA5 * y0 + KA3) * y0 * s0; 778*25c28e83SPiotr Jasiukajtis s1 = (KA5 * y1 + KA3) * y1 * s1; 779*25c28e83SPiotr Jasiukajtis s2 = (KA5 * y2 + KA3) * y2 * s2; 780*25c28e83SPiotr Jasiukajtis 781*25c28e83SPiotr Jasiukajtis s_l0 = (x0 - (ux0 - ax0)); 782*25c28e83SPiotr Jasiukajtis s_l1 = (x1 - (ux1 - ax1)); 783*25c28e83SPiotr Jasiukajtis s_l2 = (x2 - (ux2 - ax2)); 784*25c28e83SPiotr Jasiukajtis 785*25c28e83SPiotr Jasiukajtis s_l0 = u0 - s_h0 * ux0 - s_h0 * s_l0; 786*25c28e83SPiotr Jasiukajtis s_l1 = u1 - s_h1 * ux1 - s_h1 * s_l1; 787*25c28e83SPiotr Jasiukajtis s_l2 = u2 - s_h2 * ux2 - s_h2 * s_l2; 788*25c28e83SPiotr Jasiukajtis 789*25c28e83SPiotr Jasiukajtis s_l0 = KA1 * yd0 * s_l0; 790*25c28e83SPiotr Jasiukajtis i0 = (hx0 >> 8) & 0xff0; 791*25c28e83SPiotr Jasiukajtis exp0 += (hx0 >> 20); 792*25c28e83SPiotr Jasiukajtis 793*25c28e83SPiotr Jasiukajtis s_l1 = KA1 * yd1 * s_l1; 794*25c28e83SPiotr Jasiukajtis i1 = (hx1 >> 8) & 0xff0; 795*25c28e83SPiotr Jasiukajtis exp1 += (hx1 >> 20); 796*25c28e83SPiotr Jasiukajtis 797*25c28e83SPiotr Jasiukajtis s_l2 = KA1 * yd2 * s_l2; 798*25c28e83SPiotr Jasiukajtis i2 = (hx2 >> 8) & 0xff0; 799*25c28e83SPiotr Jasiukajtis exp2 += (hx2 >> 20); 800*25c28e83SPiotr Jasiukajtis 801*25c28e83SPiotr Jasiukajtis yd0 = KA1_HI * s_h0; 802*25c28e83SPiotr Jasiukajtis yd1 = KA1_HI * s_h1; 803*25c28e83SPiotr Jasiukajtis yd2 = KA1_HI * s_h2; 804*25c28e83SPiotr Jasiukajtis 805*25c28e83SPiotr Jasiukajtis y0 = *(double *)((char*)__TBL_log2 + i0); 806*25c28e83SPiotr Jasiukajtis y1 = *(double *)((char*)__TBL_log2 + i1); 807*25c28e83SPiotr Jasiukajtis y2 = *(double *)((char*)__TBL_log2 + i2); 808*25c28e83SPiotr Jasiukajtis 809*25c28e83SPiotr Jasiukajtis y0 += (double)(exp0 << 8); 810*25c28e83SPiotr Jasiukajtis y1 += (double)(exp1 << 8); 811*25c28e83SPiotr Jasiukajtis y2 += (double)(exp2 << 8); 812*25c28e83SPiotr Jasiukajtis 813*25c28e83SPiotr Jasiukajtis m_h0 = y0 + yd0; 814*25c28e83SPiotr Jasiukajtis m_h1 = y1 + yd1; 815*25c28e83SPiotr Jasiukajtis m_h2 = y2 + yd2; 816*25c28e83SPiotr Jasiukajtis 817*25c28e83SPiotr Jasiukajtis y0 = s0 - ((m_h0 - y0 - yd0) - s_l0); 818*25c28e83SPiotr Jasiukajtis y1 = s1 - ((m_h1 - y1 - yd1) - s_l1); 819*25c28e83SPiotr Jasiukajtis y2 = s2 - ((m_h2 - y2 - yd2) - s_l2); 820*25c28e83SPiotr Jasiukajtis 821*25c28e83SPiotr Jasiukajtis y0 += *(double *)((char*)__TBL_log2 + i0 + 8) + KA1_LO * s_h0; 822*25c28e83SPiotr Jasiukajtis y1 += *(double *)((char*)__TBL_log2 + i1 + 8) + KA1_LO * s_h1; 823*25c28e83SPiotr Jasiukajtis y2 += *(double *)((char*)__TBL_log2 + i2 + 8) + KA1_LO * s_h2; 824*25c28e83SPiotr Jasiukajtis 825*25c28e83SPiotr Jasiukajtis s_h0 = y0 + m_h0; 826*25c28e83SPiotr Jasiukajtis s_h1 = y1 + m_h1; 827*25c28e83SPiotr Jasiukajtis s_h2 = y2 + m_h2; 828*25c28e83SPiotr Jasiukajtis 829*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 830*25c28e83SPiotr Jasiukajtis LO(&s_h1) = 0; 831*25c28e83SPiotr Jasiukajtis LO(&s_h2) = 0; 832*25c28e83SPiotr Jasiukajtis 833*25c28e83SPiotr Jasiukajtis yd0 = *py0; 834*25c28e83SPiotr Jasiukajtis yd1 = *py1; 835*25c28e83SPiotr Jasiukajtis yd2 = *py2; 836*25c28e83SPiotr Jasiukajtis s0 = yd0; 837*25c28e83SPiotr Jasiukajtis s1 = yd1; 838*25c28e83SPiotr Jasiukajtis s2 = yd2; 839*25c28e83SPiotr Jasiukajtis LO(&s0) = 0; 840*25c28e83SPiotr Jasiukajtis LO(&s1) = 0; 841*25c28e83SPiotr Jasiukajtis LO(&s2) = 0; 842*25c28e83SPiotr Jasiukajtis 843*25c28e83SPiotr Jasiukajtis y0 = y0 - (s_h0 - m_h0); 844*25c28e83SPiotr Jasiukajtis y1 = y1 - (s_h1 - m_h1); 845*25c28e83SPiotr Jasiukajtis y2 = y2 - (s_h2 - m_h2); 846*25c28e83SPiotr Jasiukajtis 847*25c28e83SPiotr Jasiukajtis yd0 = (yd0 - s0) * s_h0 + yd0 * y0; 848*25c28e83SPiotr Jasiukajtis yd1 = (yd1 - s1) * s_h1 + yd1 * y1; 849*25c28e83SPiotr Jasiukajtis yd2 = (yd2 - s2) * s_h2 + yd2 * y2; 850*25c28e83SPiotr Jasiukajtis 851*25c28e83SPiotr Jasiukajtis s0 = s_h0 * s0; 852*25c28e83SPiotr Jasiukajtis s1 = s_h1 * s1; 853*25c28e83SPiotr Jasiukajtis s2 = s_h2 * s2; 854*25c28e83SPiotr Jasiukajtis 855*25c28e83SPiotr Jasiukajtis /* perform 2 ** ((si+ydi)/256) */ 856*25c28e83SPiotr Jasiukajtis if (s0 > HTHRESH) 857*25c28e83SPiotr Jasiukajtis { 858*25c28e83SPiotr Jasiukajtis s0 = HTHRESH; 859*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 860*25c28e83SPiotr Jasiukajtis } 861*25c28e83SPiotr Jasiukajtis if (s1 > HTHRESH) 862*25c28e83SPiotr Jasiukajtis { 863*25c28e83SPiotr Jasiukajtis s1 = HTHRESH; 864*25c28e83SPiotr Jasiukajtis yd1 = DZERO; 865*25c28e83SPiotr Jasiukajtis } 866*25c28e83SPiotr Jasiukajtis if (s2 > HTHRESH) 867*25c28e83SPiotr Jasiukajtis { 868*25c28e83SPiotr Jasiukajtis s2 = HTHRESH; 869*25c28e83SPiotr Jasiukajtis yd2 = DZERO; 870*25c28e83SPiotr Jasiukajtis } 871*25c28e83SPiotr Jasiukajtis 872*25c28e83SPiotr Jasiukajtis if (s0 < LTHRESH) 873*25c28e83SPiotr Jasiukajtis { 874*25c28e83SPiotr Jasiukajtis s0 = LTHRESH; 875*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 876*25c28e83SPiotr Jasiukajtis } 877*25c28e83SPiotr Jasiukajtis ind0 = (int) (s0 + yd0); 878*25c28e83SPiotr Jasiukajtis if (s1 < LTHRESH) 879*25c28e83SPiotr Jasiukajtis { 880*25c28e83SPiotr Jasiukajtis s1 = LTHRESH; 881*25c28e83SPiotr Jasiukajtis yd1 = DZERO; 882*25c28e83SPiotr Jasiukajtis } 883*25c28e83SPiotr Jasiukajtis ind1 = (int) (s1 + yd1); 884*25c28e83SPiotr Jasiukajtis if (s2 < LTHRESH) 885*25c28e83SPiotr Jasiukajtis { 886*25c28e83SPiotr Jasiukajtis s2 = LTHRESH; 887*25c28e83SPiotr Jasiukajtis yd2 = DZERO; 888*25c28e83SPiotr Jasiukajtis } 889*25c28e83SPiotr Jasiukajtis ind2 = (int) (s2 + yd2); 890*25c28e83SPiotr Jasiukajtis 891*25c28e83SPiotr Jasiukajtis i0 = (ind0 & 0xff) << 4; 892*25c28e83SPiotr Jasiukajtis u0 = (double) ind0; 893*25c28e83SPiotr Jasiukajtis ind0 >>= 8; 894*25c28e83SPiotr Jasiukajtis 895*25c28e83SPiotr Jasiukajtis i1 = (ind1 & 0xff) << 4; 896*25c28e83SPiotr Jasiukajtis u1 = (double)ind1; 897*25c28e83SPiotr Jasiukajtis ind1 >>= 8; 898*25c28e83SPiotr Jasiukajtis 899*25c28e83SPiotr Jasiukajtis i2 = (ind2 & 0xff) << 4; 900*25c28e83SPiotr Jasiukajtis u2 = (double) ind2; 901*25c28e83SPiotr Jasiukajtis ind2 >>= 8; 902*25c28e83SPiotr Jasiukajtis 903*25c28e83SPiotr Jasiukajtis y0 = s0 - u0 + yd0; 904*25c28e83SPiotr Jasiukajtis y1 = s1 - u1 + yd1; 905*25c28e83SPiotr Jasiukajtis y2 = s2 - u2 + yd2; 906*25c28e83SPiotr Jasiukajtis 907*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0); 908*25c28e83SPiotr Jasiukajtis y0 = ((((KB5 * y0 + KB4) * y0 + KB3) * y0 + KB2) * y0 + KB1) * y0; 909*25c28e83SPiotr Jasiukajtis u1 = *(double*)((char*)__TBL_exp2 + i1); 910*25c28e83SPiotr Jasiukajtis y1 = ((((KB5 * y1 + KB4) * y1 + KB3) * y1 + KB2) * y1 + KB1) * y1; 911*25c28e83SPiotr Jasiukajtis u2 = *(double*)((char*)__TBL_exp2 + i2); 912*25c28e83SPiotr Jasiukajtis y2 = ((((KB5 * y2 + KB4) * y2 + KB3) * y2 + KB2) * y2 + KB1) * y2; 913*25c28e83SPiotr Jasiukajtis 914*25c28e83SPiotr Jasiukajtis eflag0 = (ind0 + 1021) >> 31; 915*25c28e83SPiotr Jasiukajtis gflag0 = (1022 - ind0) >> 31; 916*25c28e83SPiotr Jasiukajtis eflag1 = (ind1 + 1021) >> 31; 917*25c28e83SPiotr Jasiukajtis gflag1 = (1022 - ind1) >> 31; 918*25c28e83SPiotr Jasiukajtis eflag2 = (ind2 + 1021) >> 31; 919*25c28e83SPiotr Jasiukajtis gflag2 = (1022 - ind2) >> 31; 920*25c28e83SPiotr Jasiukajtis 921*25c28e83SPiotr Jasiukajtis ind0 = (yisint0 << 11) + ind0 + (54 & eflag0) - (52 & gflag0); 922*25c28e83SPiotr Jasiukajtis ind0 <<= 20; 923*25c28e83SPiotr Jasiukajtis ind1 = (yisint1 << 11) + ind1 + (54 & eflag1) - (52 & gflag1); 924*25c28e83SPiotr Jasiukajtis ind1 <<= 20; 925*25c28e83SPiotr Jasiukajtis ind2 = (yisint2 << 11) + ind2 + (54 & eflag2) - (52 & gflag2); 926*25c28e83SPiotr Jasiukajtis ind2 <<= 20; 927*25c28e83SPiotr Jasiukajtis 928*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0 + 8) + u0 * y0 + u0; 929*25c28e83SPiotr Jasiukajtis u1 = *(double*)((char*)__TBL_exp2 + i1 + 8) + u1 * y1 + u1; 930*25c28e83SPiotr Jasiukajtis u2 = *(double*)((char*)__TBL_exp2 + i2 + 8) + u2 * y2 + u2; 931*25c28e83SPiotr Jasiukajtis 932*25c28e83SPiotr Jasiukajtis ull_x0 = *(unsigned long long*)&u0; 933*25c28e83SPiotr Jasiukajtis HI(&ull_x0) += ind0; 934*25c28e83SPiotr Jasiukajtis u0 = *(double*)&ull_x0; 935*25c28e83SPiotr Jasiukajtis 936*25c28e83SPiotr Jasiukajtis ull_x1 = *(unsigned long long*)&u1; 937*25c28e83SPiotr Jasiukajtis HI(&ull_x1) += ind1; 938*25c28e83SPiotr Jasiukajtis u1 = *(double*)&ull_x1; 939*25c28e83SPiotr Jasiukajtis 940*25c28e83SPiotr Jasiukajtis ull_x2 = *(unsigned long long*)&u2; 941*25c28e83SPiotr Jasiukajtis HI(&ull_x2) += ind2; 942*25c28e83SPiotr Jasiukajtis u2 = *(double*)&ull_x2; 943*25c28e83SPiotr Jasiukajtis 944*25c28e83SPiotr Jasiukajtis *pz0 = u0 * SCALE_ARR[eflag0 - gflag0]; 945*25c28e83SPiotr Jasiukajtis *pz1 = u1 * SCALE_ARR[eflag1 - gflag1]; 946*25c28e83SPiotr Jasiukajtis *pz2 = u2 * SCALE_ARR[eflag2 - gflag2]; 947*25c28e83SPiotr Jasiukajtis 948*25c28e83SPiotr Jasiukajtis px += stridex; 949*25c28e83SPiotr Jasiukajtis py += stridey; 950*25c28e83SPiotr Jasiukajtis pz += stridez; 951*25c28e83SPiotr Jasiukajtis i = 0; 952*25c28e83SPiotr Jasiukajtis 953*25c28e83SPiotr Jasiukajtis } while (--n > 0); 954*25c28e83SPiotr Jasiukajtis 955*25c28e83SPiotr Jasiukajtis if (i > 0) 956*25c28e83SPiotr Jasiukajtis { 957*25c28e83SPiotr Jasiukajtis /* perform si + ydi = 256*log2(xi)*yi */ 958*25c28e83SPiotr Jasiukajtis u0 = x0 - ax0; 959*25c28e83SPiotr Jasiukajtis s0 = u0 * yd0; 960*25c28e83SPiotr Jasiukajtis LO(&ux0) = 0; 961*25c28e83SPiotr Jasiukajtis y0 = s0 * s0; 962*25c28e83SPiotr Jasiukajtis s_h0 = s0; 963*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 964*25c28e83SPiotr Jasiukajtis s0 = (KA5 * y0 + KA3) * y0 * s0; 965*25c28e83SPiotr Jasiukajtis s_l0 = (x0 - (ux0 - ax0)); 966*25c28e83SPiotr Jasiukajtis s_l0 = u0 - s_h0 * ux0 - s_h0 * s_l0; 967*25c28e83SPiotr Jasiukajtis s_l0 = KA1 * yd0 * s_l0; 968*25c28e83SPiotr Jasiukajtis i0 = (hx0 >> 8) & 0xff0; 969*25c28e83SPiotr Jasiukajtis exp0 += (hx0 >> 20); 970*25c28e83SPiotr Jasiukajtis yd0 = KA1_HI * s_h0; 971*25c28e83SPiotr Jasiukajtis y0 = *(double *)((char*)__TBL_log2 + i0); 972*25c28e83SPiotr Jasiukajtis y0 += (double)(exp0 << 8); 973*25c28e83SPiotr Jasiukajtis m_h0 = y0 + yd0; 974*25c28e83SPiotr Jasiukajtis y0 = s0 - ((m_h0 - y0 - yd0) - s_l0); 975*25c28e83SPiotr Jasiukajtis y0 += *(double *)((char*)__TBL_log2 + i0 + 8) + KA1_LO * s_h0; 976*25c28e83SPiotr Jasiukajtis s_h0 = y0 + m_h0; 977*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 978*25c28e83SPiotr Jasiukajtis y0 = y0 - (s_h0 - m_h0); 979*25c28e83SPiotr Jasiukajtis s0 = yd0 = *py0; 980*25c28e83SPiotr Jasiukajtis LO(&s0) = 0; 981*25c28e83SPiotr Jasiukajtis yd0 = (yd0 - s0) * s_h0 + yd0 * y0; 982*25c28e83SPiotr Jasiukajtis s0 = s_h0 * s0; 983*25c28e83SPiotr Jasiukajtis 984*25c28e83SPiotr Jasiukajtis /* perform 2 ** ((si+ydi)/256) */ 985*25c28e83SPiotr Jasiukajtis if (s0 > HTHRESH) 986*25c28e83SPiotr Jasiukajtis { 987*25c28e83SPiotr Jasiukajtis s0 = HTHRESH; 988*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 989*25c28e83SPiotr Jasiukajtis } 990*25c28e83SPiotr Jasiukajtis if (s0 < LTHRESH) 991*25c28e83SPiotr Jasiukajtis { 992*25c28e83SPiotr Jasiukajtis s0 = LTHRESH; 993*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 994*25c28e83SPiotr Jasiukajtis } 995*25c28e83SPiotr Jasiukajtis ind0 = (int) (s0 + yd0); 996*25c28e83SPiotr Jasiukajtis i0 = (ind0 & 0xff) << 4; 997*25c28e83SPiotr Jasiukajtis u0 = (double) ind0; 998*25c28e83SPiotr Jasiukajtis ind0 >>= 8; 999*25c28e83SPiotr Jasiukajtis y0 = s0 - u0 + yd0; 1000*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0); 1001*25c28e83SPiotr Jasiukajtis y0 = ((((KB5 * y0 + KB4) * y0 + KB3) * y0 + KB2) * y0 + KB1) * y0; 1002*25c28e83SPiotr Jasiukajtis eflag0 = (ind0 + 1021) >> 31; 1003*25c28e83SPiotr Jasiukajtis gflag0 = (1022 - ind0) >> 31; 1004*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0 + 8) + u0 * y0 + u0; 1005*25c28e83SPiotr Jasiukajtis ind0 = (yisint0 << 11) + ind0 + (54 & eflag0) - (52 & gflag0); 1006*25c28e83SPiotr Jasiukajtis ind0 <<= 20; 1007*25c28e83SPiotr Jasiukajtis ull_x0 = *(unsigned long long*)&u0; 1008*25c28e83SPiotr Jasiukajtis HI(&ull_x0) += ind0; 1009*25c28e83SPiotr Jasiukajtis u0 = *(double*)&ull_x0; 1010*25c28e83SPiotr Jasiukajtis 1011*25c28e83SPiotr Jasiukajtis *pz0 = u0 * SCALE_ARR[eflag0 - gflag0]; 1012*25c28e83SPiotr Jasiukajtis 1013*25c28e83SPiotr Jasiukajtis if (i > 1) 1014*25c28e83SPiotr Jasiukajtis { 1015*25c28e83SPiotr Jasiukajtis /* perform si + ydi = 256*log2(xi)*yi */ 1016*25c28e83SPiotr Jasiukajtis u0 = x1 - ax1; 1017*25c28e83SPiotr Jasiukajtis s0 = u0 * yd1; 1018*25c28e83SPiotr Jasiukajtis LO(&ux1) = 0; 1019*25c28e83SPiotr Jasiukajtis y0 = s0 * s0; 1020*25c28e83SPiotr Jasiukajtis s_h0 = s0; 1021*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 1022*25c28e83SPiotr Jasiukajtis s0 = (KA5 * y0 + KA3) * y0 * s0; 1023*25c28e83SPiotr Jasiukajtis s_l0 = (x1 - (ux1 - ax1)); 1024*25c28e83SPiotr Jasiukajtis s_l0 = u0 - s_h0 * ux1 - s_h0 * s_l0; 1025*25c28e83SPiotr Jasiukajtis s_l0 = KA1 * yd1 * s_l0; 1026*25c28e83SPiotr Jasiukajtis i0 = (hx1 >> 8) & 0xff0; 1027*25c28e83SPiotr Jasiukajtis exp1 += (hx1 >> 20); 1028*25c28e83SPiotr Jasiukajtis yd0 = KA1_HI * s_h0; 1029*25c28e83SPiotr Jasiukajtis y0 = *(double *)((char*)__TBL_log2 + i0); 1030*25c28e83SPiotr Jasiukajtis y0 += (double)(exp1 << 8); 1031*25c28e83SPiotr Jasiukajtis m_h0 = y0 + yd0; 1032*25c28e83SPiotr Jasiukajtis y0 = s0 - ((m_h0 - y0 - yd0) - s_l0); 1033*25c28e83SPiotr Jasiukajtis y0 += *(double *)((char*)__TBL_log2 + i0 + 8) + KA1_LO * s_h0; 1034*25c28e83SPiotr Jasiukajtis s_h0 = y0 + m_h0; 1035*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 1036*25c28e83SPiotr Jasiukajtis y0 = y0 - (s_h0 - m_h0); 1037*25c28e83SPiotr Jasiukajtis s0 = yd0 = *py1; 1038*25c28e83SPiotr Jasiukajtis LO(&s0) = 0; 1039*25c28e83SPiotr Jasiukajtis yd0 = (yd0 - s0) * s_h0 + yd0 * y0; 1040*25c28e83SPiotr Jasiukajtis s0 = s_h0 * s0; 1041*25c28e83SPiotr Jasiukajtis /* perform 2 ** ((si+ydi)/256) */ 1042*25c28e83SPiotr Jasiukajtis if (s0 > HTHRESH) 1043*25c28e83SPiotr Jasiukajtis { 1044*25c28e83SPiotr Jasiukajtis s0 = HTHRESH; 1045*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1046*25c28e83SPiotr Jasiukajtis } 1047*25c28e83SPiotr Jasiukajtis if (s0 < LTHRESH) 1048*25c28e83SPiotr Jasiukajtis { 1049*25c28e83SPiotr Jasiukajtis s0 = LTHRESH; 1050*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1051*25c28e83SPiotr Jasiukajtis } 1052*25c28e83SPiotr Jasiukajtis ind0 = (int) (s0 + yd0); 1053*25c28e83SPiotr Jasiukajtis i0 = (ind0 & 0xff) << 4; 1054*25c28e83SPiotr Jasiukajtis u0 = (double) ind0; 1055*25c28e83SPiotr Jasiukajtis ind0 >>= 8; 1056*25c28e83SPiotr Jasiukajtis y0 = s0 - u0 + yd0; 1057*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0); 1058*25c28e83SPiotr Jasiukajtis y0 = ((((KB5 * y0 + KB4) * y0 + KB3) * y0 + KB2) * y0 + KB1) * y0; 1059*25c28e83SPiotr Jasiukajtis eflag0 = (ind0 + 1021) >> 31; 1060*25c28e83SPiotr Jasiukajtis gflag0 = (1022 - ind0) >> 31; 1061*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0 + 8) + u0 * y0 + u0; 1062*25c28e83SPiotr Jasiukajtis ind0 = (yisint1 << 11) + ind0 + (54 & eflag0) - (52 & gflag0); 1063*25c28e83SPiotr Jasiukajtis ind0 <<= 20; 1064*25c28e83SPiotr Jasiukajtis ull_x0 = *(unsigned long long*)&u0; 1065*25c28e83SPiotr Jasiukajtis HI(&ull_x0) += ind0; 1066*25c28e83SPiotr Jasiukajtis u0 = *(double*)&ull_x0; 1067*25c28e83SPiotr Jasiukajtis *pz1 = u0 * SCALE_ARR[eflag0 - gflag0]; 1068*25c28e83SPiotr Jasiukajtis } 1069*25c28e83SPiotr Jasiukajtis } 1070*25c28e83SPiotr Jasiukajtis } 1071*25c28e83SPiotr Jasiukajtis 1072*25c28e83SPiotr Jasiukajtis #undef RET_SC 1073*25c28e83SPiotr Jasiukajtis #define RET_SC(I) \ 1074*25c28e83SPiotr Jasiukajtis py += stridey; \ 1075*25c28e83SPiotr Jasiukajtis pz += stridez; \ 1076*25c28e83SPiotr Jasiukajtis if (--n <= 0) \ 1077*25c28e83SPiotr Jasiukajtis break; \ 1078*25c28e83SPiotr Jasiukajtis goto start##I; 1079*25c28e83SPiotr Jasiukajtis 1080*25c28e83SPiotr Jasiukajtis #define PREP_X(I) \ 1081*25c28e83SPiotr Jasiukajtis hy = HI(py); \ 1082*25c28e83SPiotr Jasiukajtis ly = LO(py); \ 1083*25c28e83SPiotr Jasiukajtis sy = hy >> 31; \ 1084*25c28e83SPiotr Jasiukajtis hy &= 0x7fffffff; \ 1085*25c28e83SPiotr Jasiukajtis py##I = py; \ 1086*25c28e83SPiotr Jasiukajtis \ 1087*25c28e83SPiotr Jasiukajtis if (hy < 0x3bf00000) /* |Y| < 2^(-64) */ \ 1088*25c28e83SPiotr Jasiukajtis RETURN (I, DONE) \ 1089*25c28e83SPiotr Jasiukajtis pz##I = pz; \ 1090*25c28e83SPiotr Jasiukajtis if (hy >= 0x43e00000) /* |Y|>2^63,Inf,Nan */ \ 1091*25c28e83SPiotr Jasiukajtis { \ 1092*25c28e83SPiotr Jasiukajtis if (hy >= 0x7ff00000) /* |Y|=Inf,Nan */ \ 1093*25c28e83SPiotr Jasiukajtis { \ 1094*25c28e83SPiotr Jasiukajtis if (hy == 0x7ff00000 && ly == 0) /* |Y|=Inf */ \ 1095*25c28e83SPiotr Jasiukajtis { \ 1096*25c28e83SPiotr Jasiukajtis if ((hx < 0x3ff00000) != sy) \ 1097*25c28e83SPiotr Jasiukajtis *pz = DZERO; \ 1098*25c28e83SPiotr Jasiukajtis else \ 1099*25c28e83SPiotr Jasiukajtis { \ 1100*25c28e83SPiotr Jasiukajtis HI(pz) = hy; \ 1101*25c28e83SPiotr Jasiukajtis LO(pz) = ly; \ 1102*25c28e83SPiotr Jasiukajtis } \ 1103*25c28e83SPiotr Jasiukajtis } \ 1104*25c28e83SPiotr Jasiukajtis else \ 1105*25c28e83SPiotr Jasiukajtis *pz = *px + *py; /* |Y|=Nan */ \ 1106*25c28e83SPiotr Jasiukajtis } \ 1107*25c28e83SPiotr Jasiukajtis else /* |Y|>2^63 */ \ 1108*25c28e83SPiotr Jasiukajtis { \ 1109*25c28e83SPiotr Jasiukajtis y0 = ((hx < 0x3ff00000) != sy) ? _TINY : _HUGE; \ 1110*25c28e83SPiotr Jasiukajtis *pz = y0 * y0; \ 1111*25c28e83SPiotr Jasiukajtis } \ 1112*25c28e83SPiotr Jasiukajtis RET_SC(I) \ 1113*25c28e83SPiotr Jasiukajtis } \ 1114*25c28e83SPiotr Jasiukajtis 1115*25c28e83SPiotr Jasiukajtis #define LMMANT ((unsigned long long*)LCONST)[4] /* 0x000fffffffffffff */ 1116*25c28e83SPiotr Jasiukajtis #define LMROUND ((unsigned long long*)LCONST)[5] /* 0x0000080000000000 */ 1117*25c28e83SPiotr Jasiukajtis #define LMHI20 ((unsigned long long*)LCONST)[6] /* 0xfffff00000000000 */ 1118*25c28e83SPiotr Jasiukajtis #define MMANT ((double*)LCONST)[4] /* 0x000fffffffffffff */ 1119*25c28e83SPiotr Jasiukajtis #define MROUND ((double*)LCONST)[5] /* 0x0000080000000000 */ 1120*25c28e83SPiotr Jasiukajtis #define MHI20 ((double*)LCONST)[6] /* 0xfffff00000000000 */ 1121*25c28e83SPiotr Jasiukajtis #define KA5 ((double*)LCONST)[8] /* 5.77078604860893737986e-01*256 */ 1122*25c28e83SPiotr Jasiukajtis #define KA3 ((double*)LCONST)[9] /* 9.61796693925765549423e-01*256 */ 1123*25c28e83SPiotr Jasiukajtis #define KA1_LO ((double*)LCONST)[10] /* 1.41052154268147309568e-05*256 */ 1124*25c28e83SPiotr Jasiukajtis #define KA1_HI ((double*)LCONST)[11] /* 2.8853759765625e+00*256 */ 1125*25c28e83SPiotr Jasiukajtis #define KA1 ((double*)LCONST)[12] /* 2.885390081777926774e+00*256 */ 1126*25c28e83SPiotr Jasiukajtis 1127*25c28e83SPiotr Jasiukajtis 1128*25c28e83SPiotr Jasiukajtis static void 1129*25c28e83SPiotr Jasiukajtis __vpowx(int n, double * restrict px, double * restrict py, 1130*25c28e83SPiotr Jasiukajtis int stridey, double * restrict pz, int stridez) 1131*25c28e83SPiotr Jasiukajtis { 1132*25c28e83SPiotr Jasiukajtis double *py0, *py1 = 0, *py2; 1133*25c28e83SPiotr Jasiukajtis double *pz0, *pz1 = 0, *pz2; 1134*25c28e83SPiotr Jasiukajtis double ux0, y0, yd0, u0, s0; 1135*25c28e83SPiotr Jasiukajtis double y1, yd1, u1, s1; 1136*25c28e83SPiotr Jasiukajtis double y2, yd2, u2, s2; 1137*25c28e83SPiotr Jasiukajtis double yr, s_h0, s_l0, m_h0, x0, ax0; 1138*25c28e83SPiotr Jasiukajtis unsigned long long ull_y0, ull_x0, ull_x1, ull_x2, ull_ax0; 1139*25c28e83SPiotr Jasiukajtis int eflag0, gflag0, ind0, i0, exp0; 1140*25c28e83SPiotr Jasiukajtis int eflag1, gflag1, ind1, i1; 1141*25c28e83SPiotr Jasiukajtis int eflag2, gflag2, ind2, i2; 1142*25c28e83SPiotr Jasiukajtis int i = 0; 1143*25c28e83SPiotr Jasiukajtis unsigned hx, hx0, hy, ly, sy; 1144*25c28e83SPiotr Jasiukajtis double DONE = ((double*)LCONST)[1]; /* 1.0 */ 1145*25c28e83SPiotr Jasiukajtis unsigned long long LDONE = ((unsigned long long*)LCONST)[1]; /* 1.0 */ 1146*25c28e83SPiotr Jasiukajtis double DZERO = ((double*)LCONST)[7]; /* 0.0 */ 1147*25c28e83SPiotr Jasiukajtis double HTHRESH = ((double*)LCONST)[13]; /* 262144.0 */ 1148*25c28e83SPiotr Jasiukajtis double LTHRESH = ((double*)LCONST)[14]; /* -275200.0 */ 1149*25c28e83SPiotr Jasiukajtis double KB5 = ((double*)LCONST)[15]; /* 1.21195555854068860923e-15 */ 1150*25c28e83SPiotr Jasiukajtis double KB4 = ((double*)LCONST)[16]; /* 2.23939573811855104311e-12 */ 1151*25c28e83SPiotr Jasiukajtis double KB3 = ((double*)LCONST)[17]; /* 3.30830268126604677436e-09 */ 1152*25c28e83SPiotr Jasiukajtis double KB2 = ((double*)LCONST)[18]; /* 3.66556559691003767877e-06 */ 1153*25c28e83SPiotr Jasiukajtis double KB1 = ((double*)LCONST)[19]; /* 2.70760617406228636578e-03 */ 1154*25c28e83SPiotr Jasiukajtis 1155*25c28e83SPiotr Jasiukajtis /* perform s_h + yr = 256*log2(x) */ 1156*25c28e83SPiotr Jasiukajtis ull_y0 = *(unsigned long long*)px; 1157*25c28e83SPiotr Jasiukajtis hx = HI(px); 1158*25c28e83SPiotr Jasiukajtis ull_x0 = (ull_y0 & LMMANT) | LDONE; 1159*25c28e83SPiotr Jasiukajtis x0 = *(double*)&ull_x0; 1160*25c28e83SPiotr Jasiukajtis exp0 = (hx >> 20) - 2046; 1161*25c28e83SPiotr Jasiukajtis ull_ax0 = ull_x0 + (LMROUND & LMHI20); 1162*25c28e83SPiotr Jasiukajtis ax0 = *(double*)&ull_ax0; 1163*25c28e83SPiotr Jasiukajtis hx0 = HI(&ax0); 1164*25c28e83SPiotr Jasiukajtis ux0 = x0 + ax0; 1165*25c28e83SPiotr Jasiukajtis yd0 = DONE / ux0; 1166*25c28e83SPiotr Jasiukajtis u0 = x0 - ax0; 1167*25c28e83SPiotr Jasiukajtis s0 = u0 * yd0; 1168*25c28e83SPiotr Jasiukajtis LO(&ux0) = 0; 1169*25c28e83SPiotr Jasiukajtis y0 = s0 * s0; 1170*25c28e83SPiotr Jasiukajtis s_h0 = s0; 1171*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 1172*25c28e83SPiotr Jasiukajtis s0 = (KA5 * y0 + KA3) * y0 * s0; 1173*25c28e83SPiotr Jasiukajtis s_l0 = (x0 - (ux0 - ax0)); 1174*25c28e83SPiotr Jasiukajtis s_l0 = u0 - s_h0 * ux0 - s_h0 * s_l0; 1175*25c28e83SPiotr Jasiukajtis s_l0 = KA1 * yd0 * s_l0; 1176*25c28e83SPiotr Jasiukajtis i0 = (hx0 >> 8) & 0xff0; 1177*25c28e83SPiotr Jasiukajtis exp0 += (hx0 >> 20); 1178*25c28e83SPiotr Jasiukajtis yd0 = KA1_HI * s_h0; 1179*25c28e83SPiotr Jasiukajtis y0 = *(double *)((char*)__TBL_log2 + i0); 1180*25c28e83SPiotr Jasiukajtis y0 += (double)(exp0 << 8); 1181*25c28e83SPiotr Jasiukajtis m_h0 = y0 + yd0; 1182*25c28e83SPiotr Jasiukajtis y0 = s0 - ((m_h0 - y0 - yd0) - s_l0); 1183*25c28e83SPiotr Jasiukajtis y0 += *(double *)((char*)__TBL_log2 + i0 + 8) + KA1_LO * s_h0; 1184*25c28e83SPiotr Jasiukajtis s_h0 = y0 + m_h0; 1185*25c28e83SPiotr Jasiukajtis LO(&s_h0) = 0; 1186*25c28e83SPiotr Jasiukajtis yr = y0 - (s_h0 - m_h0); 1187*25c28e83SPiotr Jasiukajtis 1188*25c28e83SPiotr Jasiukajtis do 1189*25c28e83SPiotr Jasiukajtis { 1190*25c28e83SPiotr Jasiukajtis /* perform 2 ** ((s_h0+yr)*yi/256) */ 1191*25c28e83SPiotr Jasiukajtis start0: 1192*25c28e83SPiotr Jasiukajtis PREP_X(0) 1193*25c28e83SPiotr Jasiukajtis py += stridey; 1194*25c28e83SPiotr Jasiukajtis pz += stridez; 1195*25c28e83SPiotr Jasiukajtis i = 1; 1196*25c28e83SPiotr Jasiukajtis if (--n <= 0) 1197*25c28e83SPiotr Jasiukajtis break; 1198*25c28e83SPiotr Jasiukajtis 1199*25c28e83SPiotr Jasiukajtis start1: 1200*25c28e83SPiotr Jasiukajtis PREP_X(1) 1201*25c28e83SPiotr Jasiukajtis py += stridey; 1202*25c28e83SPiotr Jasiukajtis pz += stridez; 1203*25c28e83SPiotr Jasiukajtis i = 2; 1204*25c28e83SPiotr Jasiukajtis if (--n <= 0) 1205*25c28e83SPiotr Jasiukajtis break; 1206*25c28e83SPiotr Jasiukajtis 1207*25c28e83SPiotr Jasiukajtis start2: 1208*25c28e83SPiotr Jasiukajtis PREP_X(2) 1209*25c28e83SPiotr Jasiukajtis 1210*25c28e83SPiotr Jasiukajtis s0 = yd0 = *py0; 1211*25c28e83SPiotr Jasiukajtis s1 = yd1 = *py1; 1212*25c28e83SPiotr Jasiukajtis s2 = yd2 = *py2; 1213*25c28e83SPiotr Jasiukajtis 1214*25c28e83SPiotr Jasiukajtis LO(&s0) = 0; 1215*25c28e83SPiotr Jasiukajtis LO(&s1) = 0; 1216*25c28e83SPiotr Jasiukajtis LO(&s2) = 0; 1217*25c28e83SPiotr Jasiukajtis 1218*25c28e83SPiotr Jasiukajtis yd0 = (yd0 - s0) * s_h0 + yd0 * yr; 1219*25c28e83SPiotr Jasiukajtis yd1 = (yd1 - s1) * s_h0 + yd1 * yr; 1220*25c28e83SPiotr Jasiukajtis yd2 = (yd2 - s2) * s_h0 + yd2 * yr; 1221*25c28e83SPiotr Jasiukajtis 1222*25c28e83SPiotr Jasiukajtis s0 = s_h0 * s0; 1223*25c28e83SPiotr Jasiukajtis s1 = s_h0 * s1; 1224*25c28e83SPiotr Jasiukajtis s2 = s_h0 * s2; 1225*25c28e83SPiotr Jasiukajtis 1226*25c28e83SPiotr Jasiukajtis if (s0 > HTHRESH) 1227*25c28e83SPiotr Jasiukajtis { 1228*25c28e83SPiotr Jasiukajtis s0 = HTHRESH; 1229*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1230*25c28e83SPiotr Jasiukajtis } 1231*25c28e83SPiotr Jasiukajtis if (s1 > HTHRESH) 1232*25c28e83SPiotr Jasiukajtis { 1233*25c28e83SPiotr Jasiukajtis s1 = HTHRESH; 1234*25c28e83SPiotr Jasiukajtis yd1 = DZERO; 1235*25c28e83SPiotr Jasiukajtis } 1236*25c28e83SPiotr Jasiukajtis if (s2 > HTHRESH) 1237*25c28e83SPiotr Jasiukajtis { 1238*25c28e83SPiotr Jasiukajtis s2 = HTHRESH; 1239*25c28e83SPiotr Jasiukajtis yd2 = DZERO; 1240*25c28e83SPiotr Jasiukajtis } 1241*25c28e83SPiotr Jasiukajtis 1242*25c28e83SPiotr Jasiukajtis if (s0 < LTHRESH) 1243*25c28e83SPiotr Jasiukajtis { 1244*25c28e83SPiotr Jasiukajtis s0 = LTHRESH; 1245*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1246*25c28e83SPiotr Jasiukajtis } 1247*25c28e83SPiotr Jasiukajtis ind0 = (int) (s0 + yd0); 1248*25c28e83SPiotr Jasiukajtis if (s1 < LTHRESH) 1249*25c28e83SPiotr Jasiukajtis { 1250*25c28e83SPiotr Jasiukajtis s1 = LTHRESH; 1251*25c28e83SPiotr Jasiukajtis yd1 = DZERO; 1252*25c28e83SPiotr Jasiukajtis } 1253*25c28e83SPiotr Jasiukajtis ind1 = (int) (s1 + yd1); 1254*25c28e83SPiotr Jasiukajtis if (s2 < LTHRESH) 1255*25c28e83SPiotr Jasiukajtis { 1256*25c28e83SPiotr Jasiukajtis s2 = LTHRESH; 1257*25c28e83SPiotr Jasiukajtis yd2 = DZERO; 1258*25c28e83SPiotr Jasiukajtis } 1259*25c28e83SPiotr Jasiukajtis ind2 = (int) (s2 + yd2); 1260*25c28e83SPiotr Jasiukajtis 1261*25c28e83SPiotr Jasiukajtis i0 = (ind0 & 0xff) << 4; 1262*25c28e83SPiotr Jasiukajtis u0 = (double) ind0; 1263*25c28e83SPiotr Jasiukajtis ind0 >>= 8; 1264*25c28e83SPiotr Jasiukajtis 1265*25c28e83SPiotr Jasiukajtis i1 = (ind1 & 0xff) << 4; 1266*25c28e83SPiotr Jasiukajtis u1 = (double) ind1; 1267*25c28e83SPiotr Jasiukajtis ind1 >>= 8; 1268*25c28e83SPiotr Jasiukajtis 1269*25c28e83SPiotr Jasiukajtis i2 = (ind2 & 0xff) << 4; 1270*25c28e83SPiotr Jasiukajtis u2 = (double) ind2; 1271*25c28e83SPiotr Jasiukajtis ind2 >>= 8; 1272*25c28e83SPiotr Jasiukajtis 1273*25c28e83SPiotr Jasiukajtis y0 = s0 - u0 + yd0; 1274*25c28e83SPiotr Jasiukajtis y1 = s1 - u1 + yd1; 1275*25c28e83SPiotr Jasiukajtis y2 = s2 - u2 + yd2; 1276*25c28e83SPiotr Jasiukajtis 1277*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0); 1278*25c28e83SPiotr Jasiukajtis y0 = ((((KB5 * y0 + KB4) * y0 + KB3) * y0 + KB2) * y0 + KB1) * y0; 1279*25c28e83SPiotr Jasiukajtis u1 = *(double*)((char*)__TBL_exp2 + i1); 1280*25c28e83SPiotr Jasiukajtis y1 = ((((KB5 * y1 + KB4) * y1 + KB3) * y1 + KB2) * y1 + KB1) * y1; 1281*25c28e83SPiotr Jasiukajtis u2 = *(double*)((char*)__TBL_exp2 + i2); 1282*25c28e83SPiotr Jasiukajtis y2 = ((((KB5 * y2 + KB4) * y2 + KB3) * y2 + KB2) * y2 + KB1) * y2; 1283*25c28e83SPiotr Jasiukajtis 1284*25c28e83SPiotr Jasiukajtis eflag0 = (ind0 + 1021) >> 31; 1285*25c28e83SPiotr Jasiukajtis gflag0 = (1022 - ind0) >> 31; 1286*25c28e83SPiotr Jasiukajtis eflag1 = (ind1 + 1021) >> 31; 1287*25c28e83SPiotr Jasiukajtis gflag1 = (1022 - ind1) >> 31; 1288*25c28e83SPiotr Jasiukajtis eflag2 = (ind2 + 1021) >> 31; 1289*25c28e83SPiotr Jasiukajtis gflag2 = (1022 - ind2) >> 31; 1290*25c28e83SPiotr Jasiukajtis 1291*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0 + 8) + u0 * y0 + u0; 1292*25c28e83SPiotr Jasiukajtis ind0 = ind0 + (54 & eflag0) - (52 & gflag0); 1293*25c28e83SPiotr Jasiukajtis ind0 <<= 20; 1294*25c28e83SPiotr Jasiukajtis ull_x0 = *(unsigned long long*)&u0; 1295*25c28e83SPiotr Jasiukajtis HI(&ull_x0) += ind0; 1296*25c28e83SPiotr Jasiukajtis u0 = *(double*)&ull_x0; 1297*25c28e83SPiotr Jasiukajtis 1298*25c28e83SPiotr Jasiukajtis u1 = *(double*)((char*)__TBL_exp2 + i1 + 8) + u1 * y1 + u1; 1299*25c28e83SPiotr Jasiukajtis ind1 = ind1 + (54 & eflag1) - (52 & gflag1); 1300*25c28e83SPiotr Jasiukajtis ind1 <<= 20; 1301*25c28e83SPiotr Jasiukajtis ull_x1 = *(unsigned long long*)&u1; 1302*25c28e83SPiotr Jasiukajtis HI(&ull_x1) += ind1; 1303*25c28e83SPiotr Jasiukajtis u1 = *(double*)&ull_x1; 1304*25c28e83SPiotr Jasiukajtis 1305*25c28e83SPiotr Jasiukajtis u2 = *(double*)((char*)__TBL_exp2 + i2 + 8) + u2 * y2 + u2; 1306*25c28e83SPiotr Jasiukajtis ind2 = ind2 + (54 & eflag2) - (52 & gflag2); 1307*25c28e83SPiotr Jasiukajtis ind2 <<= 20; 1308*25c28e83SPiotr Jasiukajtis ull_x2 = *(unsigned long long*)&u2; 1309*25c28e83SPiotr Jasiukajtis HI(&ull_x2) += ind2; 1310*25c28e83SPiotr Jasiukajtis u2 = *(double*)&ull_x2; 1311*25c28e83SPiotr Jasiukajtis 1312*25c28e83SPiotr Jasiukajtis *pz0 = u0 * SCALE_ARR[eflag0 - gflag0]; 1313*25c28e83SPiotr Jasiukajtis *pz1 = u1 * SCALE_ARR[eflag1 - gflag1]; 1314*25c28e83SPiotr Jasiukajtis *pz2 = u2 * SCALE_ARR[eflag2 - gflag2]; 1315*25c28e83SPiotr Jasiukajtis 1316*25c28e83SPiotr Jasiukajtis py += stridey; 1317*25c28e83SPiotr Jasiukajtis pz += stridez; 1318*25c28e83SPiotr Jasiukajtis i = 0; 1319*25c28e83SPiotr Jasiukajtis 1320*25c28e83SPiotr Jasiukajtis } while (--n > 0); 1321*25c28e83SPiotr Jasiukajtis 1322*25c28e83SPiotr Jasiukajtis if (i > 0) 1323*25c28e83SPiotr Jasiukajtis { 1324*25c28e83SPiotr Jasiukajtis /* perform 2 ** ((s_h0+yr)*yi/256) */ 1325*25c28e83SPiotr Jasiukajtis s0 = y0 = *py0; 1326*25c28e83SPiotr Jasiukajtis LO(&s0) = 0; 1327*25c28e83SPiotr Jasiukajtis yd0 = (y0 - s0) * s_h0 + y0 * yr; 1328*25c28e83SPiotr Jasiukajtis s0 = s_h0 * s0; 1329*25c28e83SPiotr Jasiukajtis if (s0 > HTHRESH) 1330*25c28e83SPiotr Jasiukajtis { 1331*25c28e83SPiotr Jasiukajtis s0 = HTHRESH; 1332*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1333*25c28e83SPiotr Jasiukajtis } 1334*25c28e83SPiotr Jasiukajtis if (s0 < LTHRESH) 1335*25c28e83SPiotr Jasiukajtis { 1336*25c28e83SPiotr Jasiukajtis s0 = LTHRESH; 1337*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1338*25c28e83SPiotr Jasiukajtis } 1339*25c28e83SPiotr Jasiukajtis ind0 = (int) (s0 + yd0); 1340*25c28e83SPiotr Jasiukajtis i0 = (ind0 & 0xff) << 4; 1341*25c28e83SPiotr Jasiukajtis u0 = (double) ind0; 1342*25c28e83SPiotr Jasiukajtis ind0 >>= 8; 1343*25c28e83SPiotr Jasiukajtis y0 = s0 - u0 + yd0; 1344*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0); 1345*25c28e83SPiotr Jasiukajtis y0 = ((((KB5 * y0 + KB4) * y0 + KB3) * y0 + KB2) * y0 + KB1) * y0; 1346*25c28e83SPiotr Jasiukajtis eflag0 = (ind0 + 1021) >> 31; 1347*25c28e83SPiotr Jasiukajtis gflag0 = (1022 - ind0) >> 31; 1348*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0 + 8) + u0 * y0 + u0; 1349*25c28e83SPiotr Jasiukajtis ind0 = ind0 + (54 & eflag0) - (52 & gflag0); 1350*25c28e83SPiotr Jasiukajtis ind0 <<= 20; 1351*25c28e83SPiotr Jasiukajtis ull_x0 = *(unsigned long long*)&u0; 1352*25c28e83SPiotr Jasiukajtis HI(&ull_x0) += ind0; 1353*25c28e83SPiotr Jasiukajtis u0 = *(double*)&ull_x0; 1354*25c28e83SPiotr Jasiukajtis *pz0 = u0 * SCALE_ARR[eflag0 - gflag0]; 1355*25c28e83SPiotr Jasiukajtis 1356*25c28e83SPiotr Jasiukajtis if (i > 1) 1357*25c28e83SPiotr Jasiukajtis { 1358*25c28e83SPiotr Jasiukajtis /* perform 2 ** ((s_h0+yr)*yi/256) */ 1359*25c28e83SPiotr Jasiukajtis s0 = y0 = *py1; 1360*25c28e83SPiotr Jasiukajtis LO(&s0) = 0; 1361*25c28e83SPiotr Jasiukajtis yd0 = (y0 - s0) * s_h0 + y0 * yr; 1362*25c28e83SPiotr Jasiukajtis s0 = s_h0 * s0; 1363*25c28e83SPiotr Jasiukajtis if (s0 > HTHRESH) 1364*25c28e83SPiotr Jasiukajtis { 1365*25c28e83SPiotr Jasiukajtis s0 = HTHRESH; 1366*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1367*25c28e83SPiotr Jasiukajtis } 1368*25c28e83SPiotr Jasiukajtis if (s0 < LTHRESH) 1369*25c28e83SPiotr Jasiukajtis { 1370*25c28e83SPiotr Jasiukajtis s0 = LTHRESH; 1371*25c28e83SPiotr Jasiukajtis yd0 = DZERO; 1372*25c28e83SPiotr Jasiukajtis } 1373*25c28e83SPiotr Jasiukajtis ind0 = (int) (s0 + yd0); 1374*25c28e83SPiotr Jasiukajtis i0 = (ind0 & 0xff) << 4; 1375*25c28e83SPiotr Jasiukajtis u0 = (double) ind0; 1376*25c28e83SPiotr Jasiukajtis ind0 >>= 8; 1377*25c28e83SPiotr Jasiukajtis y0 = s0 - u0 + yd0; 1378*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0); 1379*25c28e83SPiotr Jasiukajtis y0 = ((((KB5 * y0 + KB4) * y0 + KB3) * y0 + KB2) * y0 + KB1) * y0; 1380*25c28e83SPiotr Jasiukajtis eflag0 = (ind0 + 1021) >> 31; 1381*25c28e83SPiotr Jasiukajtis gflag0 = (1022 - ind0) >> 31; 1382*25c28e83SPiotr Jasiukajtis u0 = *(double*)((char*)__TBL_exp2 + i0 + 8) + u0 * y0 + u0; 1383*25c28e83SPiotr Jasiukajtis ind0 = ind0 + (54 & eflag0) - (52 & gflag0); 1384*25c28e83SPiotr Jasiukajtis ind0 <<= 20; 1385*25c28e83SPiotr Jasiukajtis ull_x0 = *(unsigned long long*)&u0; 1386*25c28e83SPiotr Jasiukajtis HI(&ull_x0) += ind0; 1387*25c28e83SPiotr Jasiukajtis u0 = *(double*)&ull_x0; 1388*25c28e83SPiotr Jasiukajtis *pz1 = u0 * SCALE_ARR[eflag0 - gflag0]; 1389*25c28e83SPiotr Jasiukajtis } 1390*25c28e83SPiotr Jasiukajtis } 1391*25c28e83SPiotr Jasiukajtis } 1392