Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/complex/k_clog_rl.c
+++ new/usr/src/lib/libm/common/complex/k_clog_rl.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
|
↓ open down ↓ |
10 lines elided |
↑ open up ↑ |
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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 24 */
25 +
24 26 /*
25 27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 28 * Use is subject to license terms.
27 29 */
28 30
29 -#include "libm.h" /* __k_clog_rl */
31 +#include "libm.h" /* __k_clog_rl */
30 32 #include "complex_wrapper.h"
31 33 #include "longdouble.h"
32 34
33 -/* INDENT OFF */
35 +
34 36 /*
35 37 * long double __k_clog_rl(long double x, long double y, long double *e);
36 38 *
37 39 * Compute real part of complex natural logarithm of x+iy in extra precision
38 40 *
39 41 * __k_clog_rl returns log(hypot(x, y)) with a correction term e.
40 42 *
41 43 * Accuracy: quad 140 bits, intel extended 91 bits.
42 44 *
43 45 * Method.
44 46 * Assume X > Y >= 0 . Let X = 2**nx * x, Y = 2**nx * y, where 1 <= x < 2.
45 47 * Let Z = X*X + Y*Y. Then Z = 2**(nx+nx) * z, where z = x*x + y*y.
46 48 * Note that z < 8.
47 49 * Let Z = x*x + y*y. Z can be normalized as Z = 2**N * z, 1 <= z < 2.
48 50 * We further break down z into 1 + zk + zh + zt, where
49 51 * zk = K*(2**-7) matches z to 7.5 significant bits, 0 <= K <= 2**(-7)-1
50 52 * zh = (z-zk) rounded to half of the current significant bits
51 53 * zt = (z-zk-zh) rounded.
52 54 *
53 55 * z - (1+zk) (zh+zt)
54 56 * Let s = ------------ = ---------------, then
55 57 * z + (1+zk) 2(1+zk)+zh+zt
56 58 * z
57 59 * log(Z) = N*log2 + log(z) = N*log2 + log(1+zk) + log(------)
58 60 * 1+zk
59 61 * 1+s
60 62 * = N * log2 + log(1 +zk) + log(---)
61 63 * 1-s
62 64 *
63 65 * 3 5
64 66 * = N*log2 + log(1+zk) + 2s + 1/12(2s) + 1/80(2s) + ...
65 67 *
66 68 *
|
↓ open down ↓ |
23 lines elided |
↑ open up ↑ |
67 69 * Note 1. For IEEE double precision, a fifteen degree odd polynomial
68 70 * 2s + P1*(2s)^3 + P2*(2s)^5 + P3*(2s)^7 + ... + P7*(2s)^15
69 71 * is generated by a special remez algorithm to
70 72 * approx log((1+s)/(1-s)) accurte to 145 bits.
71 73 * Note 2. 2s can be computed accurately as s2h+s2t by
72 74 * r = 2/((zh+zt)+2(1+zk))
73 75 * s2 = r*(zh+zt)
74 76 * s2h = s2 rounded to double; v = 0.5*s2h;
75 77 * s2t = r*((((zh-s2h*(1+zk))-v*zh)+zt)-v*zt)
76 78 */
77 -/* INDENT ON */
78 79
79 -static const long double
80 -zero = 0.0L,
81 -half = 0.5L,
82 -two = 2.0L,
83 -two240 = 1.7668470647783843295832975007429185158274839e+72L, /* 2^240 */
80 +static const long double zero = 0.0L,
81 + half = 0.5L,
82 + two = 2.0L,
83 + two240 = 1.7668470647783843295832975007429185158274839e+72L; /* 2^240 */
84 84
85 85 /* first 48 bits of ln2 */
86 -ln2_h = 0.693147180559943620892227045260369777679443359375L,
87 -ln2_t = 1.68852500507619780679039605677498525525412068e-15L,
88 -P1 = .083333333333333333333333333333333333341023785768375L,
89 -P2 = .01249999999999999999999999999999679085402075766159375L,
90 -P3 = .002232142857142857142857143310092047621284490564671875L,
91 -P4 = .00043402777777777777774746781319264872413156956512109375L,
92 -P5 = .0000887784090909101756336594019277185263940665468935546875L,
93 -P6 = .000018780048055589639895360927834628371268354778446533203125L,
94 -P7 = .000004069227854328982921366736003458838031087153635406494140625L;
86 +static const long double
87 + ln2_h = 0.693147180559943620892227045260369777679443359375L,
88 + ln2_t = 1.68852500507619780679039605677498525525412068e-15L,
89 + P1 = .083333333333333333333333333333333333341023785768375L,
90 + P2 = .01249999999999999999999999999999679085402075766159375L,
91 + P3 = .002232142857142857142857143310092047621284490564671875L,
92 + P4 = .00043402777777777777774746781319264872413156956512109375L,
93 + P5 = .0000887784090909101756336594019277185263940665468935546875L,
94 + P6 = .000018780048055589639895360927834628371268354778446533203125L,
95 + P7 = .000004069227854328982921366736003458838031087153635406494140625L;
95 96
96 97 /*
97 98 * T[2k, 2k+1] = log(1+k*2**-7) for k = 0, ..., 2**7 - 1,
98 99 * with T[2k] * 2^48 is an int
99 100 */
100 -
101 101 static const long double TBL_log1k[] = {
102 -0.0000000000000000000000000000000000000000e+00L,
103 -0.0000000000000000000000000000000000000000e+00L,
104 -7.7821404420532758194894995540380477905273e-03L,
105 -1.6731279734005070987158875984584325351222e-15L,
106 -1.5504186535963526694104075431823730468750e-02L,
107 -1.7274567499706106231054091184928671990316e-15L,
108 -2.3167059281533397552266251295804977416992e-02L,
109 -9.8067653290966648493916241687661877474892e-16L,
110 -3.0771658666751022792595904320478439331055e-02L,
111 -2.6655784323032762937247606420524589813624e-15L,
112 -3.8318864302134159061097307130694389343262e-02L,
113 -2.4401326580179931029010027013316092332340e-15L,
114 -4.5809536031292452662455616518855094909668e-02L,
115 -1.7505042236510958082472042641283104263139e-15L,
116 -5.3244514518809182845870964229106903076172e-02L,
117 -3.1000199992295574218738634002122149891138e-15L,
118 -6.0624621816433688081815489567816257476807e-02L,
119 -1.1544987906424726040058093958345197512800e-15L,
120 -6.7950661908504628172522643581032752990723e-02L,
121 -3.1212220426341915966610439115772728417386e-15L,
122 -7.5223421237584631171557703055441379547119e-02L,
123 -2.8945270476369282210350897509258766743153e-15L,
124 -8.2443669211073711267090402543544769287109e-02L,
125 -8.8000106966612476303662698634483335676886e-16L,
126 -8.9612158689686083334891009144484996795654e-02L,
127 -1.0492850604602339995319895311151740799226e-15L,
128 -9.6729626458550654888313147239387035369873e-02L,
129 -4.5740725790924807640164516707244620870662e-16L,
130 -1.0379679368164218544734467286616563796997e-01L,
131 -1.3793787171308978090503366050174239822054e-15L,
132 -1.1081436634028918319927470292896032333374e-01L,
133 -9.3099553146639425160476473362380086036919e-16L,
134 -1.1778303565638026384476688690483570098877e-01L,
135 -3.1906940272225656860040797111813146690890e-15L,
136 -1.2470347850095464536934741772711277008057e-01L,
137 -2.5904940590976537504984110469214193890052e-15L,
138 -1.3157635778871679121948545798659324645996e-01L,
139 -2.4813692306707028899159917911012100567219e-15L,
140 -1.3840232285911824305912887211889028549194e-01L,
141 -8.9262619700148275890190121571708972000380e-16L,
142 -1.4518200984449691759436973370611667633057e-01L,
143 -9.7968756533003444764719201050911636480025e-16L,
144 -1.5191604202583874894116888754069805145264e-01L,
145 -3.2261306345373561864598749471119213018106e-15L,
146 -1.5860503017663774016909883357584476470947e-01L,
147 -8.4392427234104999681053621980394827998735e-16L,
148 -1.6524957289530561865831259638071060180664e-01L,
149 -1.5442172988528965297119225948270579746101e-15L,
150 -1.7185025692665689689420105423778295516968e-01L,
151 -2.3254458978918173643097657009894831132739e-15L,
152 -1.7840765747281750464026117697358131408691e-01L,
153 -7.9247913906453736065426776912520942036896e-16L,
154 -1.8492233849401173984006163664162158966064e-01L,
155 -2.5282384195601762803134514624610774126020e-16L,
156 -1.9139485299962899489401024766266345977783e-01L,
157 -4.5971528855989864541366920731297729269228e-16L,
158 -1.9782574332991842425144568551331758499146e-01L,
159 -1.4561111263856836438840838027526567191527e-15L,
160 -2.0421554142868814096800633706152439117432e-01L,
161 -2.7505358140491347148810394262840919337709e-15L,
162 -2.1056476910734645002776233013719320297241e-01L,
163 -3.1876417904825951583107481283088861928977e-15L,
164 -2.1687393830061196808856038842350244522095e-01L,
165 -2.3915305291373208450532580201045871599499e-15L,
166 -2.2314355131420882116799475625157356262207e-01L,
167 -9.3459830033405826094075253077304795996257e-16L,
168 -2.2937410106484534821902343537658452987671e-01L,
169 -4.8177245728966955534167425511952551974164e-16L,
170 -2.3556607131276408040321257431060075759888e-01L,
171 -2.8286743756446304426525380844720043381780e-15L,
172 -2.4171993688714366044223424978554248809814e-01L,
173 -1.5077020732661279714120052415509585052975e-15L,
174 -2.4783616390458007572306087240576744079590e-01L,
175 -1.1810575418933407573072030113600980623171e-15L,
176 -2.5391520998096339667426946107298135757446e-01L,
177 -4.7463053836833625309891834934881898560705e-17L,
178 -2.5995752443692410338371701072901487350464e-01L,
179 -1.9635883624838132961710716735786266795913e-15L,
180 -2.6596354849713677026556979399174451828003e-01L,
181 -1.1710735561325457988709887923652142233351e-15L,
182 -2.7193371548364098089223261922597885131836e-01L,
183 -7.7793943687530702031066421537496360004376e-16L,
184 -2.7786845100345303194444568362087011337280e-01L,
185 -3.2742419043493025311197092322146237692165e-15L,
186 -2.8376817313064250924981024581938982009888e-01L,
187 -2.0890970909765308649465619266075677112425e-15L,
188 -2.8963329258304071345264674164354801177979e-01L,
189 -1.9634262463138821209582240742801727823629e-15L,
190 -2.9546421289383317798638017848134040832520e-01L,
191 -2.6984003017275736237868564402005801750600e-15L,
192 -3.0126133057816062432721082586795091629028e-01L,
193 -1.1566856647123658045763670687640673680383e-15L,
194 -3.0702503529490954292668902780860662460327e-01L,
195 -2.3191484355127267712770857311812090801833e-15L,
196 -3.1275571000389490450288576539605855941772e-01L,
197 -1.9838833607942922604727420618882220398852e-15L,
198 -3.1845373111853447767316538374871015548706e-01L,
199 -1.3813708182984188944010814590398164268227e-16L,
200 -3.2411946865421015218089451082050800323486e-01L,
201 -1.8239097762496144793489474731253815376404e-15L,
202 -3.2975328637246548169059678912162780761719e-01L,
203 -2.5001238260227991620033344720809714552230e-15L,
204 -3.3535554192113536942088103387504816055298e-01L,
205 -2.4608362985459391180385214539620341910962e-15L,
206 -3.4092658697059263772644044365733861923218e-01L,
207 -5.7257864875612301758921090406373771458003e-16L,
208 -3.4646676734620740489845047704875469207764e-01L,
209 -1.1760200117113770182586341947822306069951e-15L,
210 -3.5197642315717558858523261733353137969971e-01L,
211 -2.5960702148389259075462896448369304790506e-15L,
212 -3.5745588892180180096147523727267980575562e-01L,
213 -1.9732645342528682246686790561260072184839e-15L,
214 -3.6290549368936808605212718248367309570312e-01L,
215 -3.6708569716349381675043725477739939978160e-16L,
216 -3.6832556115870573876236448995769023895264e-01L,
217 -1.9142858656640927085879445412821643247628e-15L,
218 -3.7371640979358389245135185774415731430054e-01L,
219 -1.8836966497497166619234389157276681281343e-16L,
220 -3.7907835293496816575498087331652641296387e-01L,
221 -1.2926358724723144934459175417385013725801e-15L,
222 -3.8441169891033055705520382616668939590454e-01L,
223 -1.4826795862363146014726140088145939341729e-15L,
224 -3.8971675114002479745067830663174390792847e-01L,
225 -4.1591978529737177695912258866565331189698e-16L,
226 -3.9499380824086571806219581048935651779175e-01L,
227 -3.2600441982258756252505182317625310732365e-15L,
228 -4.0024316412701210765590076334774494171143e-01L,
229 -5.9927342433864738622836851475469574662703e-16L,
230 -4.0546510810816371872533636633306741714478e-01L,
231 -6.6325267674913128171942721503283748008372e-16L,
232 -4.1065992498526782128465129062533378601074e-01L,
233 -5.6464965491255048900165082436455718077885e-16L,
234 -4.1582789514371043537721561733633279800415e-01L,
235 -5.3023611327561856950735176370587227509442e-16L,
236 -4.2096929464412724541944044176489114761353e-01L,
237 -2.3907094267197419048248363335257046791153e-15L,
238 -4.2608439531089814522601955104619264602661e-01L,
239 -1.9178985253285492839728700574592375309985e-15L,
240 -4.3117346481836804628073878120630979537964e-01L,
241 -3.2945784336977492852031005044499611665595e-15L,
242 -4.3623676677491474151793227065354585647583e-01L,
243 -3.3288311090524075754441878570852962903891e-15L,
244 -4.4127456080487448275562201160937547683716e-01L,
245 -7.4673387443005192574852544613692268411229e-16L,
246 -4.4628710262841764233598951250314712524414e-01L,
247 -1.8691966006681165218815050615460959199251e-15L,
248 -4.5127464413945617138779198285192251205444e-01L,
249 -2.4137569004002270899666314791611479063976e-15L,
250 -4.5623743348158640742440184112638235092163e-01L,
251 -1.1869564036970375473975162509216610120281e-15L,
252 -4.6117571512216670726047595962882041931152e-01L,
253 -3.4591075239659690349392915732654828400811e-15L,
254 -4.6608972992459740680715185590088367462158e-01L,
255 -1.8177514673916038857252366108673570603067e-15L,
256 -4.7097971521878889689105562865734100341797e-01L,
257 -2.1156558422273990182479555421331461933366e-15L,
258 -4.7584590486996347635795245878398418426514e-01L,
259 -4.3790725712752039722791012358345927696967e-16L,
260 -4.8068852934575190261057286988943815231323e-01L,
261 -5.0660455855585733988956280680891477171499e-18L,
262 -4.8550781578169832641833636444061994552612e-01L,
263 -2.4813834547127501689550526444948043590905e-15L,
264 -4.9030398804519137456736643798649311065674e-01L,
265 -2.4635829797216592537498738468934647345741e-15L,
266 -4.9507726679784980206022737547755241394043e-01L,
267 -1.7125377372093652812514167461480115600063e-15L,
268 -4.9982786955644797899367404170334339141846e-01L,
269 -1.3508276573735437007500942002018098437396e-15L,
270 -5.0455601075239187025545106735080480575562e-01L,
271 -3.4168028574643873701242268618467347998876e-15L,
272 -5.0926190178980590417268103919923305511475e-01L,
273 -2.0426313938800290907697638200502614622891e-15L,
274 -5.1394575110223428282552049495279788970947e-01L,
275 -3.3975485593321419703400672813719873194659e-17L,
276 -5.1860776420804555186805373523384332656860e-01L,
277 -8.0284923261130955371987633083003284697416e-17L,
278 -5.2324814376454753528378205373883247375488e-01L,
279 -3.0123302517119603836788558832352723470118e-16L,
280 -5.2786708962084105678513878956437110900879e-01L,
281 -1.3283287534282139298545497336570406582397e-15L,
282 -5.3246479886946929127589100971817970275879e-01L,
283 -2.5525980327137419625398485590148417041921e-15L,
284 -5.3704146589688050994482182431966066360474e-01L,
285 -3.1446219074198341716354190061340477078626e-15L,
286 -5.4159728243274329884116014000028371810913e-01L,
287 -1.0727353821639001503808606766770295812627e-15L,
288 -5.4613243759813556721383065450936555862427e-01L,
289 -8.3168566554721843605240702438699163825794e-17L,
290 -5.5064711795266063631970610003918409347534e-01L,
291 -1.6429402420791657293666192255419538448840e-15L,
292 -5.5514150754050106684189813677221536636353e-01L,
293 -5.2587358222274368868380660194332415847228e-16L,
294 -5.5961578793542088305912329815328121185303e-01L,
295 -1.8032117652023735453816330571171114110385e-15L,
296 -5.6407013828480145889443519990891218185425e-01L,
297 -1.5071769490901812785299634348367857600711e-15L,
298 -5.6850473535266843327917740680277347564697e-01L,
299 -2.7879956135806418878792935692629147550413e-16L,
300 -5.7291975356178426181941176764667034149170e-01L,
301 -1.2472733449589795907271346997596471822345e-15L,
302 -5.7731536503482061561953742057085037231445e-01L,
303 -2.9886985746409486460291929160223207644146e-15L,
304 -5.8169173963462128540413687005639076232910e-01L,
305 -1.1971164738836689815783808674399742176950e-15L,
306 -5.8604904500357690722012193873524665832520e-01L,
307 -1.3016839974975520776911897855504474452726e-15L,
308 -5.9038744660217545856539800297468900680542e-01L,
309 -9.1607651870514890975077236127894522134392e-16L,
310 -5.9470710774668944509357970673590898513794e-01L,
311 -3.3444207638397932963480545729233567201211e-15L,
312 -5.9900818964608149030937056522816419601440e-01L,
313 -1.9090722294592334873060460706130642200729e-15L,
314 -6.0329085143808214297678205184638500213623e-01L,
315 -2.1193638031348149256035110177854940281795e-15L,
316 -6.0755525022453937822319858241826295852661e-01L,
317 -2.4172778865703728624133665395876418941354e-15L,
318 -6.1180154110599005434778518974781036376953e-01L,
319 -2.8491821045766810044199163148675291775782e-15L,
320 -6.1602987721551372146677749697118997573853e-01L,
321 -2.9818078843122551067455400545109858745295e-16L,
322 -6.2024040975185457114093878772109746932983e-01L,
323 -2.9577105558448461493874424529516311623184e-15L,
324 -6.2443328801189323939979658462107181549072e-01L,
325 -2.6164274215943360130441858075903119505815e-16L,
326 -6.2860865942237253989333112258464097976685e-01L,
327 -1.5978509770831895426601797458058854400463e-15L,
328 -6.3276666957103699928666173946112394332886e-01L,
329 -8.3025912472904245581515990140161946934461e-16L,
330 -6.3690746223706895534633076749742031097412e-01L,
331 -2.7627416365968377888021629180796328536455e-16L,
332 -6.4103117942092779912854894064366817474365e-01L,
333 -3.4919270523937617243719652995048419893186e-15L,
334 -6.4513796137358170312836591619998216629028e-01L,
335 -2.9985368625799347497396478978681548584217e-15L,
336 -6.4922794662510696639401430729776620864868e-01L,
337 -2.8524968256626075449136225882322854909611e-15L,
338 -6.5330127201274379444839723873883485794067e-01L,
339 -1.8443102186424720390266302263929355424008e-15L,
340 -6.5735807270835877602621621917933225631714e-01L,
341 -1.2541156738040666039091970075936624723645e-15L,
342 -6.6139848224536379461824253667145967483521e-01L,
343 -1.2136419933020381912633127333149145382797e-15L,
344 -6.6542263254508782210905337706208229064941e-01L,
345 -2.6268410392329445778904988886114643307320e-15L,
346 -6.6943065394262646350398426875472068786621e-01L,
347 -2.8037949010021747828222575923191438798877e-15L,
348 -6.7342267521216570003161905333399772644043e-01L,
349 -1.0202663413354670195383104149875619397268e-15L,
350 -6.7739882359180469961756898555904626846313e-01L,
351 -1.4411921136244383020300914304078010801275e-15L,
352 -6.8135922480790256372529256623238325119019e-01L,
353 -5.0522277899333570619054540068138110661023e-16L,
354 -6.8530400309891703614084690343588590621948e-01L,
355 -2.3804032011755313470802014258958896193599e-15L,
356 -6.8923328123880622797514661215245723724365e-01L,
357 -2.7523497677256621466659891416404053623832e-15L,
102 + 0.0000000000000000000000000000000000000000e+00L,
103 + 0.0000000000000000000000000000000000000000e+00L,
104 + 7.7821404420532758194894995540380477905273e-03L,
105 + 1.6731279734005070987158875984584325351222e-15L,
106 + 1.5504186535963526694104075431823730468750e-02L,
107 + 1.7274567499706106231054091184928671990316e-15L,
108 + 2.3167059281533397552266251295804977416992e-02L,
109 + 9.8067653290966648493916241687661877474892e-16L,
110 + 3.0771658666751022792595904320478439331055e-02L,
111 + 2.6655784323032762937247606420524589813624e-15L,
112 + 3.8318864302134159061097307130694389343262e-02L,
113 + 2.4401326580179931029010027013316092332340e-15L,
114 + 4.5809536031292452662455616518855094909668e-02L,
115 + 1.7505042236510958082472042641283104263139e-15L,
116 + 5.3244514518809182845870964229106903076172e-02L,
117 + 3.1000199992295574218738634002122149891138e-15L,
118 + 6.0624621816433688081815489567816257476807e-02L,
119 + 1.1544987906424726040058093958345197512800e-15L,
120 + 6.7950661908504628172522643581032752990723e-02L,
121 + 3.1212220426341915966610439115772728417386e-15L,
122 + 7.5223421237584631171557703055441379547119e-02L,
123 + 2.8945270476369282210350897509258766743153e-15L,
124 + 8.2443669211073711267090402543544769287109e-02L,
125 + 8.8000106966612476303662698634483335676886e-16L,
126 + 8.9612158689686083334891009144484996795654e-02L,
127 + 1.0492850604602339995319895311151740799226e-15L,
128 + 9.6729626458550654888313147239387035369873e-02L,
129 + 4.5740725790924807640164516707244620870662e-16L,
130 + 1.0379679368164218544734467286616563796997e-01L,
131 + 1.3793787171308978090503366050174239822054e-15L,
132 + 1.1081436634028918319927470292896032333374e-01L,
133 + 9.3099553146639425160476473362380086036919e-16L,
134 + 1.1778303565638026384476688690483570098877e-01L,
135 + 3.1906940272225656860040797111813146690890e-15L,
136 + 1.2470347850095464536934741772711277008057e-01L,
137 + 2.5904940590976537504984110469214193890052e-15L,
138 + 1.3157635778871679121948545798659324645996e-01L,
139 + 2.4813692306707028899159917911012100567219e-15L,
140 + 1.3840232285911824305912887211889028549194e-01L,
141 + 8.9262619700148275890190121571708972000380e-16L,
142 + 1.4518200984449691759436973370611667633057e-01L,
143 + 9.7968756533003444764719201050911636480025e-16L,
144 + 1.5191604202583874894116888754069805145264e-01L,
145 + 3.2261306345373561864598749471119213018106e-15L,
146 + 1.5860503017663774016909883357584476470947e-01L,
147 + 8.4392427234104999681053621980394827998735e-16L,
148 + 1.6524957289530561865831259638071060180664e-01L,
149 + 1.5442172988528965297119225948270579746101e-15L,
150 + 1.7185025692665689689420105423778295516968e-01L,
151 + 2.3254458978918173643097657009894831132739e-15L,
152 + 1.7840765747281750464026117697358131408691e-01L,
153 + 7.9247913906453736065426776912520942036896e-16L,
154 + 1.8492233849401173984006163664162158966064e-01L,
155 + 2.5282384195601762803134514624610774126020e-16L,
156 + 1.9139485299962899489401024766266345977783e-01L,
157 + 4.5971528855989864541366920731297729269228e-16L,
158 + 1.9782574332991842425144568551331758499146e-01L,
159 + 1.4561111263856836438840838027526567191527e-15L,
160 + 2.0421554142868814096800633706152439117432e-01L,
161 + 2.7505358140491347148810394262840919337709e-15L,
162 + 2.1056476910734645002776233013719320297241e-01L,
163 + 3.1876417904825951583107481283088861928977e-15L,
164 + 2.1687393830061196808856038842350244522095e-01L,
165 + 2.3915305291373208450532580201045871599499e-15L,
166 + 2.2314355131420882116799475625157356262207e-01L,
167 + 9.3459830033405826094075253077304795996257e-16L,
168 + 2.2937410106484534821902343537658452987671e-01L,
169 + 4.8177245728966955534167425511952551974164e-16L,
170 + 2.3556607131276408040321257431060075759888e-01L,
171 + 2.8286743756446304426525380844720043381780e-15L,
172 + 2.4171993688714366044223424978554248809814e-01L,
173 + 1.5077020732661279714120052415509585052975e-15L,
174 + 2.4783616390458007572306087240576744079590e-01L,
175 + 1.1810575418933407573072030113600980623171e-15L,
176 + 2.5391520998096339667426946107298135757446e-01L,
177 + 4.7463053836833625309891834934881898560705e-17L,
178 + 2.5995752443692410338371701072901487350464e-01L,
179 + 1.9635883624838132961710716735786266795913e-15L,
180 + 2.6596354849713677026556979399174451828003e-01L,
181 + 1.1710735561325457988709887923652142233351e-15L,
182 + 2.7193371548364098089223261922597885131836e-01L,
183 + 7.7793943687530702031066421537496360004376e-16L,
184 + 2.7786845100345303194444568362087011337280e-01L,
185 + 3.2742419043493025311197092322146237692165e-15L,
186 + 2.8376817313064250924981024581938982009888e-01L,
187 + 2.0890970909765308649465619266075677112425e-15L,
188 + 2.8963329258304071345264674164354801177979e-01L,
189 + 1.9634262463138821209582240742801727823629e-15L,
190 + 2.9546421289383317798638017848134040832520e-01L,
191 + 2.6984003017275736237868564402005801750600e-15L,
192 + 3.0126133057816062432721082586795091629028e-01L,
193 + 1.1566856647123658045763670687640673680383e-15L,
194 + 3.0702503529490954292668902780860662460327e-01L,
195 + 2.3191484355127267712770857311812090801833e-15L,
196 + 3.1275571000389490450288576539605855941772e-01L,
197 + 1.9838833607942922604727420618882220398852e-15L,
198 + 3.1845373111853447767316538374871015548706e-01L,
199 + 1.3813708182984188944010814590398164268227e-16L,
200 + 3.2411946865421015218089451082050800323486e-01L,
201 + 1.8239097762496144793489474731253815376404e-15L,
202 + 3.2975328637246548169059678912162780761719e-01L,
203 + 2.5001238260227991620033344720809714552230e-15L,
204 + 3.3535554192113536942088103387504816055298e-01L,
205 + 2.4608362985459391180385214539620341910962e-15L,
206 + 3.4092658697059263772644044365733861923218e-01L,
207 + 5.7257864875612301758921090406373771458003e-16L,
208 + 3.4646676734620740489845047704875469207764e-01L,
209 + 1.1760200117113770182586341947822306069951e-15L,
210 + 3.5197642315717558858523261733353137969971e-01L,
211 + 2.5960702148389259075462896448369304790506e-15L,
212 + 3.5745588892180180096147523727267980575562e-01L,
213 + 1.9732645342528682246686790561260072184839e-15L,
214 + 3.6290549368936808605212718248367309570312e-01L,
215 + 3.6708569716349381675043725477739939978160e-16L,
216 + 3.6832556115870573876236448995769023895264e-01L,
217 + 1.9142858656640927085879445412821643247628e-15L,
218 + 3.7371640979358389245135185774415731430054e-01L,
219 + 1.8836966497497166619234389157276681281343e-16L,
220 + 3.7907835293496816575498087331652641296387e-01L,
221 + 1.2926358724723144934459175417385013725801e-15L,
222 + 3.8441169891033055705520382616668939590454e-01L,
223 + 1.4826795862363146014726140088145939341729e-15L,
224 + 3.8971675114002479745067830663174390792847e-01L,
225 + 4.1591978529737177695912258866565331189698e-16L,
226 + 3.9499380824086571806219581048935651779175e-01L,
227 + 3.2600441982258756252505182317625310732365e-15L,
228 + 4.0024316412701210765590076334774494171143e-01L,
229 + 5.9927342433864738622836851475469574662703e-16L,
230 + 4.0546510810816371872533636633306741714478e-01L,
231 + 6.6325267674913128171942721503283748008372e-16L,
232 + 4.1065992498526782128465129062533378601074e-01L,
233 + 5.6464965491255048900165082436455718077885e-16L,
234 + 4.1582789514371043537721561733633279800415e-01L,
235 + 5.3023611327561856950735176370587227509442e-16L,
236 + 4.2096929464412724541944044176489114761353e-01L,
237 + 2.3907094267197419048248363335257046791153e-15L,
238 + 4.2608439531089814522601955104619264602661e-01L,
239 + 1.9178985253285492839728700574592375309985e-15L,
240 + 4.3117346481836804628073878120630979537964e-01L,
241 + 3.2945784336977492852031005044499611665595e-15L,
242 + 4.3623676677491474151793227065354585647583e-01L,
243 + 3.3288311090524075754441878570852962903891e-15L,
244 + 4.4127456080487448275562201160937547683716e-01L,
245 + 7.4673387443005192574852544613692268411229e-16L,
246 + 4.4628710262841764233598951250314712524414e-01L,
247 + 1.8691966006681165218815050615460959199251e-15L,
248 + 4.5127464413945617138779198285192251205444e-01L,
249 + 2.4137569004002270899666314791611479063976e-15L,
250 + 4.5623743348158640742440184112638235092163e-01L,
251 + 1.1869564036970375473975162509216610120281e-15L,
252 + 4.6117571512216670726047595962882041931152e-01L,
253 + 3.4591075239659690349392915732654828400811e-15L,
254 + 4.6608972992459740680715185590088367462158e-01L,
255 + 1.8177514673916038857252366108673570603067e-15L,
256 + 4.7097971521878889689105562865734100341797e-01L,
257 + 2.1156558422273990182479555421331461933366e-15L,
258 + 4.7584590486996347635795245878398418426514e-01L,
259 + 4.3790725712752039722791012358345927696967e-16L,
260 + 4.8068852934575190261057286988943815231323e-01L,
261 + 5.0660455855585733988956280680891477171499e-18L,
262 + 4.8550781578169832641833636444061994552612e-01L,
263 + 2.4813834547127501689550526444948043590905e-15L,
264 + 4.9030398804519137456736643798649311065674e-01L,
265 + 2.4635829797216592537498738468934647345741e-15L,
266 + 4.9507726679784980206022737547755241394043e-01L,
267 + 1.7125377372093652812514167461480115600063e-15L,
268 + 4.9982786955644797899367404170334339141846e-01L,
269 + 1.3508276573735437007500942002018098437396e-15L,
270 + 5.0455601075239187025545106735080480575562e-01L,
271 + 3.4168028574643873701242268618467347998876e-15L,
272 + 5.0926190178980590417268103919923305511475e-01L,
273 + 2.0426313938800290907697638200502614622891e-15L,
274 + 5.1394575110223428282552049495279788970947e-01L,
275 + 3.3975485593321419703400672813719873194659e-17L,
276 + 5.1860776420804555186805373523384332656860e-01L,
277 + 8.0284923261130955371987633083003284697416e-17L,
278 + 5.2324814376454753528378205373883247375488e-01L,
279 + 3.0123302517119603836788558832352723470118e-16L,
280 + 5.2786708962084105678513878956437110900879e-01L,
281 + 1.3283287534282139298545497336570406582397e-15L,
282 + 5.3246479886946929127589100971817970275879e-01L,
283 + 2.5525980327137419625398485590148417041921e-15L,
284 + 5.3704146589688050994482182431966066360474e-01L,
285 + 3.1446219074198341716354190061340477078626e-15L,
286 + 5.4159728243274329884116014000028371810913e-01L,
287 + 1.0727353821639001503808606766770295812627e-15L,
288 + 5.4613243759813556721383065450936555862427e-01L,
289 + 8.3168566554721843605240702438699163825794e-17L,
290 + 5.5064711795266063631970610003918409347534e-01L,
291 + 1.6429402420791657293666192255419538448840e-15L,
292 + 5.5514150754050106684189813677221536636353e-01L,
293 + 5.2587358222274368868380660194332415847228e-16L,
294 + 5.5961578793542088305912329815328121185303e-01L,
295 + 1.8032117652023735453816330571171114110385e-15L,
296 + 5.6407013828480145889443519990891218185425e-01L,
297 + 1.5071769490901812785299634348367857600711e-15L,
298 + 5.6850473535266843327917740680277347564697e-01L,
299 + 2.7879956135806418878792935692629147550413e-16L,
300 + 5.7291975356178426181941176764667034149170e-01L,
301 + 1.2472733449589795907271346997596471822345e-15L,
302 + 5.7731536503482061561953742057085037231445e-01L,
303 + 2.9886985746409486460291929160223207644146e-15L,
304 + 5.8169173963462128540413687005639076232910e-01L,
305 + 1.1971164738836689815783808674399742176950e-15L,
306 + 5.8604904500357690722012193873524665832520e-01L,
307 + 1.3016839974975520776911897855504474452726e-15L,
308 + 5.9038744660217545856539800297468900680542e-01L,
309 + 9.1607651870514890975077236127894522134392e-16L,
310 + 5.9470710774668944509357970673590898513794e-01L,
311 + 3.3444207638397932963480545729233567201211e-15L,
312 + 5.9900818964608149030937056522816419601440e-01L,
313 + 1.9090722294592334873060460706130642200729e-15L,
314 + 6.0329085143808214297678205184638500213623e-01L,
315 + 2.1193638031348149256035110177854940281795e-15L,
316 + 6.0755525022453937822319858241826295852661e-01L,
317 + 2.4172778865703728624133665395876418941354e-15L,
318 + 6.1180154110599005434778518974781036376953e-01L,
319 + 2.8491821045766810044199163148675291775782e-15L,
320 + 6.1602987721551372146677749697118997573853e-01L,
321 + 2.9818078843122551067455400545109858745295e-16L,
322 + 6.2024040975185457114093878772109746932983e-01L,
323 + 2.9577105558448461493874424529516311623184e-15L,
324 + 6.2443328801189323939979658462107181549072e-01L,
325 + 2.6164274215943360130441858075903119505815e-16L,
326 + 6.2860865942237253989333112258464097976685e-01L,
327 + 1.5978509770831895426601797458058854400463e-15L,
328 + 6.3276666957103699928666173946112394332886e-01L,
329 + 8.3025912472904245581515990140161946934461e-16L,
330 + 6.3690746223706895534633076749742031097412e-01L,
331 + 2.7627416365968377888021629180796328536455e-16L,
332 + 6.4103117942092779912854894064366817474365e-01L,
333 + 3.4919270523937617243719652995048419893186e-15L,
334 + 6.4513796137358170312836591619998216629028e-01L,
335 + 2.9985368625799347497396478978681548584217e-15L,
336 + 6.4922794662510696639401430729776620864868e-01L,
337 + 2.8524968256626075449136225882322854909611e-15L,
338 + 6.5330127201274379444839723873883485794067e-01L,
339 + 1.8443102186424720390266302263929355424008e-15L,
340 + 6.5735807270835877602621621917933225631714e-01L,
341 + 1.2541156738040666039091970075936624723645e-15L,
342 + 6.6139848224536379461824253667145967483521e-01L,
343 + 1.2136419933020381912633127333149145382797e-15L,
344 + 6.6542263254508782210905337706208229064941e-01L,
345 + 2.6268410392329445778904988886114643307320e-15L,
346 + 6.6943065394262646350398426875472068786621e-01L,
347 + 2.8037949010021747828222575923191438798877e-15L,
348 + 6.7342267521216570003161905333399772644043e-01L,
349 + 1.0202663413354670195383104149875619397268e-15L,
350 + 6.7739882359180469961756898555904626846313e-01L,
351 + 1.4411921136244383020300914304078010801275e-15L,
352 + 6.8135922480790256372529256623238325119019e-01L,
353 + 5.0522277899333570619054540068138110661023e-16L,
354 + 6.8530400309891703614084690343588590621948e-01L,
355 + 2.3804032011755313470802014258958896193599e-15L,
356 + 6.8923328123880622797514661215245723724365e-01L,
357 + 2.7523497677256621466659891416404053623832e-15L,
358 358 };
359 359
360 360 /*
361 361 * Compute N*log2 + log(1+zk+zh+zt) in extra precision
362 362 */
363 -static long double k_log_NKzl(int N, int K, long double zh, long double *zt)
363 +static long double
364 +k_log_NKzl(int N, int K, long double zh, long double *zt)
364 365 {
365 366 long double y, r, w, s2, s2h, s2t, t, zk, v, P;
366 367 double dzk;
367 368
368 369 #if !defined(__x86)
369 370 unsigned lx, ly;
370 371 int j;
371 372 #endif
372 373
373 374 ((int *)&dzk)[HIWORD] = 0x3ff00000 + (K << 13);
374 375 ((int *)&dzk)[LOWORD] = 0;
375 - t = zh + (*zt);
376 - zk = (long double) dzk;
376 + t = zh + (*zt);
377 + zk = (long double)dzk;
377 378 r = two / (t + two * zk);
378 379 s2h = s2 = r * t;
379 380 /* split s2 into correctly rounded half */
380 381
381 382 #if defined(__x86)
382 - ((unsigned *)&s2h)[0] = 0; /* 32 bits chopped */
383 + ((unsigned *)&s2h)[0] = 0; /* 32 bits chopped */
383 384 #else
384 -
385 - lx = ((unsigned *)&s2h)[2]; /* 56 bits rounded */
386 - j = ((lx >> 24) + 1) >> 1;
387 - ((unsigned *)&s2h)[2] = (j << 25);
388 - lx = ((unsigned *)&s2h)[1];
389 - ly = lx + (j >> 7);
390 - ((unsigned *)&s2h)[1] = ly;
391 - ((unsigned *)&s2h)[0] += (ly == 0 && lx != 0);
392 - ((unsigned *)&s2h)[3] = 0;
385 + lx = ((unsigned *)&s2h)[2]; /* 56 bits rounded */
386 + j = ((lx >> 24) + 1) >> 1;
387 + ((unsigned *)&s2h)[2] = (j << 25);
388 + lx = ((unsigned *)&s2h)[1];
389 + ly = lx + (j >> 7);
390 + ((unsigned *)&s2h)[1] = ly;
391 + ((unsigned *)&s2h)[0] += (ly == 0 && lx != 0);
392 + ((unsigned *)&s2h)[3] = 0;
393 393 #endif
394 394
395 395 v = half * s2h;
396 396 w = s2 * s2;
397 397 s2t = r * ((((zh - s2h * zk) - v * zh) + (*zt)) - v * (*zt));
398 - P = s2t + (w * s2) * ((P1 + w * P2) + (w * w) * ((P3 + w * P4)
399 - + (w * w) * (P5 + w * P6 + (w * w) * P7)));
398 + P = s2t + (w * s2) * ((P1 + w * P2) + (w * w) * ((P3 + w * P4) + (w *
399 + w) * (P5 + w * P6 + (w * w) * P7)));
400 400 P += N * ln2_t + TBL_log1k[K + K + 1];
401 - t = N*ln2_h + TBL_log1k[K+K];
401 + t = N * ln2_h + TBL_log1k[K + K];
402 402 y = t + (P + s2h);
403 403 P -= ((y - t) - s2h);
404 404 *zt = P;
405 405 return (y);
406 406 }
407 407
408 408 long double
409 409 __k_clog_rl(long double x, long double y, long double *er)
410 410 {
411 411 long double t1, t2, t3, t4, tk, z, wh, w, zh, zk;
412 412 int n, k, ix, iy, iz, nx, ny, nz, i;
413 413 double dk;
414 414
415 415 #if !defined(__x86)
416 416 int j;
417 417 unsigned lx, ly;
418 418 #endif
419 419
420 420 ix = HI_XWORD(x) & ~0x80000000;
421 421 iy = HI_XWORD(y) & ~0x80000000;
422 - y = fabsl(y); x = fabsl(x);
422 + y = fabsl(y);
423 + x = fabsl(x);
424 +
423 425 if (ix < iy || (ix < 0x7fff0000 && ix == iy && x < y)) {
424 426 /* force x >= y */
425 - tk = x; x = y; y = tk;
426 - n = ix, ix = iy; iy = n;
427 + tk = x;
428 + x = y;
429 + y = tk;
430 + n = ix, ix = iy;
431 + iy = n;
427 432 }
433 +
428 434 *er = zero;
429 - nx = ix >> 16; ny = iy >> 16;
430 - if (nx >= 0x7fff) { /* x or y is Inf or NaN */
435 + nx = ix >> 16;
436 + ny = iy >> 16;
437 +
438 + if (nx >= 0x7fff) { /* x or y is Inf or NaN */
431 439 if (isinfl(x))
432 440 return (x);
433 441 else if (isinfl(y))
434 442 return (y);
435 443 else
436 - return (x+y);
444 + return (x + y);
437 445 }
446 +
438 447 /*
439 448 * for tiny y:(double y < 2^-35, extended y < 2^-46, quad y < 2^-70)
440 449 *
441 450 * log(sqrt(1 + y**2)) = y**2 / 2 - y**4 / 8 + ... = y**2 / 2
442 451 */
443 452 #if defined(__x86)
444 453 if (x == 1.0L && ny < (0x3fff - 46)) {
445 454 #else
446 455 if (x == 1.0L && ny < (0x3fff - 70)) {
447 456 #endif
448 457
449 458 t2 = y * y;
459 +
450 460 if (ny >= 8305) { /* compute er = tail of t2 */
451 - dk = (double) y;
461 + dk = (double)y;
452 462
453 463 #if defined(__x86)
454 464 ((unsigned *)&dk)[LOWORD] &= 0xfffe0000;
455 465 #endif
456 466
457 - wh = (long double) dk;
467 + wh = (long double)dk;
458 468 *er = half * ((y - wh) * (y + wh) - (t2 - wh * wh));
459 469 }
470 +
460 471 return (half * t2);
461 472 }
473 +
462 474 /*
463 475 * x or y is subnormal or zero
464 476 */
465 477 if (nx == 0) {
466 - if (x == 0.0L)
478 + if (x == 0.0L) {
467 479 return (-1.0L / x);
468 - else {
480 + } else {
469 481 x *= two240;
470 482 y *= two240;
471 483 ix = HI_XWORD(x);
472 484 iy = HI_XWORD(y);
473 485 nx = (ix >> 16) - 240;
474 486 ny = (iy >> 16) - 240;
487 +
475 488 /* guard subnormal flush to 0 */
476 489 if (x == 0.0L)
477 490 return (-1.0L / x);
478 491 }
479 - } else if (ny == 0) { /* y subnormal, scale it */
492 + } else if (ny == 0) { /* y subnormal, scale it */
480 493 y *= two240;
481 494 iy = HI_XWORD(y);
482 495 ny = (iy >> 16) - 240;
483 496 }
497 +
484 498 n = nx - ny;
499 +
485 500 /*
486 501 * When y is zero or when x >> y, i.e., n > 62, 78, 122 for DBLE,
487 502 * EXTENDED, QUAD respectively,
488 503 * log(x) = log(sqrt(x * x + y * y)) to 27 extra bits.
489 504 */
490 505
491 506 #if defined(__x86)
492 507 if (n > 78 || y == 0.0L) {
493 508 #else
494 509 if (n > 122 || y == 0.0L) {
495 510 #endif
496 511
497 512 XFSCALE(x, (0x3fff - (ix >> 16)));
498 - i = ((ix & 0xffff) + 0x100) >> 9; /* 7.5 bits of x */
499 - zk = 1.0L + ((long double) i) * 0.0078125L;
513 + i = ((ix & 0xffff) + 0x100) >> 9; /* 7.5 bits of x */
514 + zk = 1.0L + ((long double)i) * 0.0078125L;
500 515 z = x - zk;
501 516 dk = (double)z;
502 517
503 518 #if defined(__x86)
504 519 ((unsigned *)&dk)[LOWORD] &= 0xfffe0000;
505 520 #endif
506 521
507 522 zh = (long double)dk;
508 - k = i & 0x7f; /* index of zk */
523 + k = i & 0x7f; /* index of zk */
509 524 n = nx - 0x3fff;
510 525 *er = z - zh;
526 +
511 527 if (i == 0x80) { /* if zk = 2.0, adjust scaling */
512 528 n += 1;
513 - zh *= 0.5L; *er *= 0.5L;
529 + zh *= 0.5L;
530 + *er *= 0.5L;
514 531 }
532 +
515 533 w = k_log_NKzl(n, k, zh, er);
516 534 } else {
517 535 /*
518 536 * compute z = x*x + y*y
519 537 */
520 538 XFSCALE(x, (0x3fff - (ix >> 16)));
521 539 XFSCALE(y, (0x3fff - n - (iy >> 16)));
522 540 ix = (ix & 0xffff) | 0x3fff0000;
523 541 iy = (iy & 0xffff) | (0x3fff0000 - (n << 16));
524 542 nx -= 0x3fff;
525 - t1 = x * x; t2 = y * y;
543 + t1 = x * x;
544 + t2 = y * y;
526 545 wh = x;
527 546
528 547 /* split x into correctly rounded half */
529 548 #if defined(__x86)
530 549 ((unsigned *)&wh)[0] = 0; /* 32 bits chopped */
531 550 #else
532 551 lx = ((unsigned *)&wh)[2]; /* 56 rounded */
533 - j = ((lx >> 24) + 1) >> 1;
552 + j = ((lx >> 24) + 1) >> 1;
534 553 ((unsigned *)&wh)[2] = (j << 25);
535 554 lx = ((unsigned *)&wh)[1];
536 555 ly = lx + (j >> 7);
537 556 ((unsigned *)&wh)[1] = ly;
538 557 ((unsigned *)&wh)[0] += (ly == 0 && lx != 0);
539 558 ((unsigned *)&wh)[3] = 0;
540 559 #endif
541 560
542 - z = t1+t2;
561 + z = t1 + t2;
562 +
543 563 /*
544 564 * higher precision simulation x*x = t1 + t3, y*y = t2 + t4
545 565 */
546 566 tk = wh - x;
547 567 t3 = tk * tk - (two * wh * tk - (wh * wh - t1));
548 568 wh = y;
549 569
550 570 /* split y into correctly rounded half */
551 571 #if defined(__x86)
552 572 ((unsigned *)&wh)[0] = 0; /* 32 bits chopped */
553 573 #else
554 574 ly = ((unsigned *)&wh)[2]; /* 56 bits rounded */
555 - j = ((ly >> 24) + 1) >> 1;
575 + j = ((ly >> 24) + 1) >> 1;
556 576 ((unsigned *)&wh)[2] = (j << 25);
557 577 lx = ((unsigned *)&wh)[1];
558 578 ly = lx + (j >> 7);
559 579 ((unsigned *)&wh)[1] = ly;
560 580 ((unsigned *)&wh)[0] += (ly == 0 && lx != 0);
561 581 ((unsigned *)&wh)[3] = 0;
562 582 #endif
563 583
564 584 tk = wh - y;
565 585 t4 = tk * tk - (two * wh * tk - (wh * wh - t2));
586 +
566 587 /*
567 588 * find zk matches z to 7.5 bits
568 589 */
569 590 iz = HI_XWORD(z);
570 - k = ((iz & 0xffff) + 0x100) >> 9; /* 7.5 bits of x */
591 + k = ((iz & 0xffff) + 0x100) >> 9; /* 7.5 bits of x */
571 592 nz = (iz >> 16) - 0x3fff + (k >> 7);
572 593 k &= 0x7f;
573 - zk = 1.0L + ((long double) k) * 0.0078125L;
574 - if (nz == 1) zk += zk;
575 - else if (nz == 2) zk *= 4.0L;
576 - else if (nz == 3) zk *= 8.0L;
594 + zk = 1.0L + ((long double)k) * 0.0078125L;
595 +
596 + if (nz == 1)
597 + zk += zk;
598 + else if (nz == 2)
599 + zk *= 4.0L;
600 + else if (nz == 3)
601 + zk *= 8.0L;
602 +
577 603 /*
578 604 * order t1, t2, t3, t4 according to their size
579 605 */
580 606 if (t2 >= fabsl(t3)) {
581 607 if (fabsl(t3) < fabsl(t4)) {
582 - wh = t3; t3 = t4; t4 = wh;
608 + wh = t3;
609 + t3 = t4;
610 + t4 = wh;
583 611 }
584 612 } else {
585 - wh = t2; t2 = t3; t3 = wh;
613 + wh = t2;
614 + t2 = t3;
615 + t3 = wh;
586 616 }
617 +
587 618 /*
588 619 * higher precision simulation: x * x + y * y = t1 + t2 + t3 + t4
589 620 * = zk(7 bits) + zh(24 bits) + *er(tail) and call k_log_NKz
590 621 */
591 622 tk = t1 - zk;
592 623 zh = ((tk + t2) + t3) + t4;
593 624
594 625 /* split zh into correctly rounded half */
595 626 #if defined(__x86)
596 627 ((unsigned *)&zh)[0] = 0;
597 628 #else
598 629 ly = ((unsigned *)&zh)[2];
599 - j = ((ly >> 24) + 1) >> 1;
630 + j = ((ly >> 24) + 1) >> 1;
600 631 ((unsigned *)&zh)[2] = (j << 25);
601 632 lx = ((unsigned *)&zh)[1];
602 633 ly = lx + (j >> 7);
603 634 ((unsigned *)&zh)[1] = ly;
604 635 ((unsigned *)&zh)[0] += (ly == 0 && lx != 0);
605 636 ((unsigned *)&zh)[3] = 0;
606 637 #endif
607 638
608 639 w = fabsl(zh);
609 - if (w >= fabsl(t2))
610 -{
611 - *er = (((tk - zh) + t2) + t3) + t4;
612 -}
613 -
614 - else {
615 640
641 + if (w >= fabsl(t2)) {
642 + *er = (((tk - zh) + t2) + t3) + t4;
643 + } else {
616 644 if (n == 0) {
617 645 wh = half * zk;
618 646 wh = (t1 - wh) - (wh - t2);
619 - } else
647 + } else {
620 648 wh = tk + t2;
621 - if (w >= fabsl(t3))
649 + }
650 +
651 + if (w >= fabsl(t3)) {
622 652 *er = ((wh - zh) + t3) + t4;
623 - else {
653 + } else {
624 654 z = t3;
625 655 t3 += t4;
626 656 t4 -= t3 - z;
657 +
627 658 if (w >= fabsl(t3))
628 659 *er = ((wh - zh) + t3) + t4;
629 660 else
630 661 *er = ((wh + t3) - zh) + t4;
631 662 }
632 663 }
664 +
633 665 if (nz == 3) {
634 - zh *= 0.125L; *er *= 0.125L;
666 + zh *= 0.125L;
667 + *er *= 0.125L;
635 668 } else if (nz == 2) {
636 - zh *= 0.25L; *er *= 0.25L;
669 + zh *= 0.25L;
670 + *er *= 0.25L;
637 671 } else if (nz == 1) {
638 - zh *= half; *er *= half;
672 + zh *= half;
673 + *er *= half;
639 674 }
675 +
640 676 nz += nx + nx;
641 677 w = half * k_log_NKzl(nz, k, zh, er);
642 678 *er *= half;
643 679 }
680 +
644 681 return (w);
645 682 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX