| blob | e353b08658d22cb9be5fe0af1be7bb8883e22caf |
1 .file "powf.s"
4 // Copyright (c) 2000 - 2005, Intel Corporation
5 // All rights reserved.
6 //
7 //
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
10 // met:
11 //
12 // * Redistributions of source code must retain the above copyright
13 // notice, this list of conditions and the following disclaimer.
14 //
15 // * Redistributions in binary form must reproduce the above copyright
16 // notice, this list of conditions and the following disclaimer in the
17 // documentation and/or other materials provided with the distribution.
18 //
19 // * The name of Intel Corporation may not be used to endorse or promote
20 // products derived from this software without specific prior written
21 // permission.
23 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
26 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
27 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
28 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
29 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
30 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
31 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
32 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34 //
35 // Intel Corporation is the author of this code, and requests that all
36 // problem reports or change requests be submitted to it directly at
37 // http://www.intel.com/software/products/opensource/libraries/num.htm.
38 //
39 // History
40 //==============================================================
41 // 02/02/00 Initial version
42 // 02/03/00 Added p12 to definite over/under path. With odd power we did not
43 // maintain the sign of x in this path.
44 // 04/04/00 Unwind support added
45 // 04/19/00 pow(+-1,inf) now returns NaN
46 // pow(+-val, +-inf) returns 0 or inf, but now does not call error
47 // support
48 // Added s1 to fcvt.fx because invalid flag was incorrectly set.
49 // 08/15/00 Bundle added after call to __libm_error_support to properly
50 // set [the previously overwritten] GR_Parameter_RESULT.
51 // 09/07/00 Improved performance by eliminating bank conflicts and other stalls,
52 // and tweaking the critical path
53 // 09/08/00 Per c99, pow(+-1,inf) now returns 1, and pow(+1,nan) returns 1
54 // 09/28/00 Updated NaN**0 path
55 // 01/20/01 Fixed denormal flag settings.
56 // 02/13/01 Improved speed.
57 // 03/19/01 Reordered exp polynomial to improve speed and eliminate monotonicity
58 // problem in round up, down, and to zero modes. Also corrected
59 // overflow result when x negative, y odd in round up, down, zero.
60 // 06/14/01 Added brace missing from bundle
61 // 12/10/01 Corrected case where x negative, 2^23 <= |y| < 2^24, y odd integer.
62 // 02/08/02 Fixed overflow/underflow cases that were not calling error support.
63 // 05/20/02 Cleaned up namespace and sf0 syntax
64 // 08/29/02 Improved Itanium 2 performance
65 // 02/10/03 Reordered header: .section, .global, .proc, .align
66 // 10/09/03 Modified algorithm to improve performance, reduce table size, and
67 // fix boundary case powf(2.0,-150.0)
68 // 03/31/05 Reformatted delimiters between data tables
69 //
70 // API
71 //==============================================================
72 // float powf(float x, float y)
73 //
74 // Overview of operation
75 //==============================================================
76 //
77 // Three steps...
78 // 1. Log(x)
79 // 2. y Log(x)
80 // 3. exp(y log(x))
81 //
82 // This means we work with the absolute value of x and merge in the sign later.
83 // Log(x) = G + delta + r -rsq/2 + p
84 // G,delta depend on the exponent of x and table entries. The table entries are
85 // indexed by the exponent of x, called K.
86 //
87 // The G and delta come out of the reduction; r is the reduced x.
88 //
89 // B = frcpa(x)
90 // xB-1 is small means that B is the approximate inverse of x.
91 //
92 // Log(x) = Log( (1/B)(Bx) )
93 // = Log(1/B) + Log(Bx)
94 // = Log(1/B) + Log( 1 + (Bx-1))
95 //
96 // x = 2^K 1.x_1x_2.....x_52
97 // B= frcpa(x) = 2^-k Cm
98 // Log(1/B) = Log(1/(2^-K Cm))
99 // Log(1/B) = Log((2^K/ Cm))
100 // Log(1/B) = K Log(2) + Log(1/Cm)
101 //
102 // Log(x) = K Log(2) + Log(1/Cm) + Log( 1 + (Bx-1))
103 //
104 // If you take the significand of x, set the exponent to true 0, then Cm is
105 // the frcpa. We tabulate the Log(1/Cm) values. There are 256 of them.
106 // The frcpa table is indexed by 8 bits, the x_1 thru x_8.
107 // m = x_1x_2...x_8 is an 8-bit index.
108 //
109 // Log(1/Cm) = log(1/frcpa(1+m/256)) where m goes from 0 to 255.
110 //
111 // We tabulate as one double, T for single precision power
112 //
113 // Log(x) = (K Log(2)_hi + T) + (K Log(2)_lo) + Log( 1 + (Bx-1))
114 // Log(x) = G + delta + Log( 1 + (Bx-1))
115 //
116 // The Log( 1 + (Bx-1)) can be calculated as a series in r = Bx-1.
117 //
118 // Log( 1 + (Bx-1)) = r - rsq/2 + p
119 // where p = r^3(P0 + P1*r + P2*r^2)
120 //
121 // Then,
122 //
123 // yLog(x) = yG + y delta + y(r-rsq/2) + yp
124 // yLog(x) = Z1 + e3 + Z2 + Z3
125 //
126 //
127 // exp(yLog(x)) = exp(Z1 + Z2) exp(Z3) exp(e3)
128 //
129 //
130 // exp(Z3) is another series.
131 // exp(e3) is approximated as f3 = 1 + e3
132 //
133 // exp(Z1 + Z2) = exp(Z)
134 // Z (128/log2) = number of log2/128 in Z is N
135 //
136 // s = Z - N log2/128
137 //
138 // exp(Z) = exp(s) exp(N log2/128)
139 //
140 // exp(r) = exp(Z - N log2/128)
141 //
142 // r = s + d = (Z - N (log2/128)_hi) -N (log2/128)_lo
143 // = Z - N (log2/128)
144 //
145 // Z = s+d +N (log2/128)
146 //
147 // exp(Z) = exp(s) (1+d) exp(N log2/128)
148 //
149 // N = M 128 + n
150 //
151 // N log2/128 = M log2 + n log2/128
152 //
153 // n is 8 binary digits = n_7n_6...n_1
154 //
155 // n log2/128 = n_7n_6n_5 16 log2/128 + n_4n_3n_2n_1 log2/128
156 // n log2/128 = n_7n_6n_5 log2/8 + n_4n_3n_2n_1 log2/128
157 // n log2/128 = I2 log2/8 + I1 log2/128
158 //
159 // N log2/128 = M log2 + I2 log2/8 + I1 log2/128
160 //
161 // exp(Z) = exp(s) (1+d) exp(log(2^M) + log(2^I2/8) + log(2^I1/128))
162 // exp(Z) = exp(s) f12 (2^M) 2^I2/8 2^I1/128
163 //
164 // I1, I2 are table indices. Use a series for exp(s).
165 // Then get exp(Z)
166 //
167 // exp(yLog(x)) = exp(Z) exp(Z3) f3
168 // exp(yLog(x)) = exp(Z)f3 exp(Z3)
169 // exp(yLog(x)) = A exp(Z3)
170 //
171 // We actually calculate exp(Z3) -1.
172 // Then,
173 // exp(yLog(x)) = A + A( exp(Z3) -1)
174 //
176 // Table Generation
177 //==============================================================
179 // The log values
180 // ==============
181 // The operation (K*log2_hi) must be exact. K is the true exponent of x.
182 // If we allow gradual underflow (denormals), K can be represented in 12 bits
183 // (as a two's complement number). We assume 13 bits as an engineering
184 // precaution.
185 //
186 // +------------+----------------+-+
187 // | 13 bits | 50 bits | |
188 // +------------+----------------+-+
189 // 0 1 66
190 // 2 34
191 //
192 // So we want the lsb(log2_hi) to be 2^-50
193 // We get log2 as a quad-extended (15-bit exponent, 128-bit significand)
194 //
195 // 0 fffe b17217f7d1cf79ab c9e3b39803f2f6af (4...)
196 //
197 // Consider numbering the bits left to right, starting at 0 thru 127.
198 // Bit 0 is the 2^-1 bit; bit 49 is the 2^-50 bit.
199 //
200 // ...79ab
201 // 0111 1001 1010 1011
202 // 44
203 // 89
204 //
205 // So if we shift off the rightmost 14 bits, then (shift back only
206 // the top half) we get
207 //
208 // 0 fffe b17217f7d1cf4000 e6af278ece600fcb dabc000000000000
209 //
210 // Put the right 64-bit signficand in an FR register, convert to double;
211 // it is exact. Put the next 128 bits into a quad register and round to double.
212 // The true exponent of the low part is -51.
213 //
214 // hi is 0 fffe b17217f7d1cf4000
215 // lo is 0 ffcc e6af278ece601000
216 //
217 // Convert to double memory format and get
218 //
219 // hi is 0x3fe62e42fefa39e8
220 // lo is 0x3cccd5e4f1d9cc02
221 //
222 // log2_hi + log2_lo is an accurate value for log2.
223 //
224 //
225 // The T and t values
226 // ==================
227 // A similar method is used to generate the T and t values.
228 //
229 // K * log2_hi + T must be exact.
230 //
231 // Smallest T,t
232 // ----------
233 // The smallest T,t is
234 // T t
235 // 0x3f60040155d58800, 0x3c93bce0ce3ddd81 log(1/frcpa(1+0/256))= +1.95503e-003
236 //
237 // The exponent is 0x3f6 (biased) or -9 (true).
238 // For the smallest T value, what we want is to clip the significand such that
239 // when it is shifted right by 9, its lsb is in the bit for 2^-51. The 9 is the
240 // specific for the first entry. In general, it is 0xffff - (biased 15-bit
241 // exponent).
243 // Independently, what we have calculated is the table value as a quad
244 // precision number.
245 // Table entry 1 is
246 // 0 fff6 80200aaeac44ef38 338f77605fdf8000
247 //
248 // We store this quad precision number in a data structure that is
249 // sign: 1
250 // exponent: 15
251 // signficand_hi: 64 (includes explicit bit)
252 // signficand_lo: 49
253 // Because the explicit bit is included, the significand is 113 bits.
254 //
255 // Consider significand_hi for table entry 1.
256 //
257 //
258 // +-+--- ... -------+--------------------+
259 // | |
260 // +-+--- ... -------+--------------------+
261 // 0 1 4444444455555555556666
262 // 2345678901234567890123
263 //
264 // Labeled as above, bit 0 is 2^0, bit 1 is 2^-1, etc.
265 // Bit 42 is 2^-42. If we shift to the right by 9, the bit in
266 // bit 42 goes in 51.
267 //
268 // So what we want to do is shift bits 43 thru 63 into significand_lo.
269 // This is shifting bit 42 into bit 63, taking care to retain shifted-off bits.
270 // Then shifting (just with signficaand_hi) back into bit 42.
271 //
272 // The shift_value is 63-42 = 21. In general, this is
273 // 63 - (51 -(0xffff - 0xfff6))
274 // For this example, it is
275 // 63 - (51 - 9) = 63 - 42 = 21
276 //
277 // This means we are shifting 21 bits into significand_lo. We must maintain more
278 // that a 128-bit signficand not to lose bits. So before the shift we put the
279 // 128-bit significand into a 256-bit signficand and then shift.
280 // The 256-bit significand has four parts: hh, hl, lh, and ll.
281 //
282 // Start off with
283 // hh hl lh ll
284 // <64> <49><15_0> <64_0> <64_0>
285 //
286 // After shift by 21 (then return for significand_hi),
287 // <43><21_0> <21><43> <6><58_0> <64_0>
288 //
289 // Take the hh part and convert to a double. There is no rounding here.
290 // The conversion is exact. The true exponent of the high part is the same as
291 // the true exponent of the input quad.
292 //
293 // We have some 64 plus significand bits for the low part. In this example, we
294 // have 70 bits. We want to round this to a double. Put them in a quad and then
295 // do a quad fnorm.
296 // For this example the true exponent of the low part is
297 // true_exponent_of_high - 43 = true_exponent_of_high - (64-21)
298 // In general, this is
299 // true_exponent_of_high - (64 - shift_value)
300 //
301 //
302 // Largest T,t
303 // ----------
304 // The largest T,t is
305 // 0x3fe62643fecf9742, 0x3c9e3147684bd37d log(1/frcpa(1+255/256))=+6.92171e-001
306 //
307 // Table entry 256 is
308 // 0 fffe b1321ff67cba178c 51da12f4df5a0000
309 //
310 // The shift value is
311 // 63 - (51 -(0xffff - 0xfffe)) = 13
312 //
313 // The true exponent of the low part is
314 // true_exponent_of_high - (64 - shift_value)
315 // -1 - (64-13) = -52
316 // Biased as a double, this is 0x3cb
317 //
318 //
319 //
320 // So then lsb(T) must be >= 2^-51
321 // msb(Klog2_hi) <= 2^12
322 //
323 // +--------+---------+
324 // | 51 bits | <== largest T
325 // +--------+---------+
326 // | 9 bits | 42 bits | <== smallest T
327 // +------------+----------------+-+
328 // | 13 bits | 50 bits | |
329 // +------------+----------------+-+
330 //
331 // Note: For powf only the table of T is needed
334 // Special Cases
335 //==============================================================
337 // double float
338 // overflow error 24 30
340 // underflow error 25 31
342 // X zero Y zero
343 // +0 +0 +1 error 26 32
344 // -0 +0 +1 error 26 32
345 // +0 -0 +1 error 26 32
346 // -0 -0 +1 error 26 32
348 // X zero Y negative
349 // +0 -odd integer +inf error 27 33 divide-by-zero
350 // -0 -odd integer -inf error 27 33 divide-by-zero
351 // +0 !-odd integer +inf error 27 33 divide-by-zero
352 // -0 !-odd integer +inf error 27 33 divide-by-zero
353 // +0 -inf +inf error 27 33 divide-by-zero
354 // -0 -inf +inf error 27 33 divide-by-zero
356 // X zero Y positive
357 // +0 +odd integer +0
358 // -0 +odd integer -0
359 // +0 !+odd integer +0
360 // -0 !+odd integer +0
361 // +0 +inf +0
362 // -0 +inf +0
363 // +0 Y NaN quiet Y invalid if Y SNaN
364 // -0 Y NaN quiet Y invalid if Y SNaN
366 // X one
367 // -1 Y inf +1
368 // -1 Y NaN quiet Y invalid if Y SNaN
369 // +1 Y NaN +1 invalid if Y SNaN
370 // +1 Y any else +1
372 // X - Y not integer QNAN error 28 34 invalid
374 // X NaN Y 0 +1 error 29 35
375 // X NaN Y NaN quiet X invalid if X or Y SNaN
376 // X NaN Y any else quiet X invalid if X SNaN
377 // X !+1 Y NaN quiet Y invalid if Y SNaN
380 // X +inf Y >0 +inf
381 // X -inf Y >0, !odd integer +inf
382 // X -inf Y >0, odd integer -inf
384 // X +inf Y <0 +0
385 // X -inf Y <0, !odd integer +0
386 // X -inf Y <0, odd integer -0
388 // X +inf Y =0 +1
389 // X -inf Y =0 +1
391 // |X|<1 Y +inf +0
392 // |X|<1 Y -inf +inf
393 // |X|>1 Y +inf +inf
394 // |X|>1 Y -inf +0
396 // X any Y =0 +1
398 // Assembly macros
399 //==============================================================
401 // integer registers used
403 pow_GR_exp_half = r10
404 pow_GR_signexp_Xm1 = r11
405 pow_GR_tmp = r11
407 pow_GR_signexp_X = r14
408 pow_GR_17ones = r15
409 pow_GR_Fpsr = r15
410 pow_AD_P = r16
411 pow_GR_rcs0_mask = r16
412 pow_GR_exp_2tom8 = r17
413 pow_GR_rcs0 = r17
414 pow_GR_sig_X = r18
415 pow_GR_10033 = r19
416 pow_GR_16ones = r20
418 pow_AD_Tt = r21
419 pow_GR_exp_X = r22
420 pow_AD_Q = r23
421 pow_GR_true_exp_X = r24
422 pow_GR_y_zero = r25
424 pow_GR_exp_Y = r26
425 pow_AD_tbl1 = r27
426 pow_AD_tbl2 = r28
427 pow_GR_offset = r29
428 pow_GR_exp_Xm1 = r30
429 pow_GR_xneg_yodd = r31
431 pow_GR_int_N = r38
432 pow_GR_index1 = r39
433 pow_GR_index2 = r40
435 pow_AD_T1 = r41
436 pow_AD_T2 = r42
437 pow_int_GR_M = r43
438 pow_GR_sig_int_Y = r44
439 pow_GR_sign_Y_Gpr = r45
441 pow_GR_17ones_m1 = r46
442 pow_GR_one = r47
443 pow_GR_sign_Y = r48
444 pow_GR_signexp_Y_Gpr = r49
445 pow_GR_exp_Y_Gpr = r50
447 pow_GR_true_exp_Y_Gpr = r51
448 pow_GR_signexp_Y = r52
449 pow_GR_x_one = r53
450 pow_GR_big_pos = r55
452 pow_GR_big_neg = r56
454 GR_SAVE_B0 = r50
455 GR_SAVE_GP = r51
456 GR_SAVE_PFS = r52
458 GR_Parameter_X = r53
459 GR_Parameter_Y = r54
460 GR_Parameter_RESULT = r55
461 pow_GR_tag = r56
464 // floating point registers used
466 POW_B = f32
467 POW_NORM_X = f33
468 POW_Xm1 = f34
469 POW_r1 = f34
471 POW_NORM_Y = f37
472 POW_Q2 = f38
473 POW_eps = f39
474 POW_P2 = f40
476 POW_P0 = f42
477 POW_log2_lo = f43
478 POW_r = f44
479 POW_Q0_half = f45
481 POW_tmp = f47
482 POW_log2_hi = f48
483 POW_Q1 = f49
484 POW_P1 = f50
486 POW_log2_by_128_hi = f51
487 POW_inv_log2_by_128 = f52
488 POW_rsq = f53
489 POW_Yrcub = f54
490 POW_log2_by_128_lo = f55
492 POW_xsq = f57
493 POW_v2 = f59
494 POW_T = f60
496 POW_RSHF = f62
497 POW_v210 = f63
498 POW_twoV = f65
500 POW_U = f66
501 POW_G = f67
502 POW_delta = f68
503 POW_V = f70
505 POW_p = f71
506 POW_Z = f72
507 POW_e3 = f73
508 POW_Z2 = f75
510 POW_W1 = f77
511 POW_Z3 = f80
513 POW_Z3sq = f85
515 POW_Nfloat = f87
516 POW_f3 = f89
517 POW_q = f90
519 POW_T1 = f96
520 POW_T2 = f97
521 POW_2M = f98
522 POW_s = f99
523 POW_f12 = f100
525 POW_ssq = f101
526 POW_T1T2 = f102
527 POW_1ps = f103
528 POW_A = f104
529 POW_es = f105
531 POW_Xp1 = f106
532 POW_int_K = f107
533 POW_K = f108
534 POW_f123 = f109
535 POW_Gpr = f110
537 POW_Y_Gpr = f111
538 POW_int_Y = f112
539 POW_2Mqp1 = f113
541 POW_float_int_Y = f116
542 POW_ftz_urm_f8 = f117
543 POW_wre_urm_f8 = f118
544 POW_big_neg = f119
545 POW_big_pos = f120
547 // Data tables
548 //==============================================================
550 RODATA
552 .align 16
554 LOCAL_OBJECT_START(pow_table_P)
555 data8 0x80000000000018E5, 0x0000BFFD // P_1
556 data8 0xb8aa3b295c17f0bc, 0x00004006 // inv_ln2_by_128
557 //
558 //
559 data8 0x3FA5555555554A9E // Q_2
560 data8 0x0000000000000000 // Pad
561 data8 0x3FC5555555554733 // Q_1
562 data8 0x43e8000000000000 // Right shift constant for exp
563 data8 0xc9e3b39803f2f6af, 0x00003fb7 // ln2_by_128_lo
564 LOCAL_OBJECT_END(pow_table_P)
566 LOCAL_OBJECT_START(pow_table_Q)
567 data8 0xCCCCCCCC4ED2BA7F, 0x00003FFC // P_2
568 data8 0xAAAAAAAAAAAAB505, 0x00003FFD // P_0
569 data8 0x3fe62e42fefa39e8, 0x3cccd5e4f1d9cc02 // log2 hi lo = +6.93147e-001
570 data8 0xb17217f7d1cf79ab, 0x00003ff7 // ln2_by_128_hi
571 LOCAL_OBJECT_END(pow_table_Q)
574 LOCAL_OBJECT_START(pow_Tt)
575 data8 0x3f60040155d58800 // log(1/frcpa(1+0/256))= +1.95503e-003
576 data8 0x3f78121214586a00 // log(1/frcpa(1+1/256))= +5.87661e-003
577 data8 0x3f841929f9683200 // log(1/frcpa(1+2/256))= +9.81362e-003
578 data8 0x3f8c317384c75f00 // log(1/frcpa(1+3/256))= +1.37662e-002
579 data8 0x3f91a6b91ac73380 // log(1/frcpa(1+4/256))= +1.72376e-002
580 data8 0x3f95ba9a5d9ac000 // log(1/frcpa(1+5/256))= +2.12196e-002
581 data8 0x3f99d2a807432580 // log(1/frcpa(1+6/256))= +2.52177e-002
582 data8 0x3f9d6b2725979800 // log(1/frcpa(1+7/256))= +2.87291e-002
583 data8 0x3fa0c58fa19dfa80 // log(1/frcpa(1+8/256))= +3.27573e-002
584 data8 0x3fa2954c78cbce00 // log(1/frcpa(1+9/256))= +3.62953e-002
585 data8 0x3fa4a94d2da96c40 // log(1/frcpa(1+10/256))= +4.03542e-002
586 data8 0x3fa67c94f2d4bb40 // log(1/frcpa(1+11/256))= +4.39192e-002
587 data8 0x3fa85188b630f040 // log(1/frcpa(1+12/256))= +4.74971e-002
588 data8 0x3faa6b8abe73af40 // log(1/frcpa(1+13/256))= +5.16017e-002
589 data8 0x3fac441e06f72a80 // log(1/frcpa(1+14/256))= +5.52072e-002
590 data8 0x3fae1e6713606d00 // log(1/frcpa(1+15/256))= +5.88257e-002
591 data8 0x3faffa6911ab9300 // log(1/frcpa(1+16/256))= +6.24574e-002
592 data8 0x3fb0ec139c5da600 // log(1/frcpa(1+17/256))= +6.61022e-002
593 data8 0x3fb1dbd2643d1900 // log(1/frcpa(1+18/256))= +6.97605e-002
594 data8 0x3fb2cc7284fe5f00 // log(1/frcpa(1+19/256))= +7.34321e-002
595 data8 0x3fb3bdf5a7d1ee60 // log(1/frcpa(1+20/256))= +7.71173e-002
596 data8 0x3fb4b05d7aa012e0 // log(1/frcpa(1+21/256))= +8.08161e-002
597 data8 0x3fb580db7ceb5700 // log(1/frcpa(1+22/256))= +8.39975e-002
598 data8 0x3fb674f089365a60 // log(1/frcpa(1+23/256))= +8.77219e-002
599 data8 0x3fb769ef2c6b5680 // log(1/frcpa(1+24/256))= +9.14602e-002
600 data8 0x3fb85fd927506a40 // log(1/frcpa(1+25/256))= +9.52125e-002
601 data8 0x3fb9335e5d594980 // log(1/frcpa(1+26/256))= +9.84401e-002
602 data8 0x3fba2b0220c8e5e0 // log(1/frcpa(1+27/256))= +1.02219e-001
603 data8 0x3fbb0004ac1a86a0 // log(1/frcpa(1+28/256))= +1.05469e-001
604 data8 0x3fbbf968769fca00 // log(1/frcpa(1+29/256))= +1.09274e-001
605 data8 0x3fbccfedbfee13a0 // log(1/frcpa(1+30/256))= +1.12548e-001
606 data8 0x3fbda727638446a0 // log(1/frcpa(1+31/256))= +1.15832e-001
607 data8 0x3fbea3257fe10f60 // log(1/frcpa(1+32/256))= +1.19677e-001
608 data8 0x3fbf7be9fedbfde0 // log(1/frcpa(1+33/256))= +1.22985e-001
609 data8 0x3fc02ab352ff25f0 // log(1/frcpa(1+34/256))= +1.26303e-001
610 data8 0x3fc097ce579d2040 // log(1/frcpa(1+35/256))= +1.29633e-001
611 data8 0x3fc1178e8227e470 // log(1/frcpa(1+36/256))= +1.33531e-001
612 data8 0x3fc185747dbecf30 // log(1/frcpa(1+37/256))= +1.36885e-001
613 data8 0x3fc1f3b925f25d40 // log(1/frcpa(1+38/256))= +1.40250e-001
614 data8 0x3fc2625d1e6ddf50 // log(1/frcpa(1+39/256))= +1.43627e-001
615 data8 0x3fc2d1610c868130 // log(1/frcpa(1+40/256))= +1.47015e-001
616 data8 0x3fc340c597411420 // log(1/frcpa(1+41/256))= +1.50414e-001
617 data8 0x3fc3b08b6757f2a0 // log(1/frcpa(1+42/256))= +1.53825e-001
618 data8 0x3fc40dfb08378000 // log(1/frcpa(1+43/256))= +1.56677e-001
619 data8 0x3fc47e74e8ca5f70 // log(1/frcpa(1+44/256))= +1.60109e-001
620 data8 0x3fc4ef51f6466de0 // log(1/frcpa(1+45/256))= +1.63553e-001
621 data8 0x3fc56092e02ba510 // log(1/frcpa(1+46/256))= +1.67010e-001
622 data8 0x3fc5d23857cd74d0 // log(1/frcpa(1+47/256))= +1.70478e-001
623 data8 0x3fc6313a37335d70 // log(1/frcpa(1+48/256))= +1.73377e-001
624 data8 0x3fc6a399dabbd380 // log(1/frcpa(1+49/256))= +1.76868e-001
625 data8 0x3fc70337dd3ce410 // log(1/frcpa(1+50/256))= +1.79786e-001
626 data8 0x3fc77654128f6120 // log(1/frcpa(1+51/256))= +1.83299e-001
627 data8 0x3fc7e9d82a0b0220 // log(1/frcpa(1+52/256))= +1.86824e-001
628 data8 0x3fc84a6b759f5120 // log(1/frcpa(1+53/256))= +1.89771e-001
629 data8 0x3fc8ab47d5f5a300 // log(1/frcpa(1+54/256))= +1.92727e-001
630 data8 0x3fc91fe490965810 // log(1/frcpa(1+55/256))= +1.96286e-001
631 data8 0x3fc981634011aa70 // log(1/frcpa(1+56/256))= +1.99261e-001
632 data8 0x3fc9f6c407089660 // log(1/frcpa(1+57/256))= +2.02843e-001
633 data8 0x3fca58e729348f40 // log(1/frcpa(1+58/256))= +2.05838e-001
634 data8 0x3fcabb55c31693a0 // log(1/frcpa(1+59/256))= +2.08842e-001
635 data8 0x3fcb1e104919efd0 // log(1/frcpa(1+60/256))= +2.11855e-001
636 data8 0x3fcb94ee93e367c0 // log(1/frcpa(1+61/256))= +2.15483e-001
637 data8 0x3fcbf851c0675550 // log(1/frcpa(1+62/256))= +2.18516e-001
638 data8 0x3fcc5c0254bf23a0 // log(1/frcpa(1+63/256))= +2.21558e-001
639 data8 0x3fccc000c9db3c50 // log(1/frcpa(1+64/256))= +2.24609e-001
640 data8 0x3fcd244d99c85670 // log(1/frcpa(1+65/256))= +2.27670e-001
641 data8 0x3fcd88e93fb2f450 // log(1/frcpa(1+66/256))= +2.30741e-001
642 data8 0x3fcdedd437eaef00 // log(1/frcpa(1+67/256))= +2.33820e-001
643 data8 0x3fce530effe71010 // log(1/frcpa(1+68/256))= +2.36910e-001
644 data8 0x3fceb89a1648b970 // log(1/frcpa(1+69/256))= +2.40009e-001
645 data8 0x3fcf1e75fadf9bd0 // log(1/frcpa(1+70/256))= +2.43117e-001
646 data8 0x3fcf84a32ead7c30 // log(1/frcpa(1+71/256))= +2.46235e-001
647 data8 0x3fcfeb2233ea07c0 // log(1/frcpa(1+72/256))= +2.49363e-001
648 data8 0x3fd028f9c7035c18 // log(1/frcpa(1+73/256))= +2.52501e-001
649 data8 0x3fd05c8be0d96358 // log(1/frcpa(1+74/256))= +2.55649e-001
650 data8 0x3fd085eb8f8ae790 // log(1/frcpa(1+75/256))= +2.58174e-001
651 data8 0x3fd0b9c8e32d1910 // log(1/frcpa(1+76/256))= +2.61339e-001
652 data8 0x3fd0edd060b78080 // log(1/frcpa(1+77/256))= +2.64515e-001
653 data8 0x3fd122024cf00638 // log(1/frcpa(1+78/256))= +2.67701e-001
654 data8 0x3fd14be2927aecd0 // log(1/frcpa(1+79/256))= +2.70257e-001
655 data8 0x3fd180618ef18ad8 // log(1/frcpa(1+80/256))= +2.73461e-001
656 data8 0x3fd1b50bbe2fc638 // log(1/frcpa(1+81/256))= +2.76675e-001
657 data8 0x3fd1df4cc7cf2428 // log(1/frcpa(1+82/256))= +2.79254e-001
658 data8 0x3fd214456d0eb8d0 // log(1/frcpa(1+83/256))= +2.82487e-001
659 data8 0x3fd23ec5991eba48 // log(1/frcpa(1+84/256))= +2.85081e-001
660 data8 0x3fd2740d9f870af8 // log(1/frcpa(1+85/256))= +2.88333e-001
661 data8 0x3fd29ecdabcdfa00 // log(1/frcpa(1+86/256))= +2.90943e-001
662 data8 0x3fd2d46602adcce8 // log(1/frcpa(1+87/256))= +2.94214e-001
663 data8 0x3fd2ff66b04ea9d0 // log(1/frcpa(1+88/256))= +2.96838e-001
664 data8 0x3fd335504b355a30 // log(1/frcpa(1+89/256))= +3.00129e-001
665 data8 0x3fd360925ec44f58 // log(1/frcpa(1+90/256))= +3.02769e-001
666 data8 0x3fd38bf1c3337e70 // log(1/frcpa(1+91/256))= +3.05417e-001
667 data8 0x3fd3c25277333180 // log(1/frcpa(1+92/256))= +3.08735e-001
668 data8 0x3fd3edf463c16838 // log(1/frcpa(1+93/256))= +3.11399e-001
669 data8 0x3fd419b423d5e8c0 // log(1/frcpa(1+94/256))= +3.14069e-001
670 data8 0x3fd44591e0539f48 // log(1/frcpa(1+95/256))= +3.16746e-001
671 data8 0x3fd47c9175b6f0a8 // log(1/frcpa(1+96/256))= +3.20103e-001
672 data8 0x3fd4a8b341552b08 // log(1/frcpa(1+97/256))= +3.22797e-001
673 data8 0x3fd4d4f390890198 // log(1/frcpa(1+98/256))= +3.25498e-001
674 data8 0x3fd501528da1f960 // log(1/frcpa(1+99/256))= +3.28206e-001
675 data8 0x3fd52dd06347d4f0 // log(1/frcpa(1+100/256))= +3.30921e-001
676 data8 0x3fd55a6d3c7b8a88 // log(1/frcpa(1+101/256))= +3.33644e-001
677 data8 0x3fd5925d2b112a58 // log(1/frcpa(1+102/256))= +3.37058e-001
678 data8 0x3fd5bf406b543db0 // log(1/frcpa(1+103/256))= +3.39798e-001
679 data8 0x3fd5ec433d5c35a8 // log(1/frcpa(1+104/256))= +3.42545e-001
680 data8 0x3fd61965cdb02c18 // log(1/frcpa(1+105/256))= +3.45300e-001
681 data8 0x3fd646a84935b2a0 // log(1/frcpa(1+106/256))= +3.48063e-001
682 data8 0x3fd6740add31de90 // log(1/frcpa(1+107/256))= +3.50833e-001
683 data8 0x3fd6a18db74a58c0 // log(1/frcpa(1+108/256))= +3.53610e-001
684 data8 0x3fd6cf31058670e8 // log(1/frcpa(1+109/256))= +3.56396e-001
685 data8 0x3fd6f180e852f0b8 // log(1/frcpa(1+110/256))= +3.58490e-001
686 data8 0x3fd71f5d71b894e8 // log(1/frcpa(1+111/256))= +3.61289e-001
687 data8 0x3fd74d5aefd66d58 // log(1/frcpa(1+112/256))= +3.64096e-001
688 data8 0x3fd77b79922bd378 // log(1/frcpa(1+113/256))= +3.66911e-001
689 data8 0x3fd7a9b9889f19e0 // log(1/frcpa(1+114/256))= +3.69734e-001
690 data8 0x3fd7d81b037eb6a0 // log(1/frcpa(1+115/256))= +3.72565e-001
691 data8 0x3fd8069e33827230 // log(1/frcpa(1+116/256))= +3.75404e-001
692 data8 0x3fd82996d3ef8bc8 // log(1/frcpa(1+117/256))= +3.77538e-001
693 data8 0x3fd85855776dcbf8 // log(1/frcpa(1+118/256))= +3.80391e-001
694 data8 0x3fd8873658327cc8 // log(1/frcpa(1+119/256))= +3.83253e-001
695 data8 0x3fd8aa75973ab8c8 // log(1/frcpa(1+120/256))= +3.85404e-001
696 data8 0x3fd8d992dc8824e0 // log(1/frcpa(1+121/256))= +3.88280e-001
697 data8 0x3fd908d2ea7d9510 // log(1/frcpa(1+122/256))= +3.91164e-001
698 data8 0x3fd92c59e79c0e50 // log(1/frcpa(1+123/256))= +3.93332e-001
699 data8 0x3fd95bd750ee3ed0 // log(1/frcpa(1+124/256))= +3.96231e-001
700 data8 0x3fd98b7811a3ee58 // log(1/frcpa(1+125/256))= +3.99138e-001
701 data8 0x3fd9af47f33d4068 // log(1/frcpa(1+126/256))= +4.01323e-001
702 data8 0x3fd9df270c1914a0 // log(1/frcpa(1+127/256))= +4.04245e-001
703 data8 0x3fda0325ed14fda0 // log(1/frcpa(1+128/256))= +4.06442e-001
704 data8 0x3fda33440224fa78 // log(1/frcpa(1+129/256))= +4.09379e-001
705 data8 0x3fda57725e80c380 // log(1/frcpa(1+130/256))= +4.11587e-001
706 data8 0x3fda87d0165dd198 // log(1/frcpa(1+131/256))= +4.14539e-001
707 data8 0x3fdaac2e6c03f890 // log(1/frcpa(1+132/256))= +4.16759e-001
708 data8 0x3fdadccc6fdf6a80 // log(1/frcpa(1+133/256))= +4.19726e-001
709 data8 0x3fdb015b3eb1e790 // log(1/frcpa(1+134/256))= +4.21958e-001
710 data8 0x3fdb323a3a635948 // log(1/frcpa(1+135/256))= +4.24941e-001
711 data8 0x3fdb56fa04462908 // log(1/frcpa(1+136/256))= +4.27184e-001
712 data8 0x3fdb881aa659bc90 // log(1/frcpa(1+137/256))= +4.30182e-001
713 data8 0x3fdbad0bef3db160 // log(1/frcpa(1+138/256))= +4.32437e-001
714 data8 0x3fdbd21297781c28 // log(1/frcpa(1+139/256))= +4.34697e-001
715 data8 0x3fdc039236f08818 // log(1/frcpa(1+140/256))= +4.37718e-001
716 data8 0x3fdc28cb1e4d32f8 // log(1/frcpa(1+141/256))= +4.39990e-001
717 data8 0x3fdc4e19b84723c0 // log(1/frcpa(1+142/256))= +4.42267e-001
718 data8 0x3fdc7ff9c74554c8 // log(1/frcpa(1+143/256))= +4.45311e-001
719 data8 0x3fdca57b64e9db00 // log(1/frcpa(1+144/256))= +4.47600e-001
720 data8 0x3fdccb130a5ceba8 // log(1/frcpa(1+145/256))= +4.49895e-001
721 data8 0x3fdcf0c0d18f3268 // log(1/frcpa(1+146/256))= +4.52194e-001
722 data8 0x3fdd232075b5a200 // log(1/frcpa(1+147/256))= +4.55269e-001
723 data8 0x3fdd490246defa68 // log(1/frcpa(1+148/256))= +4.57581e-001
724 data8 0x3fdd6efa918d25c8 // log(1/frcpa(1+149/256))= +4.59899e-001
725 data8 0x3fdd9509707ae528 // log(1/frcpa(1+150/256))= +4.62221e-001
726 data8 0x3fddbb2efe92c550 // log(1/frcpa(1+151/256))= +4.64550e-001
727 data8 0x3fddee2f3445e4a8 // log(1/frcpa(1+152/256))= +4.67663e-001
728 data8 0x3fde148a1a2726c8 // log(1/frcpa(1+153/256))= +4.70004e-001
729 data8 0x3fde3afc0a49ff38 // log(1/frcpa(1+154/256))= +4.72350e-001
730 data8 0x3fde6185206d5168 // log(1/frcpa(1+155/256))= +4.74702e-001
731 data8 0x3fde882578823d50 // log(1/frcpa(1+156/256))= +4.77060e-001
732 data8 0x3fdeaedd2eac9908 // log(1/frcpa(1+157/256))= +4.79423e-001
733 data8 0x3fded5ac5f436be0 // log(1/frcpa(1+158/256))= +4.81792e-001
734 data8 0x3fdefc9326d16ab8 // log(1/frcpa(1+159/256))= +4.84166e-001
735 data8 0x3fdf2391a21575f8 // log(1/frcpa(1+160/256))= +4.86546e-001
736 data8 0x3fdf4aa7ee031928 // log(1/frcpa(1+161/256))= +4.88932e-001
737 data8 0x3fdf71d627c30bb0 // log(1/frcpa(1+162/256))= +4.91323e-001
738 data8 0x3fdf991c6cb3b378 // log(1/frcpa(1+163/256))= +4.93720e-001
739 data8 0x3fdfc07ada69a908 // log(1/frcpa(1+164/256))= +4.96123e-001
740 data8 0x3fdfe7f18eb03d38 // log(1/frcpa(1+165/256))= +4.98532e-001
741 data8 0x3fe007c053c5002c // log(1/frcpa(1+166/256))= +5.00946e-001
742 data8 0x3fe01b942198a5a0 // log(1/frcpa(1+167/256))= +5.03367e-001
743 data8 0x3fe02f74400c64e8 // log(1/frcpa(1+168/256))= +5.05793e-001
744 data8 0x3fe04360be7603ac // log(1/frcpa(1+169/256))= +5.08225e-001
745 data8 0x3fe05759ac47fe30 // log(1/frcpa(1+170/256))= +5.10663e-001
746 data8 0x3fe06b5f1911cf50 // log(1/frcpa(1+171/256))= +5.13107e-001
747 data8 0x3fe078bf0533c568 // log(1/frcpa(1+172/256))= +5.14740e-001
748 data8 0x3fe08cd9687e7b0c // log(1/frcpa(1+173/256))= +5.17194e-001
749 data8 0x3fe0a10074cf9018 // log(1/frcpa(1+174/256))= +5.19654e-001
750 data8 0x3fe0b5343a234474 // log(1/frcpa(1+175/256))= +5.22120e-001
751 data8 0x3fe0c974c89431cc // log(1/frcpa(1+176/256))= +5.24592e-001
752 data8 0x3fe0ddc2305b9884 // log(1/frcpa(1+177/256))= +5.27070e-001
753 data8 0x3fe0eb524bafc918 // log(1/frcpa(1+178/256))= +5.28726e-001
754 data8 0x3fe0ffb54213a474 // log(1/frcpa(1+179/256))= +5.31214e-001
755 data8 0x3fe114253da97d9c // log(1/frcpa(1+180/256))= +5.33709e-001
756 data8 0x3fe128a24f1d9afc // log(1/frcpa(1+181/256))= +5.36210e-001
757 data8 0x3fe1365252bf0864 // log(1/frcpa(1+182/256))= +5.37881e-001
758 data8 0x3fe14ae558b4a92c // log(1/frcpa(1+183/256))= +5.40393e-001
759 data8 0x3fe15f85a19c7658 // log(1/frcpa(1+184/256))= +5.42910e-001
760 data8 0x3fe16d4d38c119f8 // log(1/frcpa(1+185/256))= +5.44592e-001
761 data8 0x3fe18203c20dd130 // log(1/frcpa(1+186/256))= +5.47121e-001
762 data8 0x3fe196c7bc4b1f38 // log(1/frcpa(1+187/256))= +5.49656e-001
763 data8 0x3fe1a4a738b7a33c // log(1/frcpa(1+188/256))= +5.51349e-001
764 data8 0x3fe1b981c0c9653c // log(1/frcpa(1+189/256))= +5.53895e-001
765 data8 0x3fe1ce69e8bb1068 // log(1/frcpa(1+190/256))= +5.56447e-001
766 data8 0x3fe1dc619de06944 // log(1/frcpa(1+191/256))= +5.58152e-001
767 data8 0x3fe1f160a2ad0da0 // log(1/frcpa(1+192/256))= +5.60715e-001
768 data8 0x3fe2066d7740737c // log(1/frcpa(1+193/256))= +5.63285e-001
769 data8 0x3fe2147dba47a390 // log(1/frcpa(1+194/256))= +5.65001e-001
770 data8 0x3fe229a1bc5ebac0 // log(1/frcpa(1+195/256))= +5.67582e-001
771 data8 0x3fe237c1841a502c // log(1/frcpa(1+196/256))= +5.69306e-001
772 data8 0x3fe24cfce6f80d98 // log(1/frcpa(1+197/256))= +5.71898e-001
773 data8 0x3fe25b2c55cd5760 // log(1/frcpa(1+198/256))= +5.73630e-001
774 data8 0x3fe2707f4d5f7c40 // log(1/frcpa(1+199/256))= +5.76233e-001
775 data8 0x3fe285e0842ca380 // log(1/frcpa(1+200/256))= +5.78842e-001
776 data8 0x3fe294294708b770 // log(1/frcpa(1+201/256))= +5.80586e-001
777 data8 0x3fe2a9a2670aff0c // log(1/frcpa(1+202/256))= +5.83207e-001
778 data8 0x3fe2b7fb2c8d1cc0 // log(1/frcpa(1+203/256))= +5.84959e-001
779 data8 0x3fe2c65a6395f5f4 // log(1/frcpa(1+204/256))= +5.86713e-001
780 data8 0x3fe2dbf557b0df40 // log(1/frcpa(1+205/256))= +5.89350e-001
781 data8 0x3fe2ea64c3f97654 // log(1/frcpa(1+206/256))= +5.91113e-001
782 data8 0x3fe3001823684d70 // log(1/frcpa(1+207/256))= +5.93762e-001
783 data8 0x3fe30e97e9a8b5cc // log(1/frcpa(1+208/256))= +5.95531e-001
784 data8 0x3fe32463ebdd34e8 // log(1/frcpa(1+209/256))= +5.98192e-001
785 data8 0x3fe332f4314ad794 // log(1/frcpa(1+210/256))= +5.99970e-001
786 data8 0x3fe348d90e7464cc // log(1/frcpa(1+211/256))= +6.02643e-001
787 data8 0x3fe35779f8c43d6c // log(1/frcpa(1+212/256))= +6.04428e-001
788 data8 0x3fe36621961a6a98 // log(1/frcpa(1+213/256))= +6.06217e-001
789 data8 0x3fe37c299f3c3668 // log(1/frcpa(1+214/256))= +6.08907e-001
790 data8 0x3fe38ae2171976e4 // log(1/frcpa(1+215/256))= +6.10704e-001
791 data8 0x3fe399a157a603e4 // log(1/frcpa(1+216/256))= +6.12504e-001
792 data8 0x3fe3afccfe77b9d0 // log(1/frcpa(1+217/256))= +6.15210e-001
793 data8 0x3fe3be9d503533b4 // log(1/frcpa(1+218/256))= +6.17018e-001
794 data8 0x3fe3cd7480b4a8a0 // log(1/frcpa(1+219/256))= +6.18830e-001
795 data8 0x3fe3e3c43918f76c // log(1/frcpa(1+220/256))= +6.21554e-001
796 data8 0x3fe3f2acb27ed6c4 // log(1/frcpa(1+221/256))= +6.23373e-001
797 data8 0x3fe4019c2125ca90 // log(1/frcpa(1+222/256))= +6.25197e-001
798 data8 0x3fe4181061389720 // log(1/frcpa(1+223/256))= +6.27937e-001
799 data8 0x3fe42711518df544 // log(1/frcpa(1+224/256))= +6.29769e-001
800 data8 0x3fe436194e12b6bc // log(1/frcpa(1+225/256))= +6.31604e-001
801 data8 0x3fe445285d68ea68 // log(1/frcpa(1+226/256))= +6.33442e-001
802 data8 0x3fe45bcc464c8938 // log(1/frcpa(1+227/256))= +6.36206e-001
803 data8 0x3fe46aed21f117fc // log(1/frcpa(1+228/256))= +6.38053e-001
804 data8 0x3fe47a1527e8a2d0 // log(1/frcpa(1+229/256))= +6.39903e-001
805 data8 0x3fe489445efffcc8 // log(1/frcpa(1+230/256))= +6.41756e-001
806 data8 0x3fe4a018bcb69834 // log(1/frcpa(1+231/256))= +6.44543e-001
807 data8 0x3fe4af5a0c9d65d4 // log(1/frcpa(1+232/256))= +6.46405e-001
808 data8 0x3fe4bea2a5bdbe84 // log(1/frcpa(1+233/256))= +6.48271e-001
809 data8 0x3fe4cdf28f10ac44 // log(1/frcpa(1+234/256))= +6.50140e-001
810 data8 0x3fe4dd49cf994058 // log(1/frcpa(1+235/256))= +6.52013e-001
811 data8 0x3fe4eca86e64a680 // log(1/frcpa(1+236/256))= +6.53889e-001
812 data8 0x3fe503c43cd8eb68 // log(1/frcpa(1+237/256))= +6.56710e-001
813 data8 0x3fe513356667fc54 // log(1/frcpa(1+238/256))= +6.58595e-001
814 data8 0x3fe522ae0738a3d4 // log(1/frcpa(1+239/256))= +6.60483e-001
815 data8 0x3fe5322e26867854 // log(1/frcpa(1+240/256))= +6.62376e-001
816 data8 0x3fe541b5cb979808 // log(1/frcpa(1+241/256))= +6.64271e-001
817 data8 0x3fe55144fdbcbd60 // log(1/frcpa(1+242/256))= +6.66171e-001
818 data8 0x3fe560dbc45153c4 // log(1/frcpa(1+243/256))= +6.68074e-001
819 data8 0x3fe5707a26bb8c64 // log(1/frcpa(1+244/256))= +6.69980e-001
820 data8 0x3fe587f60ed5b8fc // log(1/frcpa(1+245/256))= +6.72847e-001
821 data8 0x3fe597a7977c8f30 // log(1/frcpa(1+246/256))= +6.74763e-001
822 data8 0x3fe5a760d634bb88 // log(1/frcpa(1+247/256))= +6.76682e-001
823 data8 0x3fe5b721d295f10c // log(1/frcpa(1+248/256))= +6.78605e-001
824 data8 0x3fe5c6ea94431ef8 // log(1/frcpa(1+249/256))= +6.80532e-001
825 data8 0x3fe5d6bb22ea86f4 // log(1/frcpa(1+250/256))= +6.82462e-001
826 data8 0x3fe5e6938645d38c // log(1/frcpa(1+251/256))= +6.84397e-001
827 data8 0x3fe5f673c61a2ed0 // log(1/frcpa(1+252/256))= +6.86335e-001
828 data8 0x3fe6065bea385924 // log(1/frcpa(1+253/256))= +6.88276e-001
829 data8 0x3fe6164bfa7cc068 // log(1/frcpa(1+254/256))= +6.90222e-001
830 data8 0x3fe62643fecf9740 // log(1/frcpa(1+255/256))= +6.92171e-001
831 LOCAL_OBJECT_END(pow_Tt)
834 // Table 1 is 2^(index_1/128) where
835 // index_1 goes from 0 to 15
836 LOCAL_OBJECT_START(pow_tbl1)
837 data8 0x8000000000000000 , 0x00003FFF
838 data8 0x80B1ED4FD999AB6C , 0x00003FFF
839 data8 0x8164D1F3BC030773 , 0x00003FFF
840 data8 0x8218AF4373FC25EC , 0x00003FFF
841 data8 0x82CD8698AC2BA1D7 , 0x00003FFF
842 data8 0x8383594EEFB6EE37 , 0x00003FFF
843 data8 0x843A28C3ACDE4046 , 0x00003FFF
844 data8 0x84F1F656379C1A29 , 0x00003FFF
845 data8 0x85AAC367CC487B15 , 0x00003FFF
846 data8 0x8664915B923FBA04 , 0x00003FFF
847 data8 0x871F61969E8D1010 , 0x00003FFF
848 data8 0x87DB357FF698D792 , 0x00003FFF
849 data8 0x88980E8092DA8527 , 0x00003FFF
850 data8 0x8955EE03618E5FDD , 0x00003FFF
851 data8 0x8A14D575496EFD9A , 0x00003FFF
852 data8 0x8AD4C6452C728924 , 0x00003FFF
853 LOCAL_OBJECT_END(pow_tbl1)
856 // Table 2 is 2^(index_1/8) where
857 // index_2 goes from 0 to 7
858 LOCAL_OBJECT_START(pow_tbl2)
859 data8 0x8000000000000000 , 0x00003FFF
860 data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF
861 data8 0x9837F0518DB8A96F , 0x00003FFF
862 data8 0xA5FED6A9B15138EA , 0x00003FFF
863 data8 0xB504F333F9DE6484 , 0x00003FFF
864 data8 0xC5672A115506DADD , 0x00003FFF
865 data8 0xD744FCCAD69D6AF4 , 0x00003FFF
866 data8 0xEAC0C6E7DD24392F , 0x00003FFF
867 LOCAL_OBJECT_END(pow_tbl2)
869 .section .text
870 WEAK_LIBM_ENTRY(powf)
872 // Get exponent of x. Will be used to calculate K.
873 { .mfi
874 getf.exp pow_GR_signexp_X = f8
875 fms.s1 POW_Xm1 = f8,f1,f1 // Will be used for r1 if x>0
876 mov pow_GR_17ones = 0x1FFFF
877 }
878 { .mfi
879 addl pow_AD_P = @ltoff(pow_table_P), gp
880 fma.s1 POW_Xp1 = f8,f1,f1 // Will be used for r1 if x<0
881 nop.i 999
882 }
883 ;;
885 // Get significand of x. Will be used to get index to fetch T, Tt.
886 { .mfi
887 getf.sig pow_GR_sig_X = f8
888 frcpa.s1 POW_B, p6 = f1,f8
889 mov pow_GR_exp_half = 0xFFFE // Exponent for 0.5
890 }
891 { .mfi
892 ld8 pow_AD_P = [pow_AD_P]
893 fma.s1 POW_NORM_X = f8,f1,f0
894 mov pow_GR_exp_2tom8 = 0xFFF7
895 }
896 ;;
898 // DOUBLE 0x10033 exponent limit at which y is an integer
899 { .mfi
900 nop.m 999
901 fcmp.lt.s1 p8,p9 = f8, f0 // Test for x<0
902 addl pow_GR_10033 = 0x10033, r0
903 }
904 { .mfi
905 mov pow_GR_16ones = 0xFFFF
906 fma.s1 POW_NORM_Y = f9,f1,f0
907 nop.i 999
908 }
909 ;;
911 // p13 = TRUE ==> X is unorm
912 { .mfi
913 setf.exp POW_Q0_half = pow_GR_exp_half // Form 0.5
914 fclass.m p13,p0 = f8, 0x0b // Test for x unorm
915 adds pow_AD_Tt = pow_Tt - pow_table_P, pow_AD_P
916 }
917 { .mfi
918 adds pow_AD_Q = pow_table_Q - pow_table_P, pow_AD_P
919 nop.f 999
920 nop.i 999
921 }
922 ;;
924 // p14 = TRUE ==> X is ZERO
925 { .mfi
926 ldfe POW_P2 = [pow_AD_Q], 16
927 fclass.m p14,p0 = f8, 0x07
928 nop.i 999
929 }
930 // Note POW_Xm1 and POW_r1 are used interchangeably
931 { .mfb
932 nop.m 999
933 (p8) fnma.s1 POW_Xm1 = POW_Xp1,f1,f0
934 (p13) br.cond.spnt POW_X_DENORM
935 }
936 ;;
938 // Continue normal and denormal paths here
939 POW_COMMON:
940 // p11 = TRUE ==> Y is a NAN
941 { .mfi
942 and pow_GR_exp_X = pow_GR_signexp_X, pow_GR_17ones
943 fclass.m p11,p0 = f9, 0xc3
944 nop.i 999
945 }
946 { .mfi
947 nop.m 999
948 fms.s1 POW_r = POW_B, POW_NORM_X,f1
949 mov pow_GR_y_zero = 0
950 }
951 ;;
953 // Get exponent of |x|-1 to use in comparison to 2^-8
954 { .mmi
955 getf.exp pow_GR_signexp_Xm1 = POW_Xm1
956 sub pow_GR_true_exp_X = pow_GR_exp_X, pow_GR_16ones
957 extr.u pow_GR_offset = pow_GR_sig_X, 55, 8
958 }
959 ;;
961 { .mfi
962 alloc r32=ar.pfs,2,19,4,0
963 fcvt.fx.s1 POW_int_Y = POW_NORM_Y
964 shladd pow_AD_Tt = pow_GR_offset, 3, pow_AD_Tt
965 }
966 { .mfi
967 setf.sig POW_int_K = pow_GR_true_exp_X
968 nop.f 999
969 nop.i 999
970 }
971 ;;
973 // p12 = TRUE if Y is ZERO
974 // Compute xsq to decide later if |x|=1
975 { .mfi
976 ldfe POW_P1 = [pow_AD_P], 16
977 fclass.m p12,p0 = f9, 0x07
978 nop.i 999
979 }
980 { .mfb
981 ldfe POW_P0 = [pow_AD_Q], 16
982 fma.s1 POW_xsq = POW_NORM_X, POW_NORM_X, f0
983 (p11) br.cond.spnt POW_Y_NAN // Branch if y=nan
984 }
985 ;;
987 { .mmf
988 getf.exp pow_GR_signexp_Y = POW_NORM_Y
989 ldfd POW_T = [pow_AD_Tt]
990 fma.s1 POW_rsq = POW_r, POW_r,f0
991 }
992 ;;
994 // p11 = TRUE ==> X is a NAN
995 { .mfi
996 ldfpd POW_log2_hi, POW_log2_lo = [pow_AD_Q], 16
997 fclass.m p11,p0 = POW_NORM_X, 0xc3
998 nop.i 999
999 }
1000 { .mfi
1001 ldfe POW_inv_log2_by_128 = [pow_AD_P], 16
1002 fma.s1 POW_delta = f0,f0,f0 // delta=0 in case |x| near 1
1003 (p12) mov pow_GR_y_zero = 1
1004 }
1005 ;;
1007 { .mfi
1008 ldfd POW_Q2 = [pow_AD_P], 16
1009 fnma.s1 POW_twoV = POW_r, POW_Q0_half,f1
1010 and pow_GR_exp_Xm1 = pow_GR_signexp_Xm1, pow_GR_17ones
1011 }
1012 { .mfi
1013 nop.m 999
1014 fma.s1 POW_U = POW_NORM_Y,POW_r,f0
1015 nop.i 999
1016 }
1017 ;;
1019 // Determine if we will use the |x| near 1 path (p6) or normal path (p7)
1020 { .mfi
1021 nop.m 999
1022 fcvt.xf POW_K = POW_int_K
1023 cmp.lt p6,p7 = pow_GR_exp_Xm1, pow_GR_exp_2tom8
1024 }
1025 { .mfb
1026 nop.m 999
1027 fma.s1 POW_G = f0,f0,f0 // G=0 in case |x| near 1
1028 (p11) br.cond.spnt POW_X_NAN // Branch if x=nan and y not nan
1029 }
1030 ;;
1032 // If on the x near 1 path, assign r1 to r
1033 { .mfi
1034 ldfpd POW_Q1, POW_RSHF = [pow_AD_P], 16
1035 (p6) fma.s1 POW_r = POW_r1, f1, f0
1036 nop.i 999
1037 }
1038 { .mfb
1039 nop.m 999
1040 (p6) fma.s1 POW_rsq = POW_r1, POW_r1, f0
1041 (p14) br.cond.spnt POW_X_0 // Branch if x zero and y not nan
1042 }
1043 ;;
1045 { .mfi
1046 getf.sig pow_GR_sig_int_Y = POW_int_Y
1047 (p6) fnma.s1 POW_twoV = POW_r1, POW_Q0_half,f1
1048 and pow_GR_exp_Y = pow_GR_signexp_Y, pow_GR_17ones
1049 }
1050 { .mfb
1051 andcm pow_GR_sign_Y = pow_GR_signexp_Y, pow_GR_17ones
1052 (p6) fma.s1 POW_U = POW_NORM_Y,POW_r1,f0
1053 (p12) br.cond.spnt POW_Y_0 // Branch if y=zero, x not zero or nan
1054 }
1055 ;;
1057 { .mfi
1058 ldfe POW_log2_by_128_lo = [pow_AD_P], 16
1059 (p7) fma.s1 POW_Z2 = POW_twoV, POW_U, f0
1060 nop.i 999
1061 }
1062 { .mfi
1063 ldfe POW_log2_by_128_hi = [pow_AD_Q], 16
1064 nop.f 999
1065 nop.i 999
1066 }
1067 ;;
1069 { .mfi
1070 nop.m 999
1071 fcvt.xf POW_float_int_Y = POW_int_Y
1072 nop.i 999
1073 }
1074 { .mfi
1075 nop.m 999
1076 (p7) fma.s1 POW_G = POW_K, POW_log2_hi, POW_T
1077 adds pow_AD_tbl1 = pow_tbl1 - pow_Tt, pow_AD_Q
1078 }
1079 ;;
1081 // p11 = TRUE ==> X is NEGATIVE but not inf
1082 { .mfi
1083 nop.m 999
1084 fclass.m p11,p0 = POW_NORM_X, 0x1a
1085 nop.i 999
1086 }
1087 { .mfi
1088 nop.m 999
1089 (p7) fma.s1 POW_delta = POW_K, POW_log2_lo, f0
1090 adds pow_AD_tbl2 = pow_tbl2 - pow_tbl1, pow_AD_tbl1
1091 }
1092 ;;
1094 { .mfi
1095 nop.m 999
1096 (p6) fma.s1 POW_Z = POW_twoV, POW_U, f0
1097 nop.i 999
1098 }
1099 { .mfi
1100 nop.m 999
1101 fma.s1 POW_v2 = POW_P1, POW_r, POW_P0
1102 nop.i 999
1103 }
1104 ;;
1106 // p11 = TRUE ==> X is NEGATIVE but not inf
1107 // p12 = TRUE ==> X is NEGATIVE AND Y already even int
1108 // p13 = TRUE ==> X is NEGATIVE AND Y possible int
1109 { .mfi
1110 nop.m 999
1111 (p7) fma.s1 POW_Z = POW_NORM_Y, POW_G, POW_Z2
1112 (p11) cmp.gt.unc p12,p13 = pow_GR_exp_Y, pow_GR_10033
1113 }
1114 { .mfi
1115 nop.m 999
1116 fma.s1 POW_Gpr = POW_G, f1, POW_r
1117 nop.i 999
1118 }
1119 ;;
1121 { .mfi
1122 nop.m 999
1123 fma.s1 POW_Yrcub = POW_rsq, POW_U, f0
1124 nop.i 999
1125 }
1126 { .mfi
1127 nop.m 999
1128 fma.s1 POW_p = POW_rsq, POW_P2, POW_v2
1129 nop.i 999
1130 }
1131 ;;
1133 // Test if x inf
1134 { .mfi
1135 nop.m 999
1136 fclass.m p15,p0 = POW_NORM_X, 0x23
1137 nop.i 999
1138 }
1139 // By adding RSHF (1.1000...*2^63) we put integer part in rightmost significand
1140 { .mfi
1141 nop.m 999
1142 fma.s1 POW_W1 = POW_Z, POW_inv_log2_by_128, POW_RSHF
1143 nop.i 999
1144 }
1145 ;;
1147 // p13 = TRUE ==> X is NEGATIVE AND Y possible int
1148 // p10 = TRUE ==> X is NEG and Y is an int
1149 // p12 = TRUE ==> X is NEG and Y is not an int
1150 { .mfi
1151 nop.m 999
1152 (p13) fcmp.eq.unc.s1 p10,p12 = POW_float_int_Y, POW_NORM_Y
1153 mov pow_GR_xneg_yodd = 0
1154 }
1155 { .mfi
1156 nop.m 999
1157 fma.s1 POW_Y_Gpr = POW_NORM_Y, POW_Gpr, f0
1158 nop.i 999
1159 }
1160 ;;
1162 // p11 = TRUE ==> X is +1.0
1163 { .mfi
1164 nop.m 999
1165 fcmp.eq.s1 p11,p0 = POW_NORM_X, f1
1166 nop.i 999
1167 }
1168 ;;
1170 // Extract rounded integer from rightmost significand of POW_W1
1171 // By subtracting RSHF we get rounded integer POW_Nfloat
1172 { .mfi
1173 getf.sig pow_GR_int_N = POW_W1
1174 fms.s1 POW_Nfloat = POW_W1, f1, POW_RSHF
1175 nop.i 999
1176 }
1177 { .mfb
1178 nop.m 999
1179 fma.s1 POW_Z3 = POW_p, POW_Yrcub, f0
1180 (p12) br.cond.spnt POW_X_NEG_Y_NONINT // Branch if x neg, y not integer
1181 }
1182 ;;
1184 // p7 = TRUE ==> Y is +1.0
1185 // p12 = TRUE ==> X is NEGATIVE AND Y is an odd integer
1186 { .mfi
1187 getf.exp pow_GR_signexp_Y_Gpr = POW_Y_Gpr
1188 fcmp.eq.s1 p7,p0 = POW_NORM_Y, f1 // Test for y=1.0
1189 (p10) tbit.nz.unc p12,p0 = pow_GR_sig_int_Y,0
1190 }
1191 { .mfb
1192 nop.m 999
1193 (p11) fma.s.s0 f8 = f1,f1,f0 // If x=1, result is +1
1194 (p15) br.cond.spnt POW_X_INF
1195 }
1196 ;;
1198 // Test x and y and flag denormal
1199 { .mfi
1200 nop.m 999
1201 fcmp.eq.s0 p15,p0 = f8,f9
1202 nop.i 999
1203 }
1204 { .mfb
1205 nop.m 999
1206 fma.s1 POW_e3 = POW_NORM_Y, POW_delta, f0
1207 (p11) br.ret.spnt b0 // Early exit if x=1.0, result is +1
1208 }
1209 ;;
1211 { .mfi
1212 (p12) mov pow_GR_xneg_yodd = 1
1213 fnma.s1 POW_f12 = POW_Nfloat, POW_log2_by_128_lo, f1
1214 nop.i 999
1215 }
1216 { .mfb
1217 nop.m 999
1218 fnma.s1 POW_s = POW_Nfloat, POW_log2_by_128_hi, POW_Z
1219 (p7) br.ret.spnt b0 // Early exit if y=1.0, result is x
1220 }
1221 ;;
1223 { .mmi
1224 and pow_GR_index1 = 0x0f, pow_GR_int_N
1225 and pow_GR_index2 = 0x70, pow_GR_int_N
1226 shr pow_int_GR_M = pow_GR_int_N, 7 // M = N/128
1227 }
1228 ;;
1230 { .mfi
1231 shladd pow_AD_T1 = pow_GR_index1, 4, pow_AD_tbl1
1232 fma.s1 POW_q = POW_Z3, POW_Q1, POW_Q0_half
1233 add pow_int_GR_M = pow_GR_16ones, pow_int_GR_M
1234 }
1235 { .mfi
1236 add pow_AD_T2 = pow_AD_tbl2, pow_GR_index2
1237 fma.s1 POW_Z3sq = POW_Z3, POW_Z3, f0
1238 nop.i 999
1239 }
1240 ;;
1242 { .mmi
1243 ldfe POW_T1 = [pow_AD_T1]
1244 ldfe POW_T2 = [pow_AD_T2]
1245 nop.i 999
1246 }
1247 ;;
1249 // f123 = f12*(e3+1) = f12*e3+f12
1250 { .mfi
1251 setf.exp POW_2M = pow_int_GR_M
1252 fma.s1 POW_f123 = POW_e3,POW_f12,POW_f12
1253 nop.i 999
1254 }
1255 { .mfi
1256 nop.m 999
1257 fma.s1 POW_ssq = POW_s, POW_s, f0
1258 nop.i 999
1259 }
1260 ;;
1262 { .mfi
1263 nop.m 999
1264 fma.s1 POW_v2 = POW_s, POW_Q2, POW_Q1
1265 and pow_GR_exp_Y_Gpr = pow_GR_signexp_Y_Gpr, pow_GR_17ones
1266 }
1267 ;;
1269 { .mfi
1270 cmp.ne p12,p13 = pow_GR_xneg_yodd, r0
1271 fma.s1 POW_q = POW_Z3sq, POW_q, POW_Z3
1272 sub pow_GR_true_exp_Y_Gpr = pow_GR_exp_Y_Gpr, pow_GR_16ones
1273 }
1274 ;;
1276 // p8 TRUE ==> |Y(G + r)| >= 7
1278 // single
1279 // -2^7 -2^6 2^6 2^7
1280 // -----+-----+----+ ... +-----+-----+-----
1281 // p8 | p9 | p8
1282 // | | p10 | |
1284 // Form signexp of constants to indicate overflow
1285 { .mfi
1286 mov pow_GR_big_pos = 0x1007f
1287 nop.f 999
1288 cmp.le p8,p9 = 7, pow_GR_true_exp_Y_Gpr
1289 }
1290 { .mfi
1291 mov pow_GR_big_neg = 0x3007f
1292 nop.f 999
1293 andcm pow_GR_sign_Y_Gpr = pow_GR_signexp_Y_Gpr, pow_GR_17ones
1294 }
1295 ;;
1297 // Form big positive and negative constants to test for possible overflow
1298 // Scale both terms of the polynomial by POW_f123
1299 { .mfi
1300 setf.exp POW_big_pos = pow_GR_big_pos
1301 fma.s1 POW_ssq = POW_ssq, POW_f123, f0
1302 (p9) cmp.le.unc p0,p10 = 6, pow_GR_true_exp_Y_Gpr
1303 }
1304 { .mfb
1305 setf.exp POW_big_neg = pow_GR_big_neg
1306 fma.s1 POW_1ps = POW_s, POW_f123, POW_f123
1307 (p8) br.cond.spnt POW_OVER_UNDER_X_NOT_INF
1308 }
1309 ;;
1311 { .mfi
1312 nop.m 999
1313 (p12) fnma.s1 POW_T1T2 = POW_T1, POW_T2, f0
1314 nop.i 999
1315 }
1316 { .mfi
1317 nop.m 999
1318 (p13) fma.s1 POW_T1T2 = POW_T1, POW_T2, f0
1319 nop.i 999
1320 }
1321 ;;
1323 { .mfi
1324 nop.m 999
1325 fma.s1 POW_v210 = POW_s, POW_v2, POW_Q0_half
1326 nop.i 999
1327 }
1328 { .mfi
1329 nop.m 999
1330 fma.s1 POW_2Mqp1 = POW_2M, POW_q, POW_2M
1331 nop.i 999
1332 }
1333 ;;
1335 { .mfi
1336 nop.m 999
1337 fma.s1 POW_es = POW_ssq, POW_v210, POW_1ps
1338 nop.i 999
1339 }
1340 { .mfi
1341 nop.m 999
1342 fma.s1 POW_A = POW_T1T2, POW_2Mqp1, f0
1343 nop.i 999
1344 }
1345 ;;
1347 // Dummy op to set inexact
1348 { .mfi
1349 nop.m 999
1350 fma.s0 POW_tmp = POW_2M, POW_q, POW_2M
1351 nop.i 999
1352 }
1353 ;;
1355 { .mfb
1356 nop.m 999
1357 fma.s.s0 f8 = POW_A, POW_es, f0
1358 (p10) br.ret.sptk b0 // Exit main branch if no over/underflow
1359 }
1360 ;;
1362 // POSSIBLE_OVER_UNDER
1363 // p6 = TRUE ==> Y_Gpr negative
1364 // Result is already computed. We just need to know if over/underflow occurred.
1366 { .mfb
1367 cmp.eq p0,p6 = pow_GR_sign_Y_Gpr, r0
1368 nop.f 999
1369 (p6) br.cond.spnt POW_POSSIBLE_UNDER
1370 }
1371 ;;
1373 // POSSIBLE_OVER
1374 // We got an answer.
1375 // overflow is a possibility, not a certainty
1378 // We define an overflow when the answer with
1379 // WRE set
1380 // user-defined rounding mode
1382 // double
1383 // Largest double is 7FE (biased double)
1384 // 7FE - 3FF + FFFF = 103FE
1385 // Create + largest_double_plus_ulp
1386 // Create - largest_double_plus_ulp
1387 // Calculate answer with WRE set.
1389 // single
1390 // Largest single is FE (biased double)
1391 // FE - 7F + FFFF = 1007E
1392 // Create + largest_single_plus_ulp
1393 // Create - largest_single_plus_ulp
1394 // Calculate answer with WRE set.
1396 // Cases when answer is ldn+1 are as follows:
1397 // ldn ldn+1
1398 // --+----------|----------+------------
1399 // |
1400 // +inf +inf -inf
1401 // RN RN
1402 // RZ
1404 // Put in s2 (td set, wre set)
1405 { .mfi
1406 nop.m 999
1407 fsetc.s2 0x7F,0x42
1408 nop.i 999
1409 }
1410 ;;
1412 { .mfi
1413 nop.m 999
1414 fma.s.s2 POW_wre_urm_f8 = POW_A, POW_es, f0
1415 nop.i 999
1416 }
1417 ;;
1419 // Return s2 to default
1420 { .mfi
1421 nop.m 999
1422 fsetc.s2 0x7F,0x40
1423 nop.i 999
1424 }
1425 ;;
1427 // p7 = TRUE ==> yes, we have an overflow
1428 { .mfi
1429 nop.m 999
1430 fcmp.ge.s1 p7, p8 = POW_wre_urm_f8, POW_big_pos
1431 nop.i 999
1432 }
1433 ;;
1435 { .mfi
1436 nop.m 999
1437 (p8) fcmp.le.s1 p7, p0 = POW_wre_urm_f8, POW_big_neg
1438 nop.i 999
1439 }
1440 ;;
1442 { .mbb
1443 (p7) mov pow_GR_tag = 30
1444 (p7) br.cond.spnt __libm_error_region // Branch if overflow
1445 br.ret.sptk b0 // Exit if did not overflow
1446 }
1447 ;;
1450 POW_POSSIBLE_UNDER:
1451 // We got an answer. input was < -2^9 but > -2^10 (double)
1452 // We got an answer. input was < -2^6 but > -2^7 (float)
1453 // underflow is a possibility, not a certainty
1455 // We define an underflow when the answer with
1456 // ftz set
1457 // is zero (tiny numbers become zero)
1458 // Notice (from below) that if we have an unlimited exponent range,
1459 // then there is an extra machine number E between the largest denormal and
1460 // the smallest normal.
1461 // So if with unbounded exponent we round to E or below, then we are
1462 // tiny and underflow has occurred.
1463 // But notice that you can be in a situation where we are tiny, namely
1464 // rounded to E, but when the exponent is bounded we round to smallest
1465 // normal. So the answer can be the smallest normal with underflow.
1466 // E
1467 // -----+--------------------+--------------------+-----
1468 // | | |
1469 // 1.1...10 2^-3fff 1.1...11 2^-3fff 1.0...00 2^-3ffe
1470 // 0.1...11 2^-3ffe (biased, 1)
1471 // largest dn smallest normal
1473 // Form small constant (2^-170) to correct underflow result near region of
1474 // smallest denormal in round-nearest.
1476 // Put in s2 (td set, ftz set)
1477 .pred.rel "mutex",p12,p13
1478 { .mfi
1479 mov pow_GR_Fpsr = ar40 // Read the fpsr--need to check rc.s0
1480 fsetc.s2 0x7F,0x41
1481 mov pow_GR_rcs0_mask = 0x0c00 // Set mask for rc.s0
1482 }
1483 { .mfi
1484 (p12) mov pow_GR_tmp = 0x2ffff - 170
1485 nop.f 999
1486 (p13) mov pow_GR_tmp = 0x0ffff - 170
1487 }
1488 ;;
1490 { .mfi
1491 setf.exp POW_eps = pow_GR_tmp // Form 2^-170
1492 fma.s.s2 POW_ftz_urm_f8 = POW_A, POW_es, f0
1493 nop.i 999
1494 }
1495 ;;
1497 // Return s2 to default
1498 { .mfi
1499 nop.m 999
1500 fsetc.s2 0x7F,0x40
1501 nop.i 999
1502 }
1503 ;;
1505 // p7 = TRUE ==> yes, we have an underflow
1506 { .mfi
1507 nop.m 999
1508 fcmp.eq.s1 p7, p0 = POW_ftz_urm_f8, f0
1509 nop.i 999
1510 }
1511 ;;
1513 { .mmi
1514 (p7) and pow_GR_rcs0 = pow_GR_rcs0_mask, pow_GR_Fpsr // Isolate rc.s0
1515 ;;
1516 (p7) cmp.eq.unc p6,p0 = pow_GR_rcs0, r0 // Test for round to nearest
1517 nop.i 999
1518 }
1519 ;;
1521 // Tweak result slightly if underflow to get correct rounding near smallest
1522 // denormal if round-nearest
1523 { .mfi
1524 nop.m 999
1525 (p6) fms.s.s0 f8 = POW_A, POW_es, POW_eps
1526 nop.i 999
1527 }
1528 { .mbb
1529 (p7) mov pow_GR_tag = 31
1530 (p7) br.cond.spnt __libm_error_region // Branch if underflow
1531 br.ret.sptk b0 // Exit if did not underflow
1532 }
1533 ;;
1535 POW_X_DENORM:
1536 // Here if x unorm. Use the NORM_X for getf instructions, and then back
1537 // to normal path
1538 { .mfi
1539 getf.exp pow_GR_signexp_X = POW_NORM_X
1540 nop.f 999
1541 nop.i 999
1542 }
1543 ;;
1545 { .mib
1546 getf.sig pow_GR_sig_X = POW_NORM_X
1547 nop.i 999
1548 br.cond.sptk POW_COMMON
1549 }
1550 ;;
1552 POW_X_0:
1553 // Here if x=0 and y not nan
1554 //
1555 // We have the following cases:
1556 // p6 x=0 and y>0 and is an integer (may be even or odd)
1557 // p7 x=0 and y>0 and is NOT an integer, return +0
1558 // p8 x=0 and y>0 and so big as to always be an even integer, return +0
1559 // p9 x=0 and y>0 and may not be integer
1560 // p10 x=0 and y>0 and is an odd integer, return x
1561 // p11 x=0 and y>0 and is an even integer, return +0
1562 // p12 used in dummy fcmp to set denormal flag if y=unorm
1563 // p13 x=0 and y>0
1564 // p14 x=0 and y=0, branch to code for calling error handling
1565 // p15 x=0 and y<0, branch to code for calling error handling
1566 //
1567 { .mfi
1568 getf.sig pow_GR_sig_int_Y = POW_int_Y // Get signif of int_Y
1569 fcmp.lt.s1 p15,p13 = f9, f0 // Test for y<0
1570 and pow_GR_exp_Y = pow_GR_signexp_Y, pow_GR_17ones
1571 }
1572 { .mfb
1573 cmp.ne p14,p0 = pow_GR_y_zero,r0 // Test for y=0
1574 fcvt.xf POW_float_int_Y = POW_int_Y
1575 (p14) br.cond.spnt POW_X_0_Y_0 // Branch if x=0 and y=0
1576 }
1577 ;;
1579 // If x=0 and y>0, test y and flag denormal
1580 { .mfb
1581 (p13) cmp.gt.unc p8,p9 = pow_GR_exp_Y, pow_GR_10033 // Test y +big = even int
1582 (p13) fcmp.eq.s0 p12,p0 = f9,f0 // If x=0, y>0 dummy op to flag denormal
1583 (p15) br.cond.spnt POW_X_0_Y_NEG // Branch if x=0 and y<0
1584 }
1585 ;;
1587 // Here if x=0 and y>0
1588 { .mfi
1589 nop.m 999
1590 (p9) fcmp.eq.unc.s1 p6,p7 = POW_float_int_Y, POW_NORM_Y // Test y=int
1591 nop.i 999
1592 }
1593 { .mfi
1594 nop.m 999
1595 (p8) fma.s.s0 f8 = f0,f0,f0 // If x=0, y>0 and large even int, return +0
1596 nop.i 999
1597 }
1598 ;;
1600 { .mfi
1601 nop.m 999
1602 (p7) fma.s.s0 f8 = f0,f0,f0 // Result +0 if x=0 and y>0 and not integer
1603 (p6) tbit.nz.unc p10,p11 = pow_GR_sig_int_Y,0 // If y>0 int, test y even/odd
1604 }
1605 ;;
1607 // Note if x=0, y>0 and odd integer, just return x
1608 { .mfb
1609 nop.m 999
1610 (p11) fma.s.s0 f8 = f0,f0,f0 // Result +0 if x=0 and y even integer
1611 br.ret.sptk b0 // Exit if x=0 and y>0
1612 }
1613 ;;
1615 POW_X_0_Y_0:
1616 // When X is +-0 and Y is +-0, IEEE returns 1.0
1617 // We call error support with this value
1619 { .mfb
1620 mov pow_GR_tag = 32
1621 fma.s.s0 f8 = f1,f1,f0
1622 br.cond.sptk __libm_error_region
1623 }
1624 ;;
1626 POW_X_0_Y_NEG:
1627 // When X is +-0 and Y is negative, IEEE returns
1628 // X Y answer
1629 // +0 -odd int +inf
1630 // -0 -odd int -inf
1632 // +0 !-odd int +inf
1633 // -0 !-odd int +inf
1635 // p6 == Y is a floating point number outside the integer.
1636 // Hence it is an integer and is even.
1637 // return +inf
1639 // p7 == Y is a floating point number within the integer range.
1640 // p9 == (int_Y = NORM_Y), Y is an integer, which may be odd or even.
1641 // p11 odd
1642 // return (sign_of_x)inf
1643 // p12 even
1644 // return +inf
1645 // p10 == Y is not an integer
1646 // return +inf
1647 //
1649 { .mfi
1650 nop.m 999
1651 nop.f 999
1652 cmp.gt p6,p7 = pow_GR_exp_Y, pow_GR_10033
1653 }
1654 ;;
1656 { .mfi
1657 mov pow_GR_tag = 33
1658 (p7) fcmp.eq.unc.s1 p9,p10 = POW_float_int_Y, POW_NORM_Y
1659 nop.i 999
1660 }
1661 ;;
1663 { .mfb
1664 nop.m 999
1665 (p6) frcpa.s0 f8,p13 = f1, f0
1666 (p6) br.cond.sptk __libm_error_region // x=0, y<0, y large neg int
1667 }
1668 ;;
1670 { .mfb
1671 nop.m 999
1672 (p10) frcpa.s0 f8,p13 = f1, f0
1673 (p10) br.cond.sptk __libm_error_region // x=0, y<0, y not int
1674 }
1675 ;;
1677 // x=0, y<0, y an int
1678 { .mib
1679 nop.m 999
1680 (p9) tbit.nz.unc p11,p12 = pow_GR_sig_int_Y,0
1681 nop.b 999
1682 }
1683 ;;
1685 { .mfi
1686 nop.m 999
1687 (p12) frcpa.s0 f8,p13 = f1,f0
1688 nop.i 999
1689 }
1690 ;;
1692 { .mfb
1693 nop.m 999
1694 (p11) frcpa.s0 f8,p13 = f1,f8
1695 br.cond.sptk __libm_error_region
1696 }
1697 ;;
1700 POW_Y_0:
1701 // Here for y zero, x anything but zero and nan
1702 // Set flag if x denormal
1703 // Result is +1.0
1704 { .mfi
1705 nop.m 999
1706 fcmp.eq.s0 p6,p0 = f8,f0 // Sets flag if x denormal
1707 nop.i 999
1708 }
1709 { .mfb
1710 nop.m 999
1711 fma.s.s0 f8 = f1,f1,f0
1712 br.ret.sptk b0
1713 }
1714 ;;
1717 POW_X_INF:
1718 // Here when X is +-inf
1720 // X +inf Y +inf +inf
1721 // X -inf Y +inf +inf
1723 // X +inf Y >0 +inf
1724 // X -inf Y >0, !odd integer +inf <== (-inf)^0.5 = +inf !!
1725 // X -inf Y >0, odd integer -inf
1727 // X +inf Y -inf +0
1728 // X -inf Y -inf +0
1730 // X +inf Y <0 +0
1731 // X -inf Y <0, !odd integer +0
1732 // X -inf Y <0, odd integer -0
1734 // X + inf Y=+0 +1
1735 // X + inf Y=-0 +1
1736 // X - inf Y=+0 +1
1737 // X - inf Y=-0 +1
1739 // p13 == Y negative
1740 // p14 == Y positive
1742 // p6 == Y is a floating point number outside the integer.
1743 // Hence it is an integer and is even.
1744 // p13 == (Y negative)
1745 // return +inf
1746 // p14 == (Y positive)
1747 // return +0
1749 // p7 == Y is a floating point number within the integer range.
1750 // p9 == (int_Y = NORM_Y), Y is an integer, which may be odd or even.
1751 // p11 odd
1752 // p13 == (Y negative)
1753 // return (sign_of_x)inf
1754 // p14 == (Y positive)
1755 // return (sign_of_x)0
1756 // pxx even
1757 // p13 == (Y negative)
1758 // return +inf
1759 // p14 == (Y positive)
1760 // return +0
1762 // pxx == Y is not an integer
1763 // p13 == (Y negative)
1764 // return +inf
1765 // p14 == (Y positive)
1766 // return +0
1767 //
1769 // If x=inf, test y and flag denormal
1770 { .mfi
1771 nop.m 999
1772 fcmp.eq.s0 p10,p11 = f9,f0
1773 nop.i 999
1774 }
1775 ;;
1777 { .mfi
1778 nop.m 999
1779 fcmp.lt.s0 p13,p14 = POW_NORM_Y,f0
1780 cmp.gt p6,p7 = pow_GR_exp_Y, pow_GR_10033
1781 }
1782 { .mfi
1783 nop.m 999
1784 fclass.m p12,p0 = f9, 0x23 //@inf
1785 nop.i 999
1786 }
1787 ;;
1789 { .mfi
1790 nop.m 999
1791 fclass.m p15,p0 = f9, 0x07 //@zero
1792 nop.i 999
1793 }
1794 ;;
1796 { .mfb
1797 nop.m 999
1798 (p15) fmerge.s f8 = f1,f1 // Return +1.0 if x=inf, y=0
1799 (p15) br.ret.spnt b0 // Exit if x=inf, y=0
1800 }
1801 ;;
1803 { .mfi
1804 nop.m 999
1805 (p14) frcpa.s1 f8,p10 = f1,f0 // If x=inf, y>0, assume result +inf
1806 nop.i 999
1807 }
1808 { .mfb
1809 nop.m 999
1810 (p13) fma.s.s0 f8 = f0,f0,f0 // If x=inf, y<0, assume result +0.0
1811 (p12) br.ret.spnt b0 // Exit if x=inf, y=inf
1812 }
1813 ;;
1815 // Here if x=inf, and 0 < |y| < inf. Need to correct results if y odd integer.
1816 { .mfi
1817 nop.m 999
1818 (p7) fcmp.eq.unc.s1 p9,p0 = POW_float_int_Y, POW_NORM_Y // Is y integer?
1819 nop.i 999
1820 }
1821 ;;
1823 { .mfi
1824 nop.m 999
1825 nop.f 999
1826 (p9) tbit.nz.unc p11,p0 = pow_GR_sig_int_Y,0 // Test for y odd integer
1827 }
1828 ;;
1830 { .mfb
1831 nop.m 999
1832 (p11) fmerge.s f8 = POW_NORM_X,f8 // If y odd integer use sign of x
1833 br.ret.sptk b0 // Exit for x=inf, 0 < |y| < inf
1834 }
1835 ;;
1838 POW_X_NEG_Y_NONINT:
1839 // When X is negative and Y is a non-integer, IEEE
1840 // returns a qnan indefinite.
1841 // We call error support with this value
1843 { .mfb
1844 mov pow_GR_tag = 34
1845 frcpa.s0 f8,p6 = f0,f0
1846 br.cond.sptk __libm_error_region
1847 }
1848 ;;
1850 POW_X_NAN:
1851 // Here if x=nan, y not nan
1852 { .mfi
1853 nop.m 999
1854 fclass.m p9,p13 = f9, 0x07 // Test y=zero
1855 nop.i 999
1856 }
1857 ;;
1859 { .mfb
1860 nop.m 999
1861 (p13) fma.s.s0 f8 = f8,f1,f0
1862 (p13) br.ret.sptk b0 // Exit if x nan, y anything but zero or nan
1863 }
1864 ;;
1866 POW_X_NAN_Y_0:
1867 // When X is a NAN and Y is zero, IEEE returns 1.
1868 // We call error support with this value.
1869 { .mfi
1870 nop.m 999
1871 fcmp.eq.s0 p6,p0 = f8,f0 // Dummy op to set invalid on snan
1872 nop.i 999
1873 }
1874 { .mfb
1875 mov pow_GR_tag = 35
1876 fma.s.s0 f8 = f0,f0,f1
1877 br.cond.sptk __libm_error_region
1878 }
1879 ;;
1882 POW_OVER_UNDER_X_NOT_INF:
1884 // p8 is TRUE for overflow
1885 // p9 is TRUE for underflow
1887 // if y is infinity, we should not over/underflow
1889 { .mfi
1890 nop.m 999
1891 fcmp.eq.s1 p14, p13 = POW_xsq,f1 // Test |x|=1
1892 cmp.eq p8,p9 = pow_GR_sign_Y_Gpr, r0
1893 }
1894 ;;
1896 { .mfi
1897 nop.m 999
1898 (p14) fclass.m.unc p15, p0 = f9, 0x23 // If |x|=1, test y=inf
1899 nop.i 999
1900 }
1901 { .mfi
1902 nop.m 999
1903 (p13) fclass.m.unc p11,p0 = f9, 0x23 // If |x| not 1, test y=inf
1904 nop.i 999
1905 }
1906 ;;
1908 // p15 = TRUE if |x|=1, y=inf, return +1
1909 { .mfb
1910 nop.m 999
1911 (p15) fma.s.s0 f8 = f1,f1,f0 // If |x|=1, y=inf, result +1
1912 (p15) br.ret.spnt b0 // Exit if |x|=1, y=inf
1913 }
1914 ;;
1916 .pred.rel "mutex",p8,p9
1917 { .mfb
1918 (p8) setf.exp f8 = pow_GR_17ones // If exp(+big), result inf
1919 (p9) fmerge.s f8 = f0,f0 // If exp(-big), result 0
1920 (p11) br.ret.sptk b0 // Exit if |x| not 1, y=inf
1921 }
1922 ;;
1924 { .mfb
1925 nop.m 999
1926 nop.f 999
1927 br.cond.sptk POW_OVER_UNDER_ERROR // Branch if y not inf
1928 }
1929 ;;
1932 POW_Y_NAN:
1933 // Here if y=nan, x anything
1934 // If x = +1 then result is +1, else result is quiet Y
1935 { .mfi
1936 nop.m 999
1937 fcmp.eq.s1 p10,p9 = POW_NORM_X, f1
1938 nop.i 999
1939 }
1940 ;;
1942 { .mfi
1943 nop.m 999
1944 (p10) fcmp.eq.s0 p6,p0 = f9,f1 // Set invalid, even if x=+1
1945 nop.i 999
1946 }
1947 ;;
1949 { .mfi
1950 nop.m 999
1951 (p10) fma.s.s0 f8 = f1,f1,f0
1952 nop.i 999
1953 }
1954 { .mfb
1955 nop.m 999
1956 (p9) fma.s.s0 f8 = f9,f8,f0
1957 br.ret.sptk b0 // Exit y=nan
1958 }
1959 ;;
1962 POW_OVER_UNDER_ERROR:
1963 // Here if we have overflow or underflow.
1964 // Enter with p12 true if x negative and y odd int to force -0 or -inf
1966 { .mfi
1967 sub pow_GR_17ones_m1 = pow_GR_17ones, r0, 1
1968 nop.f 999
1969 mov pow_GR_one = 0x1
1970 }
1971 ;;
1973 // overflow, force inf with O flag
1974 { .mmb
1975 (p8) mov pow_GR_tag = 30
1976 (p8) setf.exp POW_tmp = pow_GR_17ones_m1
1977 nop.b 999
1978 }
1979 ;;
1981 // underflow, force zero with I, U flags
1982 { .mmi
1983 (p9) mov pow_GR_tag = 31
1984 (p9) setf.exp POW_tmp = pow_GR_one
1985 nop.i 999
1986 }
1987 ;;
1989 { .mfi
1990 nop.m 999
1991 fma.s.s0 f8 = POW_tmp, POW_tmp, f0
1992 nop.i 999
1993 }
1994 ;;
1996 // p12 x is negative and y is an odd integer, change sign of result
1997 { .mfi
1998 nop.m 999
1999 (p12) fnma.s.s0 f8 = POW_tmp, POW_tmp, f0
2000 nop.i 999
2001 }
2002 ;;
2004 WEAK_LIBM_END(powf)
2005 libm_alias_float_other (__pow, pow)
2006 #ifdef SHARED
2007 .symver powf,powf@@GLIBC_2.27
2008 .weak __powf_compat
2009 .set __powf_compat,__powf
2010 .symver __powf_compat,powf@GLIBC_2.2
2011 #endif
2014 LOCAL_LIBM_ENTRY(__libm_error_region)
2016 .prologue
2017 { .mfi
2018 add GR_Parameter_Y=-32,sp // Parameter 2 value
2019 nop.f 0
2020 .save ar.pfs,GR_SAVE_PFS
2021 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
2022 }
2023 { .mfi
2024 .fframe 64
2025 add sp=-64,sp // Create new stack
2026 nop.f 0
2027 mov GR_SAVE_GP=gp // Save gp
2028 };;
2030 { .mmi
2031 stfs [GR_Parameter_Y] = POW_NORM_Y,16 // STORE Parameter 2 on stack
2032 add GR_Parameter_X = 16,sp // Parameter 1 address
2033 .save b0, GR_SAVE_B0
2034 mov GR_SAVE_B0=b0 // Save b0
2035 };;
2037 .body
2038 { .mib
2039 stfs [GR_Parameter_X] = POW_NORM_X // STORE Parameter 1 on stack
2040 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address
2041 nop.b 0
2042 }
2043 { .mib
2044 stfs [GR_Parameter_Y] = f8 // STORE Parameter 3 on stack
2045 add GR_Parameter_Y = -16,GR_Parameter_Y
2046 br.call.sptk b0=__libm_error_support# // Call error handling function
2047 };;
2049 { .mmi
2050 add GR_Parameter_RESULT = 48,sp
2051 nop.m 0
2052 nop.i 0
2053 };;
2055 { .mmi
2056 ldfs f8 = [GR_Parameter_RESULT] // Get return result off stack
2057 .restore sp
2058 add sp = 64,sp // Restore stack pointer
2059 mov b0 = GR_SAVE_B0 // Restore return address
2060 };;
2062 { .mib
2063 mov gp = GR_SAVE_GP // Restore gp
2064 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
2065 br.ret.sptk b0 // Return
2066 };;
2068 LOCAL_LIBM_END(__libm_error_region)
2070 .type __libm_error_support#,@function
2071 .global __libm_error_support#