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