1.file "pow.s" 2 3 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. 22 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^52 <= |y| < 2^53, y odd integer. 62// 12/20/01 Fixed monotonity problem in round to nearest. 63// 02/08/02 Fixed overflow/underflow cases that were not calling error support. 64// 05/20/02 Cleaned up namespace and sf0 syntax 65// 08/29/02 Improved Itanium 2 performance 66// 09/21/02 Added branch for |y*log(x)|<2^-11 to fix monotonicity problems. 67// 02/10/03 Reordered header: .section, .global, .proc, .align 68// 03/31/05 Reformatted delimiters between data tables 69// 70// API 71//============================================================== 72// double pow(double x, double 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 tabluate as two doubles, T and t, where T +t is the value itself. 112// 113// Log(x) = (K Log(2)_hi + T) + (Log(2)_hi + t) + 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// 120// Then, 121// 122// yLog(x) = yG + y delta + y(r-rsq/2) + yp 123// yLog(x) = Z1 + e3 + Z2 + Z3 + (e2 + e3) 124// 125// 126// exp(yLog(x)) = exp(Z1 + Z2 + Z3) exp(e1 + e2 + e3) 127// 128// 129// exp(Z3) is another series. 130// exp(e1 + e2 + e3) is approximated as f3 = 1 + (e1 + e2 + e3) 131// 132// Z1 (128/log2) = number of log2/128 in Z1 is N1 133// Z2 (128/log2) = number of log2/128 in Z2 is N2 134// 135// s1 = Z1 - N1 log2/128 136// s2 = Z2 - N2 log2/128 137// 138// s = s1 + s2 139// N = N1 + N2 140// 141// exp(Z1 + Z2) = exp(Z) 142// exp(Z) = exp(s) exp(N log2/128) 143// 144// exp(r) = exp(Z - N log2/128) 145// 146// r = s + d = (Z - N (log2/128)_hi) -N (log2/128)_lo 147// = Z - N (log2/128) 148// 149// Z = s+d +N (log2/128) 150// 151// exp(Z) = exp(s) (1+d) exp(N log2/128) 152// 153// N = M 128 + n 154// 155// N log2/128 = M log2 + n log2/128 156// 157// n is 8 binary digits = n_7n_6...n_1 158// 159// n log2/128 = n_7n_6n_5 16 log2/128 + n_4n_3n_2n_1 log2/128 160// n log2/128 = n_7n_6n_5 log2/8 + n_4n_3n_2n_1 log2/128 161// n log2/128 = I2 log2/8 + I1 log2/128 162// 163// N log2/128 = M log2 + I2 log2/8 + I1 log2/128 164// 165// exp(Z) = exp(s) (1+d) exp(log(2^M) + log(2^I2/8) + log(2^I1/128)) 166// exp(Z) = exp(s) (1+d1) (1+d2)(2^M) 2^I2/8 2^I1/128 167// exp(Z) = exp(s) f1 f2 (2^M) 2^I2/8 2^I1/128 168// 169// I1, I2 are table indices. Use a series for exp(s). 170// Then get exp(Z) 171// 172// exp(yLog(x)) = exp(Z1 + Z2 + Z3) exp(e1 + e2 + e3) 173// exp(yLog(x)) = exp(Z) exp(Z3) f3 174// exp(yLog(x)) = exp(Z)f3 exp(Z3) 175// exp(yLog(x)) = A exp(Z3) 176// 177// We actually calculate exp(Z3) -1. 178// Then, 179// exp(yLog(x)) = A + A( exp(Z3) -1) 180// 181 182// Table Generation 183//============================================================== 184 185// The log values 186// ============== 187// The operation (K*log2_hi) must be exact. K is the true exponent of x. 188// If we allow gradual underflow (denormals), K can be represented in 12 bits 189// (as a two's complement number). We assume 13 bits as an engineering 190// precaution. 191// 192// +------------+----------------+-+ 193// | 13 bits | 50 bits | | 194// +------------+----------------+-+ 195// 0 1 66 196// 2 34 197// 198// So we want the lsb(log2_hi) to be 2^-50 199// We get log2 as a quad-extended (15-bit exponent, 128-bit significand) 200// 201// 0 fffe b17217f7d1cf79ab c9e3b39803f2f6af (4...) 202// 203// Consider numbering the bits left to right, starting at 0 thru 127. 204// Bit 0 is the 2^-1 bit; bit 49 is the 2^-50 bit. 205// 206// ...79ab 207// 0111 1001 1010 1011 208// 44 209// 89 210// 211// So if we shift off the rightmost 14 bits, then (shift back only 212// the top half) we get 213// 214// 0 fffe b17217f7d1cf4000 e6af278ece600fcb dabc000000000000 215// 216// Put the right 64-bit signficand in an FR register, convert to double; 217// it is exact. Put the next 128 bits into a quad register and round to double. 218// The true exponent of the low part is -51. 219// 220// hi is 0 fffe b17217f7d1cf4000 221// lo is 0 ffcc e6af278ece601000 222// 223// Convert to double memory format and get 224// 225// hi is 0x3fe62e42fefa39e8 226// lo is 0x3cccd5e4f1d9cc02 227// 228// log2_hi + log2_lo is an accurate value for log2. 229// 230// 231// The T and t values 232// ================== 233// A similar method is used to generate the T and t values. 234// 235// K * log2_hi + T must be exact. 236// 237// Smallest T,t 238// ---------- 239// The smallest T,t is 240// T t 241// 0x3f60040155d58800, 0x3c93bce0ce3ddd81 log(1/frcpa(1+0/256))= +1.95503e-003 242// 243// The exponent is 0x3f6 (biased) or -9 (true). 244// For the smallest T value, what we want is to clip the significand such that 245// when it is shifted right by 9, its lsb is in the bit for 2^-51. The 9 is the 246// specific for the first entry. In general, it is 0xffff - (biased 15-bit 247// exponent). 248 249// Independently, what we have calculated is the table value as a quad 250// precision number. 251// Table entry 1 is 252// 0 fff6 80200aaeac44ef38 338f77605fdf8000 253// 254// We store this quad precision number in a data structure that is 255// sign: 1 256// exponent: 15 257// signficand_hi: 64 (includes explicit bit) 258// signficand_lo: 49 259// Because the explicit bit is included, the significand is 113 bits. 260// 261// Consider significand_hi for table entry 1. 262// 263// 264// +-+--- ... -------+--------------------+ 265// | | 266// +-+--- ... -------+--------------------+ 267// 0 1 4444444455555555556666 268// 2345678901234567890123 269// 270// Labeled as above, bit 0 is 2^0, bit 1 is 2^-1, etc. 271// Bit 42 is 2^-42. If we shift to the right by 9, the bit in 272// bit 42 goes in 51. 273// 274// So what we want to do is shift bits 43 thru 63 into significand_lo. 275// This is shifting bit 42 into bit 63, taking care to retain shifted-off bits. 276// Then shifting (just with signficaand_hi) back into bit 42. 277// 278// The shift_value is 63-42 = 21. In general, this is 279// 63 - (51 -(0xffff - 0xfff6)) 280// For this example, it is 281// 63 - (51 - 9) = 63 - 42 = 21 282// 283// This means we are shifting 21 bits into significand_lo. We must maintain more 284// that a 128-bit signficand not to lose bits. So before the shift we put the 285// 128-bit significand into a 256-bit signficand and then shift. 286// The 256-bit significand has four parts: hh, hl, lh, and ll. 287// 288// Start off with 289// hh hl lh ll 290// <64> <49><15_0> <64_0> <64_0> 291// 292// After shift by 21 (then return for significand_hi), 293// <43><21_0> <21><43> <6><58_0> <64_0> 294// 295// Take the hh part and convert to a double. There is no rounding here. 296// The conversion is exact. The true exponent of the high part is the same as 297// the true exponent of the input quad. 298// 299// We have some 64 plus significand bits for the low part. In this example, we 300// have 70 bits. We want to round this to a double. Put them in a quad and then 301// do a quad fnorm. 302// For this example the true exponent of the low part is 303// true_exponent_of_high - 43 = true_exponent_of_high - (64-21) 304// In general, this is 305// true_exponent_of_high - (64 - shift_value) 306// 307// 308// Largest T,t 309// ---------- 310// The largest T,t is 311// 0x3fe62643fecf9742, 0x3c9e3147684bd37d log(1/frcpa(1+255/256))=+6.92171e-001 312// 313// Table entry 256 is 314// 0 fffe b1321ff67cba178c 51da12f4df5a0000 315// 316// The shift value is 317// 63 - (51 -(0xffff - 0xfffe)) = 13 318// 319// The true exponent of the low part is 320// true_exponent_of_high - (64 - shift_value) 321// -1 - (64-13) = -52 322// Biased as a double, this is 0x3cb 323// 324// 325// 326// So then lsb(T) must be >= 2^-51 327// msb(Klog2_hi) <= 2^12 328// 329// +--------+---------+ 330// | 51 bits | <== largest T 331// +--------+---------+ 332// | 9 bits | 42 bits | <== smallest T 333// +------------+----------------+-+ 334// | 13 bits | 50 bits | | 335// +------------+----------------+-+ 336 337 338// Special Cases 339//============================================================== 340 341// double float 342// overflow error 24 30 343 344// underflow error 25 31 345 346// X zero Y zero 347// +0 +0 +1 error 26 32 348// -0 +0 +1 error 26 32 349// +0 -0 +1 error 26 32 350// -0 -0 +1 error 26 32 351 352// X zero Y negative 353// +0 -odd integer +inf error 27 33 divide-by-zero 354// -0 -odd integer -inf error 27 33 divide-by-zero 355// +0 !-odd integer +inf error 27 33 divide-by-zero 356// -0 !-odd integer +inf error 27 33 divide-by-zero 357// +0 -inf +inf error 27 33 divide-by-zero 358// -0 -inf +inf error 27 33 divide-by-zero 359 360// X zero Y positve 361// +0 +odd integer +0 362// -0 +odd integer -0 363// +0 !+odd integer +0 364// -0 !+odd integer +0 365// +0 +inf +0 366// -0 +inf +0 367// +0 Y NaN quiet Y invalid if Y SNaN 368// -0 Y NaN quiet Y invalid if Y SNaN 369 370// X one 371// -1 Y inf +1 372// -1 Y NaN quiet Y invalid if Y SNaN 373// +1 Y NaN +1 invalid if Y SNaN 374// +1 Y any else +1 375 376// X - Y not integer QNAN error 28 34 invalid 377 378// X NaN Y 0 +1 error 29 35 379// X NaN Y NaN quiet X invalid if X or Y SNaN 380// X NaN Y any else quiet X invalid if X SNaN 381// X !+1 Y NaN quiet Y invalid if Y SNaN 382 383 384// X +inf Y >0 +inf 385// X -inf Y >0, !odd integer +inf 386// X -inf Y >0, odd integer -inf 387 388// X +inf Y <0 +0 389// X -inf Y <0, !odd integer +0 390// X -inf Y <0, odd integer -0 391 392// X +inf Y =0 +1 393// X -inf Y =0 +1 394 395// |X|<1 Y +inf +0 396// |X|<1 Y -inf +inf 397// |X|>1 Y +inf +inf 398// |X|>1 Y -inf +0 399 400// X any Y =0 +1 401 402// Assembly macros 403//============================================================== 404 405// integer registers used 406 407pow_GR_signexp_X = r14 408pow_GR_17ones = r15 409pow_AD_P = r16 410pow_GR_exp_2tom8 = r17 411pow_GR_sig_X = r18 412pow_GR_10033 = r19 413pow_GR_16ones = r20 414 415pow_AD_Tt = r21 416pow_GR_exp_X = r22 417pow_AD_Q = r23 418pow_GR_true_exp_X = r24 419pow_GR_y_zero = r25 420 421pow_GR_exp_Y = r26 422pow_AD_tbl1 = r27 423pow_AD_tbl2 = r28 424pow_GR_offset = r29 425pow_GR_exp_Xm1 = r30 426pow_GR_xneg_yodd = r31 427 428pow_GR_signexp_Xm1 = r35 429pow_GR_int_W1 = r36 430pow_GR_int_W2 = r37 431pow_GR_int_N = r38 432pow_GR_index1 = r39 433pow_GR_index2 = r40 434 435pow_AD_T1 = r41 436pow_AD_T2 = r42 437pow_int_GR_M = r43 438pow_GR_sig_int_Y = r44 439pow_GR_sign_Y_Gpr = r45 440 441pow_GR_17ones_m1 = r46 442pow_GR_one = r47 443pow_GR_sign_Y = r48 444pow_GR_signexp_Y_Gpr = r49 445pow_GR_exp_Y_Gpr = r50 446 447pow_GR_true_exp_Y_Gpr = r51 448pow_GR_signexp_Y = r52 449pow_GR_x_one = r53 450pow_GR_exp_2toM63 = r54 451pow_GR_big_pos = r55 452 453pow_GR_big_neg = r56 454 455GR_SAVE_B0 = r50 456GR_SAVE_GP = r51 457GR_SAVE_PFS = r52 458 459GR_Parameter_X = r53 460GR_Parameter_Y = r54 461GR_Parameter_RESULT = r55 462pow_GR_tag = r56 463 464 465// floating point registers used 466 467POW_B = f32 468POW_NORM_X = f33 469POW_Xm1 = f34 470POW_r1 = f34 471POW_P4 = f35 472 473POW_P5 = f36 474POW_NORM_Y = f37 475POW_Q2 = f38 476POW_Q3 = f39 477POW_P2 = f40 478 479POW_P3 = f41 480POW_P0 = f42 481POW_log2_lo = f43 482POW_r = f44 483POW_Q0_half = f45 484 485POW_Q1 = f46 486POW_tmp = f47 487POW_log2_hi = f48 488POW_Q4 = f49 489POW_P1 = f50 490 491POW_log2_by_128_hi = f51 492POW_inv_log2_by_128 = f52 493POW_rsq = f53 494POW_Yrcub = f54 495POW_log2_by_128_lo = f55 496 497POW_v6 = f56 498POW_xsq = f57 499POW_v4 = f58 500POW_v2 = f59 501POW_T = f60 502 503POW_Tt = f61 504POW_RSHF = f62 505POW_v21ps = f63 506POW_s4 = f64 507POW_twoV = f65 508 509POW_U = f66 510POW_G = f67 511POW_delta = f68 512POW_v3 = f69 513POW_V = f70 514 515POW_p = f71 516POW_Z1 = f72 517POW_e3 = f73 518POW_e2 = f74 519POW_Z2 = f75 520 521POW_e1 = f76 522POW_W1 = f77 523POW_UmZ2 = f78 524POW_W2 = f79 525POW_Z3 = f80 526 527POW_int_W1 = f81 528POW_e12 = f82 529POW_int_W2 = f83 530POW_UmZ2pV = f84 531POW_Z3sq = f85 532 533POW_e123 = f86 534POW_N1float = f87 535POW_N2float = f88 536POW_f3 = f89 537POW_q = f90 538 539POW_s1 = f91 540POW_Nfloat = f92 541POW_s2 = f93 542POW_f2 = f94 543POW_f1 = f95 544 545POW_T1 = f96 546POW_T2 = f97 547POW_2M = f98 548POW_s = f99 549POW_f12 = f100 550 551POW_ssq = f101 552POW_T1T2 = f102 553POW_1ps = f103 554POW_A = f104 555POW_es = f105 556 557POW_Xp1 = f106 558POW_int_K = f107 559POW_K = f108 560POW_f123 = f109 561POW_Gpr = f110 562 563POW_Y_Gpr = f111 564POW_int_Y = f112 565POW_abs_q = f114 566POW_2toM63 = f115 567 568POW_float_int_Y = f116 569POW_ftz_urm_f8 = f117 570POW_wre_urm_f8 = f118 571POW_big_neg = f119 572POW_big_pos = f120 573 574POW_GY_Z2 = f121 575POW_pYrcub_e3 = f122 576POW_d = f123 577POW_d2 = f124 578POW_poly_d_hi = f121 579POW_poly_d_lo = f122 580POW_poly_d = f121 581 582// Data tables 583//============================================================== 584 585RODATA 586 587.align 16 588 589LOCAL_OBJECT_START(pow_table_P) 590data8 0x8000F7B249FF332D, 0x0000BFFC // P_5 591data8 0xAAAAAAA9E7902C7F, 0x0000BFFC // P_3 592data8 0x80000000000018E5, 0x0000BFFD // P_1 593data8 0xb8aa3b295c17f0bc, 0x00004006 // inv_ln2_by_128 594// 595// 596data8 0x3FA5555555554A9E // Q_2 597data8 0x3F8111124F4DD9F9 // Q_3 598data8 0x3FE0000000000000 // Q_0 599data8 0x3FC5555555554733 // Q_1 600data8 0x3F56C16D9360FFA0 // Q_4 601data8 0x43e8000000000000 // Right shift constant for exp 602data8 0xc9e3b39803f2f6af, 0x00003fb7 // ln2_by_128_lo 603data8 0x0000000000000000 // pad to eliminate bank conflicts with pow_table_Q 604data8 0x0000000000000000 // pad to eliminate bank conflicts with pow_table_Q 605LOCAL_OBJECT_END(pow_table_P) 606 607LOCAL_OBJECT_START(pow_table_Q) 608data8 0x9249FE7F0DC423CF, 0x00003FFC // P_4 609data8 0xCCCCCCCC4ED2BA7F, 0x00003FFC // P_2 610data8 0xAAAAAAAAAAAAB505, 0x00003FFD // P_0 611data8 0x3fe62e42fefa39e8, 0x3cccd5e4f1d9cc02 // log2 hi lo = +6.93147e-001 612data8 0xb17217f7d1cf79ab, 0x00003ff7 // ln2_by_128_hi 613LOCAL_OBJECT_END(pow_table_Q) 614 615 616LOCAL_OBJECT_START(pow_Tt) 617data8 0x3f60040155d58800, 0x3c93bce0ce3ddd81 // log(1/frcpa(1+0/256))= +1.95503e-003 618data8 0x3f78121214586a00, 0x3cb540e0a5cfc9bc // log(1/frcpa(1+1/256))= +5.87661e-003 619data8 0x3f841929f9683200, 0x3cbdf1d57404da1f // log(1/frcpa(1+2/256))= +9.81362e-003 620data8 0x3f8c317384c75f00, 0x3c69806208c04c22 // log(1/frcpa(1+3/256))= +1.37662e-002 621data8 0x3f91a6b91ac73380, 0x3c7874daa716eb32 // log(1/frcpa(1+4/256))= +1.72376e-002 622data8 0x3f95ba9a5d9ac000, 0x3cacbb84e08d78ac // log(1/frcpa(1+5/256))= +2.12196e-002 623data8 0x3f99d2a807432580, 0x3cbcf80538b441e1 // log(1/frcpa(1+6/256))= +2.52177e-002 624data8 0x3f9d6b2725979800, 0x3c6095e5c8f8f359 // log(1/frcpa(1+7/256))= +2.87291e-002 625data8 0x3fa0c58fa19dfa80, 0x3cb4c5d4e9d0dda2 // log(1/frcpa(1+8/256))= +3.27573e-002 626data8 0x3fa2954c78cbce00, 0x3caa932b860ab8d6 // log(1/frcpa(1+9/256))= +3.62953e-002 627data8 0x3fa4a94d2da96c40, 0x3ca670452b76bbd5 // log(1/frcpa(1+10/256))= +4.03542e-002 628data8 0x3fa67c94f2d4bb40, 0x3ca84104f9941798 // log(1/frcpa(1+11/256))= +4.39192e-002 629data8 0x3fa85188b630f040, 0x3cb40a882cbf0153 // log(1/frcpa(1+12/256))= +4.74971e-002 630data8 0x3faa6b8abe73af40, 0x3c988d46e25c9059 // log(1/frcpa(1+13/256))= +5.16017e-002 631data8 0x3fac441e06f72a80, 0x3cae3e930a1a2a96 // log(1/frcpa(1+14/256))= +5.52072e-002 632data8 0x3fae1e6713606d00, 0x3c8a796f6283b580 // log(1/frcpa(1+15/256))= +5.88257e-002 633data8 0x3faffa6911ab9300, 0x3c5193070351e88a // log(1/frcpa(1+16/256))= +6.24574e-002 634data8 0x3fb0ec139c5da600, 0x3c623f2a75eb992d // log(1/frcpa(1+17/256))= +6.61022e-002 635data8 0x3fb1dbd2643d1900, 0x3ca649b2ef8927f0 // log(1/frcpa(1+18/256))= +6.97605e-002 636data8 0x3fb2cc7284fe5f00, 0x3cbc5e86599513e2 // log(1/frcpa(1+19/256))= +7.34321e-002 637data8 0x3fb3bdf5a7d1ee60, 0x3c90bd4bb69dada3 // log(1/frcpa(1+20/256))= +7.71173e-002 638data8 0x3fb4b05d7aa012e0, 0x3c54e377c9b8a54f // log(1/frcpa(1+21/256))= +8.08161e-002 639data8 0x3fb580db7ceb5700, 0x3c7fdb2f98354cde // log(1/frcpa(1+22/256))= +8.39975e-002 640data8 0x3fb674f089365a60, 0x3cb9994c9d3301c1 // log(1/frcpa(1+23/256))= +8.77219e-002 641data8 0x3fb769ef2c6b5680, 0x3caaec639db52a79 // log(1/frcpa(1+24/256))= +9.14602e-002 642data8 0x3fb85fd927506a40, 0x3c9f9f99a3cf8e25 // log(1/frcpa(1+25/256))= +9.52125e-002 643data8 0x3fb9335e5d594980, 0x3ca15c3abd47d99a // log(1/frcpa(1+26/256))= +9.84401e-002 644data8 0x3fba2b0220c8e5e0, 0x3cb4ca639adf6fc3 // log(1/frcpa(1+27/256))= +1.02219e-001 645data8 0x3fbb0004ac1a86a0, 0x3ca7cb81bf959a59 // log(1/frcpa(1+28/256))= +1.05469e-001 646data8 0x3fbbf968769fca00, 0x3cb0c646c121418e // log(1/frcpa(1+29/256))= +1.09274e-001 647data8 0x3fbccfedbfee13a0, 0x3ca0465fce24ab4b // log(1/frcpa(1+30/256))= +1.12548e-001 648data8 0x3fbda727638446a0, 0x3c82803f4e2e6603 // log(1/frcpa(1+31/256))= +1.15832e-001 649data8 0x3fbea3257fe10f60, 0x3cb986a3f2313d1a // log(1/frcpa(1+32/256))= +1.19677e-001 650data8 0x3fbf7be9fedbfde0, 0x3c97d16a6a621cf4 // log(1/frcpa(1+33/256))= +1.22985e-001 651data8 0x3fc02ab352ff25f0, 0x3c9cc6baad365600 // log(1/frcpa(1+34/256))= +1.26303e-001 652data8 0x3fc097ce579d2040, 0x3cb9ba16d329440b // log(1/frcpa(1+35/256))= +1.29633e-001 653data8 0x3fc1178e8227e470, 0x3cb7bc671683f8e6 // log(1/frcpa(1+36/256))= +1.33531e-001 654data8 0x3fc185747dbecf30, 0x3c9d1116f66d2345 // log(1/frcpa(1+37/256))= +1.36885e-001 655data8 0x3fc1f3b925f25d40, 0x3c8162c9ef939ac6 // log(1/frcpa(1+38/256))= +1.40250e-001 656data8 0x3fc2625d1e6ddf50, 0x3caad3a1ec384fc3 // log(1/frcpa(1+39/256))= +1.43627e-001 657data8 0x3fc2d1610c868130, 0x3cb3ad997036941b // log(1/frcpa(1+40/256))= +1.47015e-001 658data8 0x3fc340c597411420, 0x3cbc2308262c7998 // log(1/frcpa(1+41/256))= +1.50414e-001 659data8 0x3fc3b08b6757f2a0, 0x3cb2170d6cdf0526 // log(1/frcpa(1+42/256))= +1.53825e-001 660data8 0x3fc40dfb08378000, 0x3c9bb453c4f7b685 // log(1/frcpa(1+43/256))= +1.56677e-001 661data8 0x3fc47e74e8ca5f70, 0x3cb836a48fdfce9d // log(1/frcpa(1+44/256))= +1.60109e-001 662data8 0x3fc4ef51f6466de0, 0x3ca07a43919aa64b // log(1/frcpa(1+45/256))= +1.63553e-001 663data8 0x3fc56092e02ba510, 0x3ca85006899d97b0 // log(1/frcpa(1+46/256))= +1.67010e-001 664data8 0x3fc5d23857cd74d0, 0x3ca30a5ba6e7abbe // log(1/frcpa(1+47/256))= +1.70478e-001 665data8 0x3fc6313a37335d70, 0x3ca905586f0ac97e // log(1/frcpa(1+48/256))= +1.73377e-001 666data8 0x3fc6a399dabbd380, 0x3c9b2c6657a96684 // log(1/frcpa(1+49/256))= +1.76868e-001 667data8 0x3fc70337dd3ce410, 0x3cb50bc52f55cdd8 // log(1/frcpa(1+50/256))= +1.79786e-001 668data8 0x3fc77654128f6120, 0x3cad2eb7c9a39efe // log(1/frcpa(1+51/256))= +1.83299e-001 669data8 0x3fc7e9d82a0b0220, 0x3cba127e90393c01 // log(1/frcpa(1+52/256))= +1.86824e-001 670data8 0x3fc84a6b759f5120, 0x3cbd7fd52079f706 // log(1/frcpa(1+53/256))= +1.89771e-001 671data8 0x3fc8ab47d5f5a300, 0x3cbfae141751a3de // log(1/frcpa(1+54/256))= +1.92727e-001 672data8 0x3fc91fe490965810, 0x3cb69cf30a1c319e // log(1/frcpa(1+55/256))= +1.96286e-001 673data8 0x3fc981634011aa70, 0x3ca5bb3d208bc42a // log(1/frcpa(1+56/256))= +1.99261e-001 674data8 0x3fc9f6c407089660, 0x3ca04d68658179a0 // log(1/frcpa(1+57/256))= +2.02843e-001 675data8 0x3fca58e729348f40, 0x3c99f5411546c286 // log(1/frcpa(1+58/256))= +2.05838e-001 676data8 0x3fcabb55c31693a0, 0x3cb9a5350eb327d5 // log(1/frcpa(1+59/256))= +2.08842e-001 677data8 0x3fcb1e104919efd0, 0x3c18965fcce7c406 // log(1/frcpa(1+60/256))= +2.11855e-001 678data8 0x3fcb94ee93e367c0, 0x3cb503716da45184 // log(1/frcpa(1+61/256))= +2.15483e-001 679data8 0x3fcbf851c0675550, 0x3cbdf1b3f7ab5378 // log(1/frcpa(1+62/256))= +2.18516e-001 680data8 0x3fcc5c0254bf23a0, 0x3ca7aab9ed0b1d7b // log(1/frcpa(1+63/256))= +2.21558e-001 681data8 0x3fccc000c9db3c50, 0x3c92a7a2a850072a // log(1/frcpa(1+64/256))= +2.24609e-001 682data8 0x3fcd244d99c85670, 0x3c9f6019120edf4c // log(1/frcpa(1+65/256))= +2.27670e-001 683data8 0x3fcd88e93fb2f450, 0x3c6affb96815e081 // log(1/frcpa(1+66/256))= +2.30741e-001 684data8 0x3fcdedd437eaef00, 0x3c72553595897976 // log(1/frcpa(1+67/256))= +2.33820e-001 685data8 0x3fce530effe71010, 0x3c90913b020fa182 // log(1/frcpa(1+68/256))= +2.36910e-001 686data8 0x3fceb89a1648b970, 0x3c837ba4045bfd25 // log(1/frcpa(1+69/256))= +2.40009e-001 687data8 0x3fcf1e75fadf9bd0, 0x3cbcea6d13e0498d // log(1/frcpa(1+70/256))= +2.43117e-001 688data8 0x3fcf84a32ead7c30, 0x3ca5e3a67b3c6d77 // log(1/frcpa(1+71/256))= +2.46235e-001 689data8 0x3fcfeb2233ea07c0, 0x3cba0c6f0049c5a6 // log(1/frcpa(1+72/256))= +2.49363e-001 690data8 0x3fd028f9c7035c18, 0x3cb0a30b06677ff6 // log(1/frcpa(1+73/256))= +2.52501e-001 691data8 0x3fd05c8be0d96358, 0x3ca0f1c77ccb5865 // log(1/frcpa(1+74/256))= +2.55649e-001 692data8 0x3fd085eb8f8ae790, 0x3cbd513f45fe7a97 // log(1/frcpa(1+75/256))= +2.58174e-001 693data8 0x3fd0b9c8e32d1910, 0x3c927449047ca006 // log(1/frcpa(1+76/256))= +2.61339e-001 694data8 0x3fd0edd060b78080, 0x3c89b52d8435f53e // log(1/frcpa(1+77/256))= +2.64515e-001 695data8 0x3fd122024cf00638, 0x3cbdd976fabda4bd // log(1/frcpa(1+78/256))= +2.67701e-001 696data8 0x3fd14be2927aecd0, 0x3cb02f90ad0bc471 // log(1/frcpa(1+79/256))= +2.70257e-001 697data8 0x3fd180618ef18ad8, 0x3cbd003792c71a98 // log(1/frcpa(1+80/256))= +2.73461e-001 698data8 0x3fd1b50bbe2fc638, 0x3ca9ae64c6403ead // log(1/frcpa(1+81/256))= +2.76675e-001 699data8 0x3fd1df4cc7cf2428, 0x3cb43f0455f7e395 // log(1/frcpa(1+82/256))= +2.79254e-001 700data8 0x3fd214456d0eb8d0, 0x3cb0fbd748d75d30 // log(1/frcpa(1+83/256))= +2.82487e-001 701data8 0x3fd23ec5991eba48, 0x3c906edd746b77e2 // log(1/frcpa(1+84/256))= +2.85081e-001 702data8 0x3fd2740d9f870af8, 0x3ca9802e6a00a670 // log(1/frcpa(1+85/256))= +2.88333e-001 703data8 0x3fd29ecdabcdfa00, 0x3cacecef70890cfa // log(1/frcpa(1+86/256))= +2.90943e-001 704data8 0x3fd2d46602adcce8, 0x3cb97911955f3521 // log(1/frcpa(1+87/256))= +2.94214e-001 705data8 0x3fd2ff66b04ea9d0, 0x3cb12dabe191d1c9 // log(1/frcpa(1+88/256))= +2.96838e-001 706data8 0x3fd335504b355a30, 0x3cbdf9139df924ec // log(1/frcpa(1+89/256))= +3.00129e-001 707data8 0x3fd360925ec44f58, 0x3cb253e68977a1e3 // log(1/frcpa(1+90/256))= +3.02769e-001 708data8 0x3fd38bf1c3337e70, 0x3cb3d283d2a2da21 // log(1/frcpa(1+91/256))= +3.05417e-001 709data8 0x3fd3c25277333180, 0x3cadaa5b035eae27 // log(1/frcpa(1+92/256))= +3.08735e-001 710data8 0x3fd3edf463c16838, 0x3cb983d680d3c108 // log(1/frcpa(1+93/256))= +3.11399e-001 711data8 0x3fd419b423d5e8c0, 0x3cbc86dd921c139d // log(1/frcpa(1+94/256))= +3.14069e-001 712data8 0x3fd44591e0539f48, 0x3c86a76d6dc2782e // log(1/frcpa(1+95/256))= +3.16746e-001 713data8 0x3fd47c9175b6f0a8, 0x3cb59a2e013c6b5f // log(1/frcpa(1+96/256))= +3.20103e-001 714data8 0x3fd4a8b341552b08, 0x3c93f1e86e468694 // log(1/frcpa(1+97/256))= +3.22797e-001 715data8 0x3fd4d4f390890198, 0x3cbf5e4ea7c5105a // log(1/frcpa(1+98/256))= +3.25498e-001 716data8 0x3fd501528da1f960, 0x3cbf58da53e9ad10 // log(1/frcpa(1+99/256))= +3.28206e-001 717data8 0x3fd52dd06347d4f0, 0x3cb98a28cebf6eef // log(1/frcpa(1+100/256))= +3.30921e-001 718data8 0x3fd55a6d3c7b8a88, 0x3c9c76b67c2d1fd4 // log(1/frcpa(1+101/256))= +3.33644e-001 719data8 0x3fd5925d2b112a58, 0x3c9029616a4331b8 // log(1/frcpa(1+102/256))= +3.37058e-001 720data8 0x3fd5bf406b543db0, 0x3c9fb8292ecfc820 // log(1/frcpa(1+103/256))= +3.39798e-001 721data8 0x3fd5ec433d5c35a8, 0x3cb71a1229d17eec // log(1/frcpa(1+104/256))= +3.42545e-001 722data8 0x3fd61965cdb02c18, 0x3cbba94fe1dbb8d2 // log(1/frcpa(1+105/256))= +3.45300e-001 723data8 0x3fd646a84935b2a0, 0x3c9ee496d2c9ae57 // log(1/frcpa(1+106/256))= +3.48063e-001 724data8 0x3fd6740add31de90, 0x3cb1da3a6c7a9dfd // log(1/frcpa(1+107/256))= +3.50833e-001 725data8 0x3fd6a18db74a58c0, 0x3cb494c257add8dc // log(1/frcpa(1+108/256))= +3.53610e-001 726data8 0x3fd6cf31058670e8, 0x3cb0b244a70a8da9 // log(1/frcpa(1+109/256))= +3.56396e-001 727data8 0x3fd6f180e852f0b8, 0x3c9db7aefa866720 // log(1/frcpa(1+110/256))= +3.58490e-001 728data8 0x3fd71f5d71b894e8, 0x3cbe91c4bf324957 // log(1/frcpa(1+111/256))= +3.61289e-001 729data8 0x3fd74d5aefd66d58, 0x3cb06b3d9bfac023 // log(1/frcpa(1+112/256))= +3.64096e-001 730data8 0x3fd77b79922bd378, 0x3cb727d8804491f4 // log(1/frcpa(1+113/256))= +3.66911e-001 731data8 0x3fd7a9b9889f19e0, 0x3ca2ef22df5bc543 // log(1/frcpa(1+114/256))= +3.69734e-001 732data8 0x3fd7d81b037eb6a0, 0x3cb8fd3ba07a7ece // log(1/frcpa(1+115/256))= +3.72565e-001 733data8 0x3fd8069e33827230, 0x3c8bd1e25866e61a // log(1/frcpa(1+116/256))= +3.75404e-001 734data8 0x3fd82996d3ef8bc8, 0x3ca5aab9f5928928 // log(1/frcpa(1+117/256))= +3.77538e-001 735data8 0x3fd85855776dcbf8, 0x3ca56f33337789d6 // log(1/frcpa(1+118/256))= +3.80391e-001 736data8 0x3fd8873658327cc8, 0x3cbb8ef0401db49d // log(1/frcpa(1+119/256))= +3.83253e-001 737data8 0x3fd8aa75973ab8c8, 0x3cbb9961f509a680 // log(1/frcpa(1+120/256))= +3.85404e-001 738data8 0x3fd8d992dc8824e0, 0x3cb220512a53732d // log(1/frcpa(1+121/256))= +3.88280e-001 739data8 0x3fd908d2ea7d9510, 0x3c985f0e513bfb5c // log(1/frcpa(1+122/256))= +3.91164e-001 740data8 0x3fd92c59e79c0e50, 0x3cb82e073fd30d63 // log(1/frcpa(1+123/256))= +3.93332e-001 741data8 0x3fd95bd750ee3ed0, 0x3ca4aa7cdb6dd8a8 // log(1/frcpa(1+124/256))= +3.96231e-001 742data8 0x3fd98b7811a3ee58, 0x3caa93a5b660893e // log(1/frcpa(1+125/256))= +3.99138e-001 743data8 0x3fd9af47f33d4068, 0x3cac294b3b3190ba // log(1/frcpa(1+126/256))= +4.01323e-001 744data8 0x3fd9df270c1914a0, 0x3cbe1a58fd0cd67e // log(1/frcpa(1+127/256))= +4.04245e-001 745data8 0x3fda0325ed14fda0, 0x3cb1efa7950fb57e // log(1/frcpa(1+128/256))= +4.06442e-001 746data8 0x3fda33440224fa78, 0x3c8915fe75e7d477 // log(1/frcpa(1+129/256))= +4.09379e-001 747data8 0x3fda57725e80c380, 0x3ca72bd1062b1b7f // log(1/frcpa(1+130/256))= +4.11587e-001 748data8 0x3fda87d0165dd198, 0x3c91f7845f58dbad // log(1/frcpa(1+131/256))= +4.14539e-001 749data8 0x3fdaac2e6c03f890, 0x3cb6f237a911c509 // log(1/frcpa(1+132/256))= +4.16759e-001 750data8 0x3fdadccc6fdf6a80, 0x3c90ddc4b7687169 // log(1/frcpa(1+133/256))= +4.19726e-001 751data8 0x3fdb015b3eb1e790, 0x3c692dd7d90e1e8e // log(1/frcpa(1+134/256))= +4.21958e-001 752data8 0x3fdb323a3a635948, 0x3c6f85655cbe14de // log(1/frcpa(1+135/256))= +4.24941e-001 753data8 0x3fdb56fa04462908, 0x3c95252d841994de // log(1/frcpa(1+136/256))= +4.27184e-001 754data8 0x3fdb881aa659bc90, 0x3caa53a745a3642f // log(1/frcpa(1+137/256))= +4.30182e-001 755data8 0x3fdbad0bef3db160, 0x3cb32f2540dcc16a // log(1/frcpa(1+138/256))= +4.32437e-001 756data8 0x3fdbd21297781c28, 0x3cbd8e891e106f1d // log(1/frcpa(1+139/256))= +4.34697e-001 757data8 0x3fdc039236f08818, 0x3c809435af522ba7 // log(1/frcpa(1+140/256))= +4.37718e-001 758data8 0x3fdc28cb1e4d32f8, 0x3cb3944752fbd81e // log(1/frcpa(1+141/256))= +4.39990e-001 759data8 0x3fdc4e19b84723c0, 0x3c9a465260cd3fe5 // log(1/frcpa(1+142/256))= +4.42267e-001 760data8 0x3fdc7ff9c74554c8, 0x3c92447d5b6ca369 // log(1/frcpa(1+143/256))= +4.45311e-001 761data8 0x3fdca57b64e9db00, 0x3cb44344a8a00c82 // log(1/frcpa(1+144/256))= +4.47600e-001 762data8 0x3fdccb130a5ceba8, 0x3cbefaddfb97b73f // log(1/frcpa(1+145/256))= +4.49895e-001 763data8 0x3fdcf0c0d18f3268, 0x3cbd3e7bfee57898 // log(1/frcpa(1+146/256))= +4.52194e-001 764data8 0x3fdd232075b5a200, 0x3c9222599987447c // log(1/frcpa(1+147/256))= +4.55269e-001 765data8 0x3fdd490246defa68, 0x3cabafe9a767a80d // log(1/frcpa(1+148/256))= +4.57581e-001 766data8 0x3fdd6efa918d25c8, 0x3cb58a2624e1c6fd // log(1/frcpa(1+149/256))= +4.59899e-001 767data8 0x3fdd9509707ae528, 0x3cbdc3babce578e7 // log(1/frcpa(1+150/256))= +4.62221e-001 768data8 0x3fddbb2efe92c550, 0x3cb0ac0943c434a4 // log(1/frcpa(1+151/256))= +4.64550e-001 769data8 0x3fddee2f3445e4a8, 0x3cbba9d07ce820e8 // log(1/frcpa(1+152/256))= +4.67663e-001 770data8 0x3fde148a1a2726c8, 0x3cb6537e3375b205 // log(1/frcpa(1+153/256))= +4.70004e-001 771data8 0x3fde3afc0a49ff38, 0x3cbfed5518dbc20e // log(1/frcpa(1+154/256))= +4.72350e-001 772data8 0x3fde6185206d5168, 0x3cb6572601f73d5c // log(1/frcpa(1+155/256))= +4.74702e-001 773data8 0x3fde882578823d50, 0x3c9b24abd4584d1a // log(1/frcpa(1+156/256))= +4.77060e-001 774data8 0x3fdeaedd2eac9908, 0x3cb0ceb5e4d2c8f7 // log(1/frcpa(1+157/256))= +4.79423e-001 775data8 0x3fded5ac5f436be0, 0x3ca72f21f1f5238e // log(1/frcpa(1+158/256))= +4.81792e-001 776data8 0x3fdefc9326d16ab8, 0x3c85081a1639a45c // log(1/frcpa(1+159/256))= +4.84166e-001 777data8 0x3fdf2391a21575f8, 0x3cbf11015bdd297a // log(1/frcpa(1+160/256))= +4.86546e-001 778data8 0x3fdf4aa7ee031928, 0x3cb3795bc052a2d1 // log(1/frcpa(1+161/256))= +4.88932e-001 779data8 0x3fdf71d627c30bb0, 0x3c35c61f0f5a88f3 // log(1/frcpa(1+162/256))= +4.91323e-001 780data8 0x3fdf991c6cb3b378, 0x3c97d99419be6028 // log(1/frcpa(1+163/256))= +4.93720e-001 781data8 0x3fdfc07ada69a908, 0x3cbfe9341ded70b1 // log(1/frcpa(1+164/256))= +4.96123e-001 782data8 0x3fdfe7f18eb03d38, 0x3cb85718a640c33f // log(1/frcpa(1+165/256))= +4.98532e-001 783data8 0x3fe007c053c5002c, 0x3cb3addc9c065f09 // log(1/frcpa(1+166/256))= +5.00946e-001 784data8 0x3fe01b942198a5a0, 0x3c9d5aa4c77da6ac // log(1/frcpa(1+167/256))= +5.03367e-001 785data8 0x3fe02f74400c64e8, 0x3cb5a0ee4450ef52 // log(1/frcpa(1+168/256))= +5.05793e-001 786data8 0x3fe04360be7603ac, 0x3c9dd00c35630fe0 // log(1/frcpa(1+169/256))= +5.08225e-001 787data8 0x3fe05759ac47fe30, 0x3cbd063e1f0bd82c // log(1/frcpa(1+170/256))= +5.10663e-001 788data8 0x3fe06b5f1911cf50, 0x3cae8da674af5289 // log(1/frcpa(1+171/256))= +5.13107e-001 789data8 0x3fe078bf0533c568, 0x3c62241edf5fd1f7 // log(1/frcpa(1+172/256))= +5.14740e-001 790data8 0x3fe08cd9687e7b0c, 0x3cb3007febcca227 // log(1/frcpa(1+173/256))= +5.17194e-001 791data8 0x3fe0a10074cf9018, 0x3ca496e84603816b // log(1/frcpa(1+174/256))= +5.19654e-001 792data8 0x3fe0b5343a234474, 0x3cb46098d14fc90a // log(1/frcpa(1+175/256))= +5.22120e-001 793data8 0x3fe0c974c89431cc, 0x3cac0a7cdcbb86c6 // log(1/frcpa(1+176/256))= +5.24592e-001 794data8 0x3fe0ddc2305b9884, 0x3cb2f753210410ff // log(1/frcpa(1+177/256))= +5.27070e-001 795data8 0x3fe0eb524bafc918, 0x3c88affd6682229e // log(1/frcpa(1+178/256))= +5.28726e-001 796data8 0x3fe0ffb54213a474, 0x3cadeefbab9af993 // log(1/frcpa(1+179/256))= +5.31214e-001 797data8 0x3fe114253da97d9c, 0x3cbaf1c2b8bc160a // log(1/frcpa(1+180/256))= +5.33709e-001 798data8 0x3fe128a24f1d9afc, 0x3cb9cf4df375e650 // log(1/frcpa(1+181/256))= +5.36210e-001 799data8 0x3fe1365252bf0864, 0x3c985a621d4be111 // log(1/frcpa(1+182/256))= +5.37881e-001 800data8 0x3fe14ae558b4a92c, 0x3ca104c4aa8977d1 // log(1/frcpa(1+183/256))= +5.40393e-001 801data8 0x3fe15f85a19c7658, 0x3cbadf26e540f375 // log(1/frcpa(1+184/256))= +5.42910e-001 802data8 0x3fe16d4d38c119f8, 0x3cb3aea11caec416 // log(1/frcpa(1+185/256))= +5.44592e-001 803data8 0x3fe18203c20dd130, 0x3cba82d1211d1d6d // log(1/frcpa(1+186/256))= +5.47121e-001 804data8 0x3fe196c7bc4b1f38, 0x3cb6267acc4f4f4a // log(1/frcpa(1+187/256))= +5.49656e-001 805data8 0x3fe1a4a738b7a33c, 0x3c858930213c987d // log(1/frcpa(1+188/256))= +5.51349e-001 806data8 0x3fe1b981c0c9653c, 0x3c9bc2a4a30f697b // log(1/frcpa(1+189/256))= +5.53895e-001 807data8 0x3fe1ce69e8bb1068, 0x3cb7ae6199cf2a00 // log(1/frcpa(1+190/256))= +5.56447e-001 808data8 0x3fe1dc619de06944, 0x3c6b50bb38388177 // log(1/frcpa(1+191/256))= +5.58152e-001 809data8 0x3fe1f160a2ad0da0, 0x3cbd05b2778a5e1d // log(1/frcpa(1+192/256))= +5.60715e-001 810data8 0x3fe2066d7740737c, 0x3cb32e828f9c6bd6 // log(1/frcpa(1+193/256))= +5.63285e-001 811data8 0x3fe2147dba47a390, 0x3cbd579851b8b672 // log(1/frcpa(1+194/256))= +5.65001e-001 812data8 0x3fe229a1bc5ebac0, 0x3cbb321be5237ce8 // log(1/frcpa(1+195/256))= +5.67582e-001 813data8 0x3fe237c1841a502c, 0x3cb3b56e0915ea64 // log(1/frcpa(1+196/256))= +5.69306e-001 814data8 0x3fe24cfce6f80d98, 0x3cb34a4d1a422919 // log(1/frcpa(1+197/256))= +5.71898e-001 815data8 0x3fe25b2c55cd5760, 0x3cb237401ea5015e // log(1/frcpa(1+198/256))= +5.73630e-001 816data8 0x3fe2707f4d5f7c40, 0x3c9d30f20acc8341 // log(1/frcpa(1+199/256))= +5.76233e-001 817data8 0x3fe285e0842ca380, 0x3cbc4d866d5f21c0 // log(1/frcpa(1+200/256))= +5.78842e-001 818data8 0x3fe294294708b770, 0x3cb85e14d5dc54fa // log(1/frcpa(1+201/256))= +5.80586e-001 819data8 0x3fe2a9a2670aff0c, 0x3c7e6f8f468bbf91 // log(1/frcpa(1+202/256))= +5.83207e-001 820data8 0x3fe2b7fb2c8d1cc0, 0x3c930ffcf63c8b65 // log(1/frcpa(1+203/256))= +5.84959e-001 821data8 0x3fe2c65a6395f5f4, 0x3ca0afe20b53d2d2 // log(1/frcpa(1+204/256))= +5.86713e-001 822data8 0x3fe2dbf557b0df40, 0x3cb646be1188fbc9 // log(1/frcpa(1+205/256))= +5.89350e-001 823data8 0x3fe2ea64c3f97654, 0x3c96516fa8df33b2 // log(1/frcpa(1+206/256))= +5.91113e-001 824data8 0x3fe3001823684d70, 0x3cb96d64e16d1360 // log(1/frcpa(1+207/256))= +5.93762e-001 825data8 0x3fe30e97e9a8b5cc, 0x3c98ef96bc97cca0 // log(1/frcpa(1+208/256))= +5.95531e-001 826data8 0x3fe32463ebdd34e8, 0x3caef1dc9a56c1bf // log(1/frcpa(1+209/256))= +5.98192e-001 827data8 0x3fe332f4314ad794, 0x3caa4f0ac5d5fa11 // log(1/frcpa(1+210/256))= +5.99970e-001 828data8 0x3fe348d90e7464cc, 0x3cbe7889f0516acd // log(1/frcpa(1+211/256))= +6.02643e-001 829data8 0x3fe35779f8c43d6c, 0x3ca96bbab7245411 // log(1/frcpa(1+212/256))= +6.04428e-001 830data8 0x3fe36621961a6a98, 0x3ca31f32262db9fb // log(1/frcpa(1+213/256))= +6.06217e-001 831data8 0x3fe37c299f3c3668, 0x3cb15c72c107ee29 // log(1/frcpa(1+214/256))= +6.08907e-001 832data8 0x3fe38ae2171976e4, 0x3cba42a2554b2dd4 // log(1/frcpa(1+215/256))= +6.10704e-001 833data8 0x3fe399a157a603e4, 0x3cb99c62286d8919 // log(1/frcpa(1+216/256))= +6.12504e-001 834data8 0x3fe3afccfe77b9d0, 0x3ca11048f96a43bd // log(1/frcpa(1+217/256))= +6.15210e-001 835data8 0x3fe3be9d503533b4, 0x3ca4022f47588c3e // log(1/frcpa(1+218/256))= +6.17018e-001 836data8 0x3fe3cd7480b4a8a0, 0x3cb4ba7afc2dc56a // log(1/frcpa(1+219/256))= +6.18830e-001 837data8 0x3fe3e3c43918f76c, 0x3c859673d064b8ba // log(1/frcpa(1+220/256))= +6.21554e-001 838data8 0x3fe3f2acb27ed6c4, 0x3cb55c6b452a16a8 // log(1/frcpa(1+221/256))= +6.23373e-001 839data8 0x3fe4019c2125ca90, 0x3cb8c367879c5a31 // log(1/frcpa(1+222/256))= +6.25197e-001 840data8 0x3fe4181061389720, 0x3cb2c17a79c5cc6c // log(1/frcpa(1+223/256))= +6.27937e-001 841data8 0x3fe42711518df544, 0x3ca5f38d47012fc5 // log(1/frcpa(1+224/256))= +6.29769e-001 842data8 0x3fe436194e12b6bc, 0x3cb9854d65a9b426 // log(1/frcpa(1+225/256))= +6.31604e-001 843data8 0x3fe445285d68ea68, 0x3ca3ff9b3a81cd81 // log(1/frcpa(1+226/256))= +6.33442e-001 844data8 0x3fe45bcc464c8938, 0x3cb0a2d8011a6c05 // log(1/frcpa(1+227/256))= +6.36206e-001 845data8 0x3fe46aed21f117fc, 0x3c8a2be41f8e9f3d // log(1/frcpa(1+228/256))= +6.38053e-001 846data8 0x3fe47a1527e8a2d0, 0x3cba4a83594fab09 // log(1/frcpa(1+229/256))= +6.39903e-001 847data8 0x3fe489445efffcc8, 0x3cbf306a23dcbcde // log(1/frcpa(1+230/256))= +6.41756e-001 848data8 0x3fe4a018bcb69834, 0x3ca46c9285029fd1 // log(1/frcpa(1+231/256))= +6.44543e-001 849data8 0x3fe4af5a0c9d65d4, 0x3cbbc1db897580e3 // log(1/frcpa(1+232/256))= +6.46405e-001 850data8 0x3fe4bea2a5bdbe84, 0x3cb84d880d7ef775 // log(1/frcpa(1+233/256))= +6.48271e-001 851data8 0x3fe4cdf28f10ac44, 0x3cb3ec4b7893ce1f // log(1/frcpa(1+234/256))= +6.50140e-001 852data8 0x3fe4dd49cf994058, 0x3c897224d59d3408 // log(1/frcpa(1+235/256))= +6.52013e-001 853data8 0x3fe4eca86e64a680, 0x3cbccf620f24f0cd // log(1/frcpa(1+236/256))= +6.53889e-001 854data8 0x3fe503c43cd8eb68, 0x3c3f872c65971084 // log(1/frcpa(1+237/256))= +6.56710e-001 855data8 0x3fe513356667fc54, 0x3cb9ca64cc3d52c8 // log(1/frcpa(1+238/256))= +6.58595e-001 856data8 0x3fe522ae0738a3d4, 0x3cbe708164c75968 // log(1/frcpa(1+239/256))= +6.60483e-001 857data8 0x3fe5322e26867854, 0x3cb9988ba4aea615 // log(1/frcpa(1+240/256))= +6.62376e-001 858data8 0x3fe541b5cb979808, 0x3ca1662e3a6b95f5 // log(1/frcpa(1+241/256))= +6.64271e-001 859data8 0x3fe55144fdbcbd60, 0x3cb3acd4ca45c1e0 // log(1/frcpa(1+242/256))= +6.66171e-001 860data8 0x3fe560dbc45153c4, 0x3cb4988947959fed // log(1/frcpa(1+243/256))= +6.68074e-001 861data8 0x3fe5707a26bb8c64, 0x3cb3017fe6607ba9 // log(1/frcpa(1+244/256))= +6.69980e-001 862data8 0x3fe587f60ed5b8fc, 0x3cbe7a3266366ed4 // log(1/frcpa(1+245/256))= +6.72847e-001 863data8 0x3fe597a7977c8f30, 0x3ca1e12b9959a90e // log(1/frcpa(1+246/256))= +6.74763e-001 864data8 0x3fe5a760d634bb88, 0x3cb7c365e53d9602 // log(1/frcpa(1+247/256))= +6.76682e-001 865data8 0x3fe5b721d295f10c, 0x3cb716c2551ccbf0 // log(1/frcpa(1+248/256))= +6.78605e-001 866data8 0x3fe5c6ea94431ef8, 0x3ca02b2ed0e28261 // log(1/frcpa(1+249/256))= +6.80532e-001 867data8 0x3fe5d6bb22ea86f4, 0x3caf43a8bbb2f974 // log(1/frcpa(1+250/256))= +6.82462e-001 868data8 0x3fe5e6938645d38c, 0x3cbcedc98821b333 // log(1/frcpa(1+251/256))= +6.84397e-001 869data8 0x3fe5f673c61a2ed0, 0x3caa385eef5f2789 // log(1/frcpa(1+252/256))= +6.86335e-001 870data8 0x3fe6065bea385924, 0x3cb11624f165c5b4 // log(1/frcpa(1+253/256))= +6.88276e-001 871data8 0x3fe6164bfa7cc068, 0x3cbad884f87073fa // log(1/frcpa(1+254/256))= +6.90222e-001 872data8 0x3fe62643fecf9740, 0x3cb78c51da12f4df // log(1/frcpa(1+255/256))= +6.92171e-001 873LOCAL_OBJECT_END(pow_Tt) 874 875 876// Table 1 is 2^(index_1/128) where 877// index_1 goes from 0 to 15 878LOCAL_OBJECT_START(pow_tbl1) 879data8 0x8000000000000000 , 0x00003FFF 880data8 0x80B1ED4FD999AB6C , 0x00003FFF 881data8 0x8164D1F3BC030773 , 0x00003FFF 882data8 0x8218AF4373FC25EC , 0x00003FFF 883data8 0x82CD8698AC2BA1D7 , 0x00003FFF 884data8 0x8383594EEFB6EE37 , 0x00003FFF 885data8 0x843A28C3ACDE4046 , 0x00003FFF 886data8 0x84F1F656379C1A29 , 0x00003FFF 887data8 0x85AAC367CC487B15 , 0x00003FFF 888data8 0x8664915B923FBA04 , 0x00003FFF 889data8 0x871F61969E8D1010 , 0x00003FFF 890data8 0x87DB357FF698D792 , 0x00003FFF 891data8 0x88980E8092DA8527 , 0x00003FFF 892data8 0x8955EE03618E5FDD , 0x00003FFF 893data8 0x8A14D575496EFD9A , 0x00003FFF 894data8 0x8AD4C6452C728924 , 0x00003FFF 895LOCAL_OBJECT_END(pow_tbl1) 896 897 898// Table 2 is 2^(index_1/8) where 899// index_2 goes from 0 to 7 900LOCAL_OBJECT_START(pow_tbl2) 901data8 0x8000000000000000 , 0x00003FFF 902data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF 903data8 0x9837F0518DB8A96F , 0x00003FFF 904data8 0xA5FED6A9B15138EA , 0x00003FFF 905data8 0xB504F333F9DE6484 , 0x00003FFF 906data8 0xC5672A115506DADD , 0x00003FFF 907data8 0xD744FCCAD69D6AF4 , 0x00003FFF 908data8 0xEAC0C6E7DD24392F , 0x00003FFF 909LOCAL_OBJECT_END(pow_tbl2) 910 911.section .text 912WEAK_LIBM_ENTRY(pow) 913 914// Get exponent of x. Will be used to calculate K. 915{ .mfi 916 getf.exp pow_GR_signexp_X = f8 917 fms.s1 POW_Xm1 = f8,f1,f1 // Will be used for r1 if x>0 918 mov pow_GR_17ones = 0x1FFFF 919} 920{ .mfi 921 addl pow_AD_P = @ltoff(pow_table_P), gp 922 fma.s1 POW_Xp1 = f8,f1,f1 // Will be used for r1 if x<0 923 nop.i 999 924;; 925} 926 927// Get significand of x. Will be used to get index to fetch T, Tt. 928{ .mfi 929 getf.sig pow_GR_sig_X = f8 930 frcpa.s1 POW_B, p6 = f1,f8 931 nop.i 999 932} 933{ .mfi 934 ld8 pow_AD_P = [pow_AD_P] 935 fma.s1 POW_NORM_X = f8,f1,f0 936 mov pow_GR_exp_2tom8 = 0xFFF7 937} 938;; 939 940// p13 = TRUE ==> X is unorm 941// DOUBLE 0x10033 exponent limit at which y is an integer 942{ .mfi 943 nop.m 999 944 fclass.m p13,p0 = f8, 0x0b // Test for x unorm 945 addl pow_GR_10033 = 0x10033, r0 946} 947{ .mfi 948 mov pow_GR_16ones = 0xFFFF 949 fma.s1 POW_NORM_Y = f9,f1,f0 950 nop.i 999 951} 952;; 953 954// p14 = TRUE ==> X is ZERO 955{ .mfi 956 adds pow_AD_Tt = pow_Tt - pow_table_P, pow_AD_P 957 fclass.m p14,p0 = f8, 0x07 958 and pow_GR_exp_X = pow_GR_signexp_X, pow_GR_17ones 959} 960{ .mfi 961 adds pow_AD_Q = pow_table_Q - pow_table_P, pow_AD_P 962 nop.f 999 963 nop.i 999 964} 965;; 966 967{ .mfi 968 ldfe POW_P5 = [pow_AD_P], 16 969 fcmp.lt.s1 p8,p9 = f8, f0 // Test for x<0 970 nop.i 999 971} 972{ .mib 973 ldfe POW_P4 = [pow_AD_Q], 16 974 sub pow_GR_true_exp_X = pow_GR_exp_X, pow_GR_16ones 975(p13) br.cond.spnt POW_X_DENORM 976} 977;; 978 979// Continue normal and denormal paths here 980POW_COMMON: 981// p11 = TRUE ==> Y is a NAN 982{ .mfi 983 ldfe POW_P3 = [pow_AD_P], 16 984 fclass.m p11,p0 = f9, 0xc3 985 nop.i 999 986} 987{ .mfi 988 ldfe POW_P2 = [pow_AD_Q], 16 989 nop.f 999 990 mov pow_GR_y_zero = 0 991} 992;; 993 994// Note POW_Xm1 and POW_r1 are used interchangably 995{ .mfi 996 alloc r32=ar.pfs,2,19,4,0 997 fms.s1 POW_r = POW_B, POW_NORM_X,f1 998 nop.i 999 999} 1000{ .mfi 1001 setf.sig POW_int_K = pow_GR_true_exp_X 1002(p8) fnma.s1 POW_Xm1 = POW_Xp1,f1,f0 1003 nop.i 999 1004} 1005;; 1006 1007// p12 = TRUE if Y is ZERO 1008// Compute xsq to decide later if |x|=1 1009{ .mfi 1010 ldfe POW_P1 = [pow_AD_P], 16 1011 fclass.m p12,p0 = f9, 0x07 1012 shl pow_GR_offset = pow_GR_sig_X, 1 1013} 1014{ .mfb 1015 ldfe POW_P0 = [pow_AD_Q], 16 1016 fma.s1 POW_xsq = POW_NORM_X, POW_NORM_X, f0 1017(p11) br.cond.spnt POW_Y_NAN // Branch if y=nan 1018} 1019;; 1020 1021// Get exponent of |x|-1 to use in comparison to 2^-8 1022{ .mfi 1023 getf.exp pow_GR_signexp_Xm1 = POW_Xm1 1024 fcvt.fx.s1 POW_int_Y = POW_NORM_Y 1025 shr.u pow_GR_offset = pow_GR_offset,56 1026} 1027;; 1028 1029// p11 = TRUE ==> X is a NAN 1030{ .mfi 1031 ldfpd POW_log2_hi, POW_log2_lo = [pow_AD_Q], 16 1032 fclass.m p11,p0 = f8, 0xc3 1033 shladd pow_AD_Tt = pow_GR_offset, 4, pow_AD_Tt 1034} 1035{ .mfi 1036 ldfe POW_inv_log2_by_128 = [pow_AD_P], 16 1037 fma.s1 POW_delta = f0,f0,f0 // delta=0 in case |x| near 1 1038(p12) mov pow_GR_y_zero = 1 1039} 1040;; 1041 1042{ .mfi 1043 ldfpd POW_Q2, POW_Q3 = [pow_AD_P], 16 1044 fma.s1 POW_G = f0,f0,f0 // G=0 in case |x| near 1 1045 and pow_GR_exp_Xm1 = pow_GR_signexp_Xm1, pow_GR_17ones 1046} 1047;; 1048 1049// Determine if we will use the |x| near 1 path (p6) or normal path (p7) 1050{ .mfi 1051 getf.exp pow_GR_signexp_Y = POW_NORM_Y 1052 nop.f 999 1053 cmp.lt p6,p7 = pow_GR_exp_Xm1, pow_GR_exp_2tom8 1054} 1055{ .mfb 1056 ldfpd POW_T, POW_Tt = [pow_AD_Tt], 16 1057 fma.s1 POW_rsq = POW_r, POW_r,f0 1058(p11) br.cond.spnt POW_X_NAN // Branch if x=nan and y not nan 1059} 1060;; 1061 1062// If on the x near 1 path, assign r1 to r and r1*r1 to rsq 1063{ .mfi 1064 ldfpd POW_Q0_half, POW_Q1 = [pow_AD_P], 16 1065(p6) fma.s1 POW_r = POW_r1, f1, f0 1066 nop.i 999 1067} 1068{ .mfb 1069 nop.m 999 1070(p6) fma.s1 POW_rsq = POW_r1, POW_r1, f0 1071(p14) br.cond.spnt POW_X_0 // Branch if x zero and y not nan 1072} 1073;; 1074 1075{ .mfi 1076 ldfpd POW_Q4, POW_RSHF = [pow_AD_P], 16 1077(p7) fma.s1 POW_v6 = POW_r, POW_P5, POW_P4 1078 nop.i 999 1079} 1080{ .mfi 1081 mov pow_GR_exp_2toM63 = 0xffc0 // Exponent of 2^-63 1082(p6) fma.s1 POW_v6 = POW_r1, POW_P5, POW_P4 1083 nop.i 999 1084} 1085;; 1086 1087{ .mfi 1088 setf.exp POW_2toM63 = pow_GR_exp_2toM63 // Form 2^-63 for test of q 1089(p7) fma.s1 POW_v4 = POW_P3, POW_r, POW_P2 1090 nop.i 999 1091} 1092{ .mfi 1093 nop.m 999 1094(p6) fma.s1 POW_v4 = POW_P3, POW_r1, POW_P2 1095 nop.i 999 1096} 1097;; 1098 1099{ .mfi 1100 nop.m 999 1101 fcvt.xf POW_K = POW_int_K 1102 nop.i 999 1103} 1104;; 1105 1106{ .mfi 1107 getf.sig pow_GR_sig_int_Y = POW_int_Y 1108 fnma.s1 POW_twoV = POW_NORM_Y, POW_rsq,f0 1109 and pow_GR_exp_Y = pow_GR_signexp_Y, pow_GR_17ones 1110} 1111{ .mfb 1112 andcm pow_GR_sign_Y = pow_GR_signexp_Y, pow_GR_17ones 1113 fma.s1 POW_U = POW_NORM_Y,POW_r,f0 1114(p12) br.cond.spnt POW_Y_0 // Branch if y=zero, x not zero or nan 1115} 1116;; 1117 1118// p11 = TRUE ==> X is NEGATIVE but not inf 1119{ .mfi 1120 ldfe POW_log2_by_128_lo = [pow_AD_P], 16 1121 fclass.m p11,p0 = f8, 0x1a 1122 nop.i 999 1123} 1124{ .mfi 1125 ldfe POW_log2_by_128_hi = [pow_AD_Q], 16 1126 fma.s1 POW_v2 = POW_P1, POW_r, POW_P0 1127 nop.i 999 1128} 1129;; 1130 1131{ .mfi 1132 nop.m 999 1133 fcvt.xf POW_float_int_Y = POW_int_Y 1134 nop.i 999 1135} 1136{ .mfi 1137 nop.m 999 1138 fma.s1 POW_v3 = POW_v6, POW_rsq, POW_v4 1139 adds pow_AD_tbl1 = pow_tbl1 - pow_Tt, pow_AD_Q 1140} 1141;; 1142 1143{ .mfi 1144 nop.m 999 1145(p7) fma.s1 POW_delta = POW_K, POW_log2_lo, POW_Tt 1146 nop.i 999 1147} 1148{ .mfi 1149 nop.m 999 1150(p7) fma.s1 POW_G = POW_K, POW_log2_hi, POW_T 1151 adds pow_AD_tbl2 = pow_tbl2 - pow_tbl1, pow_AD_tbl1 1152} 1153;; 1154 1155{ .mfi 1156 nop.m 999 1157 fms.s1 POW_e2 = POW_NORM_Y, POW_r, POW_U 1158 nop.i 999 1159} 1160{ .mfi 1161 nop.m 999 1162 fma.s1 POW_Z2 = POW_twoV, POW_Q0_half, POW_U 1163 nop.i 999 1164} 1165;; 1166 1167{ .mfi 1168 nop.m 999 1169 fma.s1 POW_Yrcub = POW_rsq, POW_U, f0 1170 nop.i 999 1171} 1172{ .mfi 1173 nop.m 999 1174 fma.s1 POW_p = POW_rsq, POW_v3, POW_v2 1175 nop.i 999 1176} 1177;; 1178 1179// p11 = TRUE ==> X is NEGATIVE but not inf 1180// p12 = TRUE ==> X is NEGATIVE AND Y already even int 1181// p13 = TRUE ==> X is NEGATIVE AND Y possible int 1182{ .mfi 1183 nop.m 999 1184 fma.s1 POW_Z1 = POW_NORM_Y, POW_G, f0 1185(p11) cmp.gt.unc p12,p13 = pow_GR_exp_Y, pow_GR_10033 1186} 1187{ .mfi 1188 nop.m 999 1189 fma.s1 POW_Gpr = POW_G, f1, POW_r 1190 nop.i 999 1191} 1192;; 1193 1194// By adding RSHF (1.1000...*2^63) we put integer part in rightmost significand 1195{ .mfi 1196 nop.m 999 1197 fma.s1 POW_W2 = POW_Z2, POW_inv_log2_by_128, POW_RSHF 1198 nop.i 999 1199} 1200{ .mfi 1201 nop.m 999 1202 fms.s1 POW_UmZ2 = POW_U, f1, POW_Z2 1203 nop.i 999 1204} 1205;; 1206 1207{ .mfi 1208 nop.m 999 1209 fma.s1 POW_e3 = POW_NORM_Y, POW_delta, f0 1210 nop.i 999 1211} 1212;; 1213 1214{ .mfi 1215 nop.m 999 1216 fma.s1 POW_Z3 = POW_p, POW_Yrcub, f0 1217 nop.i 999 1218} 1219{ .mfi 1220 nop.m 999 1221 fma.s1 POW_GY_Z2 = POW_G, POW_NORM_Y, POW_Z2 1222 nop.i 999 1223} 1224;; 1225 1226// By adding RSHF (1.1000...*2^63) we put integer part in rightmost significand 1227{ .mfi 1228 nop.m 999 1229 fms.s1 POW_e1 = POW_NORM_Y, POW_G, POW_Z1 1230 nop.i 999 1231} 1232{ .mfi 1233 nop.m 999 1234 fma.s1 POW_W1 = POW_Z1, POW_inv_log2_by_128, POW_RSHF 1235 nop.i 999 1236} 1237;; 1238 1239// p13 = TRUE ==> X is NEGATIVE AND Y possible int 1240// p10 = TRUE ==> X is NEG and Y is an int 1241// p12 = TRUE ==> X is NEG and Y is not an int 1242{ .mfi 1243 nop.m 999 1244(p13) fcmp.eq.unc.s1 p10,p12 = POW_float_int_Y, POW_NORM_Y 1245 mov pow_GR_xneg_yodd = 0 1246} 1247{ .mfi 1248 nop.m 999 1249 fma.s1 POW_Y_Gpr = POW_NORM_Y, POW_Gpr, f0 1250 nop.i 999 1251} 1252;; 1253 1254// By subtracting RSHF we get rounded integer POW_N2float 1255{ .mfi 1256 nop.m 999 1257 fms.s1 POW_N2float = POW_W2, f1, POW_RSHF 1258 nop.i 999 1259} 1260{ .mfi 1261 nop.m 999 1262 fma.s1 POW_UmZ2pV = POW_twoV,POW_Q0_half,POW_UmZ2 1263 nop.i 999 1264} 1265;; 1266 1267{ .mfi 1268 nop.m 999 1269 fma.s1 POW_Z3sq = POW_Z3, POW_Z3, f0 1270 nop.i 999 1271} 1272{ .mfi 1273 nop.m 999 1274 fma.s1 POW_v4 = POW_Z3, POW_Q3, POW_Q2 1275 nop.i 999 1276} 1277;; 1278 1279// Extract rounded integer from rightmost significand of POW_W2 1280// By subtracting RSHF we get rounded integer POW_N1float 1281{ .mfi 1282 getf.sig pow_GR_int_W2 = POW_W2 1283 fms.s1 POW_N1float = POW_W1, f1, POW_RSHF 1284 nop.i 999 1285} 1286{ .mfi 1287 nop.m 999 1288 fma.s1 POW_v2 = POW_Z3, POW_Q1, POW_Q0_half 1289 nop.i 999 1290} 1291;; 1292 1293{ .mfi 1294 nop.m 999 1295 fnma.s1 POW_s2 = POW_N2float, POW_log2_by_128_hi, POW_Z2 1296 nop.i 999 1297} 1298{ .mfi 1299 nop.m 999 1300 fma.s1 POW_e2 = POW_e2,f1,POW_UmZ2pV 1301 nop.i 999 1302} 1303;; 1304 1305// Extract rounded integer from rightmost significand of POW_W1 1306// Test if x inf 1307{ .mfi 1308 getf.sig pow_GR_int_W1 = POW_W1 1309 fclass.m p15,p0 = POW_NORM_X, 0x23 1310 nop.i 999 1311} 1312{ .mfb 1313 nop.m 999 1314 fnma.s1 POW_f2 = POW_N2float, POW_log2_by_128_lo, f1 1315(p12) br.cond.spnt POW_X_NEG_Y_NONINT // Branch if x neg, y not integer 1316} 1317;; 1318 1319// p11 = TRUE ==> X is +1.0 1320// p12 = TRUE ==> X is NEGATIVE AND Y is an odd integer 1321{ .mfi 1322 getf.exp pow_GR_signexp_Y_Gpr = POW_Y_Gpr 1323 fcmp.eq.s1 p11,p0 = POW_NORM_X, f1 1324(p10) tbit.nz.unc p12,p0 = pow_GR_sig_int_Y,0 1325} 1326{ .mfi 1327 nop.m 999 1328 fma.s1 POW_v3 = POW_Z3sq, POW_Q4, POW_v4 1329 nop.i 999 1330} 1331;; 1332 1333{ .mfi 1334 nop.m 999 1335 fnma.s1 POW_f1 = POW_N1float, POW_log2_by_128_lo, f1 1336 nop.i 999 1337} 1338{ .mfb 1339 nop.m 999 1340 fnma.s1 POW_s1 = POW_N1float, POW_log2_by_128_hi, POW_Z1 1341(p15) br.cond.spnt POW_X_INF 1342} 1343;; 1344 1345// Test x and y and flag denormal 1346{ .mfi 1347 nop.m 999 1348 fcmp.eq.s0 p15,p0 = f8,f9 1349 nop.i 999 1350} 1351{ .mfi 1352 nop.m 999 1353 fma.s1 POW_pYrcub_e3 = POW_p, POW_Yrcub, POW_e3 1354 nop.i 999 1355} 1356;; 1357 1358{ .mfi 1359 nop.m 999 1360 fcmp.eq.s1 p7,p0 = POW_NORM_Y, f1 // Test for y=1.0 1361 nop.i 999 1362} 1363{ .mfi 1364 nop.m 999 1365 fma.s1 POW_e12 = POW_e1,f1,POW_e2 1366 nop.i 999 1367} 1368;; 1369 1370{ .mfi 1371 add pow_GR_int_N = pow_GR_int_W1, pow_GR_int_W2 1372(p11) fma.d.s0 f8 = f1,f1,f0 // If x=1, result is +1 1373 nop.i 999 1374} 1375{ .mib 1376(p12) mov pow_GR_xneg_yodd = 1 1377 nop.i 999 1378(p11) br.ret.spnt b0 // Early exit if x=1.0, result is +1 1379} 1380;; 1381 1382{ .mfi 1383 and pow_GR_index1 = 0x0f, pow_GR_int_N 1384 fma.s1 POW_q = POW_Z3sq, POW_v3, POW_v2 1385 shr pow_int_GR_M = pow_GR_int_N, 7 // M = N/128 1386} 1387{ .mib 1388 and pow_GR_index2 = 0x70, pow_GR_int_N 1389 cmp.eq p6, p0 = pow_GR_xneg_yodd, r0 1390(p7) br.ret.spnt b0 // Early exit if y=1.0, result is x 1391} 1392;; 1393 1394{ .mfi 1395 shladd pow_AD_T1 = pow_GR_index1, 4, pow_AD_tbl1 1396 fma.s1 POW_s = POW_s1, f1, POW_s2 1397 add pow_int_GR_M = pow_GR_16ones, pow_int_GR_M 1398} 1399{ .mfi 1400 add pow_AD_T2 = pow_AD_tbl2, pow_GR_index2 1401 fma.s1 POW_f12 = POW_f1, POW_f2,f0 1402 and pow_GR_exp_Y_Gpr = pow_GR_signexp_Y_Gpr, pow_GR_17ones 1403} 1404;; 1405 1406{ .mmi 1407 ldfe POW_T1 = [pow_AD_T1] 1408 ldfe POW_T2 = [pow_AD_T2] 1409 sub pow_GR_true_exp_Y_Gpr = pow_GR_exp_Y_Gpr, pow_GR_16ones 1410} 1411;; 1412 1413{ .mfi 1414 setf.exp POW_2M = pow_int_GR_M 1415 fma.s1 POW_e123 = POW_e12, f1, POW_e3 1416 nop.i 999 1417} 1418{ .mfb 1419(p6) cmp.gt p6, p0 = -11, pow_GR_true_exp_Y_Gpr 1420 fma.s1 POW_d = POW_GY_Z2, f1, POW_pYrcub_e3 1421(p6) br.cond.spnt POW_NEAR_ONE // branch if |y*log(x)| < 2^(-11) 1422} 1423;; 1424 1425{ .mfi 1426 nop.m 999 1427 fma.s1 POW_q = POW_Z3sq, POW_q, POW_Z3 1428 nop.i 999 1429} 1430;; 1431 1432// p8 TRUE ==> |Y(G + r)| >= 10 1433 1434// double 1435// -2^10 -2^9 2^9 2^10 1436// -----+-----+----+ ... +-----+-----+----- 1437// p8 | p9 | p8 1438// | | p10 | | 1439 1440// Form signexp of constants to indicate overflow 1441{ .mfi 1442 mov pow_GR_big_pos = 0x103ff 1443 fma.s1 POW_ssq = POW_s, POW_s, f0 1444 cmp.le p8,p9 = 10, pow_GR_true_exp_Y_Gpr 1445} 1446{ .mfi 1447 mov pow_GR_big_neg = 0x303ff 1448 fma.s1 POW_v4 = POW_s, POW_Q3, POW_Q2 1449 andcm pow_GR_sign_Y_Gpr = pow_GR_signexp_Y_Gpr, pow_GR_17ones 1450} 1451;; 1452 1453// Form big positive and negative constants to test for possible overflow 1454{ .mfi 1455 setf.exp POW_big_pos = pow_GR_big_pos 1456 fma.s1 POW_v2 = POW_s, POW_Q1, POW_Q0_half 1457(p9) cmp.le.unc p0,p10 = 9, pow_GR_true_exp_Y_Gpr 1458} 1459{ .mfb 1460 setf.exp POW_big_neg = pow_GR_big_neg 1461 fma.s1 POW_1ps = f1,f1,POW_s 1462(p8) br.cond.spnt POW_OVER_UNDER_X_NOT_INF 1463} 1464;; 1465 1466// f123 = f12*(e123+1) = f12*e123+f12 1467{ .mfi 1468 nop.m 999 1469 fma.s1 POW_f123 = POW_e123,POW_f12,POW_f12 1470 nop.i 999 1471} 1472;; 1473 1474{ .mfi 1475 nop.m 999 1476 fma.s1 POW_T1T2 = POW_T1, POW_T2, f0 1477 nop.i 999 1478} 1479{ .mfi 1480 nop.m 999 1481 fma.s1 POW_v3 = POW_ssq, POW_Q4, POW_v4 1482 cmp.ne p12,p13 = pow_GR_xneg_yodd, r0 1483} 1484;; 1485 1486{ .mfi 1487 nop.m 999 1488 fma.s1 POW_v21ps = POW_ssq, POW_v2, POW_1ps 1489 nop.i 999 1490} 1491{ .mfi 1492 nop.m 999 1493 fma.s1 POW_s4 = POW_ssq, POW_ssq, f0 1494 nop.i 999 1495} 1496;; 1497 1498{ .mfi 1499 nop.m 999 1500(p12) fnma.s1 POW_A = POW_2M, POW_f123, f0 1501 nop.i 999 1502} 1503{ .mfi 1504 nop.m 999 1505(p13) fma.s1 POW_A = POW_2M, POW_f123, f0 1506 cmp.eq p14,p11 = r0,r0 // Initialize p14 on, p11 off 1507} 1508;; 1509 1510{ .mfi 1511 nop.m 999 1512 fmerge.s POW_abs_q = f0, POW_q // Form |q| so can test its size 1513 nop.i 999 1514} 1515;; 1516 1517{ .mfi 1518(p10) cmp.eq p0,p14 = r0,r0 // Turn off p14 if no overflow 1519 fma.s1 POW_es = POW_s4, POW_v3, POW_v21ps 1520 nop.i 999 1521} 1522{ .mfi 1523 nop.m 999 1524 fma.s1 POW_A = POW_A, POW_T1T2, f0 1525 nop.i 999 1526} 1527;; 1528 1529{ .mfi 1530// Test for |q| < 2^-63. If so then reverse last two steps of the result 1531// to avoid monotonicity problems for results near 1.0 in round up/down/zero. 1532// p11 will be set if need to reverse the order, p14 if not. 1533 nop.m 999 1534(p10) fcmp.lt.s0 p11,p14 = POW_abs_q, POW_2toM63 // Test |q| <2^-63 1535 nop.i 999 1536} 1537;; 1538 1539.pred.rel "mutex",p11,p14 1540{ .mfi 1541 nop.m 999 1542(p14) fma.s1 POW_A = POW_A, POW_es, f0 1543 nop.i 999 1544} 1545{ .mfi 1546 nop.m 999 1547(p11) fma.s1 POW_A = POW_A, POW_q, POW_A 1548 nop.i 999 1549} 1550;; 1551 1552// Dummy op to set inexact if |q| < 2^-63 1553{ .mfi 1554 nop.m 999 1555(p11) fma.d.s0 POW_tmp = POW_A, POW_q, POW_A 1556 nop.i 999 1557} 1558;; 1559 1560{ .mfi 1561 nop.m 999 1562(p14) fma.d.s0 f8 = POW_A, POW_q, POW_A 1563 nop.i 999 1564} 1565{ .mfb 1566 nop.m 999 1567(p11) fma.d.s0 f8 = POW_A, POW_es, f0 1568(p10) br.ret.sptk b0 // Exit main branch if no over/underflow 1569} 1570;; 1571 1572// POSSIBLE_OVER_UNDER 1573// p6 = TRUE ==> Y_Gpr negative 1574// Result is already computed. We just need to know if over/underflow occurred. 1575 1576{ .mfb 1577 cmp.eq p0,p6 = pow_GR_sign_Y_Gpr, r0 1578 nop.f 999 1579(p6) br.cond.spnt POW_POSSIBLE_UNDER 1580} 1581;; 1582 1583// POSSIBLE_OVER 1584// We got an answer. 1585// overflow is a possibility, not a certainty 1586 1587 1588// We define an overflow when the answer with 1589// WRE set 1590// user-defined rounding mode 1591 1592// double 1593// Largest double is 7FE (biased double) 1594// 7FE - 3FF + FFFF = 103FE 1595// Create + largest_double_plus_ulp 1596// Create - largest_double_plus_ulp 1597// Calculate answer with WRE set. 1598 1599// single 1600// Largest single is FE (biased double) 1601// FE - 7F + FFFF = 1007E 1602// Create + largest_single_plus_ulp 1603// Create - largest_single_plus_ulp 1604// Calculate answer with WRE set. 1605 1606// Cases when answer is ldn+1 are as follows: 1607// ldn ldn+1 1608// --+----------|----------+------------ 1609// | 1610// +inf +inf -inf 1611// RN RN 1612// RZ 1613 1614// Put in s2 (td set, wre set) 1615{ .mfi 1616 nop.m 999 1617 fsetc.s2 0x7F,0x42 1618 nop.i 999 1619} 1620;; 1621 1622{ .mfi 1623 nop.m 999 1624 fma.d.s2 POW_wre_urm_f8 = POW_A, POW_q, POW_A 1625 nop.i 999 1626} 1627;; 1628 1629// Return s2 to default 1630{ .mfi 1631 nop.m 999 1632 fsetc.s2 0x7F,0x40 1633 nop.i 999 1634} 1635;; 1636 1637// p7 = TRUE ==> yes, we have an overflow 1638{ .mfi 1639 nop.m 999 1640 fcmp.ge.s1 p7, p8 = POW_wre_urm_f8, POW_big_pos 1641 nop.i 999 1642} 1643;; 1644 1645{ .mfi 1646 nop.m 999 1647(p8) fcmp.le.s1 p7, p0 = POW_wre_urm_f8, POW_big_neg 1648 nop.i 999 1649} 1650;; 1651 1652{ .mbb 1653(p7) mov pow_GR_tag = 24 1654(p7) br.cond.spnt __libm_error_region // Branch if overflow 1655 br.ret.sptk b0 // Exit if did not overflow 1656} 1657;; 1658 1659// Here if |y*log(x)| < 2^(-11) 1660// pow(x,y) ~ exp(d) ~ 1 + d + 0.5*d^2 + Q1*d^3 + Q2*d^4, where d = y*log(x) 1661.align 32 1662POW_NEAR_ONE: 1663 1664{ .mfi 1665 nop.m 999 1666 fma.s1 POW_d2 = POW_d, POW_d, f0 1667 nop.i 999 1668} 1669;; 1670 1671{ .mfi 1672 nop.m 999 1673 fma.s1 POW_poly_d_hi = POW_d, POW_Q0_half, f1 1674 nop.i 999 1675} 1676{ .mfi 1677 nop.m 999 1678 fma.s1 POW_poly_d_lo = POW_d, POW_Q2, POW_Q1 1679 nop.i 999 1680} 1681;; 1682 1683{ .mfi 1684 nop.m 999 1685 fma.s1 POW_poly_d = POW_d2, POW_poly_d_lo, POW_poly_d_hi 1686 nop.i 999 1687} 1688;; 1689 1690{ .mfb 1691 nop.m 999 1692 fma.d.s0 f8 = POW_d, POW_poly_d, f1 1693 br.ret.sptk b0 // exit function for arguments |y*log(x)| < 2^(-11) 1694} 1695;; 1696 1697POW_POSSIBLE_UNDER: 1698// We got an answer. input was < -2^9 but > -2^10 (double) 1699// We got an answer. input was < -2^6 but > -2^7 (float) 1700// underflow is a possibility, not a certainty 1701 1702// We define an underflow when the answer with 1703// ftz set 1704// is zero (tiny numbers become zero) 1705// Notice (from below) that if we have an unlimited exponent range, 1706// then there is an extra machine number E between the largest denormal and 1707// the smallest normal. 1708// So if with unbounded exponent we round to E or below, then we are 1709// tiny and underflow has occurred. 1710// But notice that you can be in a situation where we are tiny, namely 1711// rounded to E, but when the exponent is bounded we round to smallest 1712// normal. So the answer can be the smallest normal with underflow. 1713// E 1714// -----+--------------------+--------------------+----- 1715// | | | 1716// 1.1...10 2^-3fff 1.1...11 2^-3fff 1.0...00 2^-3ffe 1717// 0.1...11 2^-3ffe (biased, 1) 1718// largest dn smallest normal 1719 1720// Put in s2 (td set, ftz set) 1721{ .mfi 1722 nop.m 999 1723 fsetc.s2 0x7F,0x41 1724 nop.i 999 1725} 1726;; 1727 1728{ .mfi 1729 nop.m 999 1730 fma.d.s2 POW_ftz_urm_f8 = POW_A, POW_q, POW_A 1731 nop.i 999 1732} 1733;; 1734 1735// Return s2 to default 1736{ .mfi 1737 nop.m 999 1738 fsetc.s2 0x7F,0x40 1739 nop.i 999 1740} 1741;; 1742 1743// p7 = TRUE ==> yes, we have an underflow 1744{ .mfi 1745 nop.m 999 1746 fcmp.eq.s1 p7, p0 = POW_ftz_urm_f8, f0 1747 nop.i 999 1748} 1749;; 1750 1751{ .mbb 1752(p7) mov pow_GR_tag = 25 1753(p7) br.cond.spnt __libm_error_region // Branch if underflow 1754 br.ret.sptk b0 // Exit if did not underflow 1755} 1756;; 1757 1758POW_X_DENORM: 1759// Here if x unorm. Use the NORM_X for getf instructions, and then back 1760// to normal path 1761{ .mfi 1762 getf.exp pow_GR_signexp_X = POW_NORM_X 1763 nop.f 999 1764 nop.i 999 1765} 1766;; 1767 1768{ .mmi 1769 getf.sig pow_GR_sig_X = POW_NORM_X 1770;; 1771 and pow_GR_exp_X = pow_GR_signexp_X, pow_GR_17ones 1772 nop.i 999 1773} 1774;; 1775 1776{ .mib 1777 sub pow_GR_true_exp_X = pow_GR_exp_X, pow_GR_16ones 1778 nop.i 999 1779 br.cond.sptk POW_COMMON 1780} 1781;; 1782 1783POW_X_0: 1784// Here if x=0 and y not nan 1785// 1786// We have the following cases: 1787// p6 x=0 and y>0 and is an integer (may be even or odd) 1788// p7 x=0 and y>0 and is NOT an integer, return +0 1789// p8 x=0 and y>0 and so big as to always be an even integer, return +0 1790// p9 x=0 and y>0 and may not be integer 1791// p10 x=0 and y>0 and is an odd integer, return x 1792// p11 x=0 and y>0 and is an even integer, return +0 1793// p12 used in dummy fcmp to set denormal flag if y=unorm 1794// p13 x=0 and y>0 1795// p14 x=0 and y=0, branch to code for calling error handling 1796// p15 x=0 and y<0, branch to code for calling error handling 1797// 1798{ .mfi 1799 getf.sig pow_GR_sig_int_Y = POW_int_Y // Get signif of int_Y 1800 fcmp.lt.s1 p15,p13 = f9, f0 // Test for y<0 1801 and pow_GR_exp_Y = pow_GR_signexp_Y, pow_GR_17ones 1802} 1803{ .mfb 1804 cmp.ne p14,p0 = pow_GR_y_zero,r0 // Test for y=0 1805 fcvt.xf POW_float_int_Y = POW_int_Y 1806(p14) br.cond.spnt POW_X_0_Y_0 // Branch if x=0 and y=0 1807} 1808;; 1809 1810// If x=0 and y>0, test y and flag denormal 1811{ .mfb 1812(p13) cmp.gt.unc p8,p9 = pow_GR_exp_Y, pow_GR_10033 // Test y +big = even int 1813(p13) fcmp.eq.s0 p12,p0 = f9,f0 // If x=0, y>0 dummy op to flag denormal 1814(p15) br.cond.spnt POW_X_0_Y_NEG // Branch if x=0 and y<0 1815} 1816;; 1817 1818// Here if x=0 and y>0 1819{ .mfi 1820 nop.m 999 1821(p9) fcmp.eq.unc.s1 p6,p7 = POW_float_int_Y, POW_NORM_Y // Test y=int 1822 nop.i 999 1823} 1824{ .mfi 1825 nop.m 999 1826(p8) fma.d.s0 f8 = f0,f0,f0 // If x=0, y>0 and large even int, return +0 1827 nop.i 999 1828} 1829;; 1830 1831{ .mfi 1832 nop.m 999 1833(p7) fma.d.s0 f8 = f0,f0,f0 // Result +0 if x=0 and y>0 and not integer 1834(p6) tbit.nz.unc p10,p11 = pow_GR_sig_int_Y,0 // If y>0 int, test y even/odd 1835} 1836;; 1837 1838// Note if x=0, y>0 and odd integer, just return x 1839{ .mfb 1840 nop.m 999 1841(p11) fma.d.s0 f8 = f0,f0,f0 // Result +0 if x=0 and y even integer 1842 br.ret.sptk b0 // Exit if x=0 and y>0 1843} 1844;; 1845 1846POW_X_0_Y_0: 1847// When X is +-0 and Y is +-0, IEEE returns 1.0 1848// We call error support with this value 1849 1850{ .mfb 1851 mov pow_GR_tag = 26 1852 fma.d.s0 f8 = f1,f1,f0 1853 br.cond.sptk __libm_error_region 1854} 1855;; 1856 1857POW_X_0_Y_NEG: 1858// When X is +-0 and Y is negative, IEEE returns 1859// X Y answer 1860// +0 -odd int +inf 1861// -0 -odd int -inf 1862 1863// +0 !-odd int +inf 1864// -0 !-odd int +inf 1865 1866// p6 == Y is a floating point number outside the integer. 1867// Hence it is an integer and is even. 1868// return +inf 1869 1870// p7 == Y is a floating point number within the integer range. 1871// p9 == (int_Y = NORM_Y), Y is an integer, which may be odd or even. 1872// p11 odd 1873// return (sign_of_x)inf 1874// p12 even 1875// return +inf 1876// p10 == Y is not an integer 1877// return +inf 1878// 1879 1880{ .mfi 1881 nop.m 999 1882 nop.f 999 1883 cmp.gt p6,p7 = pow_GR_exp_Y, pow_GR_10033 1884} 1885;; 1886 1887{ .mfi 1888 mov pow_GR_tag = 27 1889(p7) fcmp.eq.unc.s1 p9,p10 = POW_float_int_Y, POW_NORM_Y 1890 nop.i 999 1891} 1892;; 1893 1894{ .mfb 1895 nop.m 999 1896(p6) frcpa.s0 f8,p13 = f1, f0 1897(p6) br.cond.sptk __libm_error_region // x=0, y<0, y large neg int 1898} 1899;; 1900 1901{ .mfb 1902 nop.m 999 1903(p10) frcpa.s0 f8,p13 = f1, f0 1904(p10) br.cond.sptk __libm_error_region // x=0, y<0, y not int 1905} 1906;; 1907 1908// x=0, y<0, y an int 1909{ .mib 1910 nop.m 999 1911(p9) tbit.nz.unc p11,p12 = pow_GR_sig_int_Y,0 1912 nop.b 999 1913} 1914;; 1915 1916{ .mfi 1917 nop.m 999 1918(p12) frcpa.s0 f8,p13 = f1,f0 1919 nop.i 999 1920} 1921;; 1922 1923{ .mfb 1924 nop.m 999 1925(p11) frcpa.s0 f8,p13 = f1,f8 1926 br.cond.sptk __libm_error_region 1927} 1928;; 1929 1930 1931POW_Y_0: 1932// Here for y zero, x anything but zero and nan 1933// Set flag if x denormal 1934// Result is +1.0 1935{ .mfi 1936 nop.m 999 1937 fcmp.eq.s0 p6,p0 = f8,f0 // Sets flag if x denormal 1938 nop.i 999 1939} 1940{ .mfb 1941 nop.m 999 1942 fma.d.s0 f8 = f1,f1,f0 1943 br.ret.sptk b0 1944} 1945;; 1946 1947 1948POW_X_INF: 1949// Here when X is +-inf 1950 1951// X +inf Y +inf +inf 1952// X -inf Y +inf +inf 1953 1954// X +inf Y >0 +inf 1955// X -inf Y >0, !odd integer +inf <== (-inf)^0.5 = +inf !! 1956// X -inf Y >0, odd integer -inf 1957 1958// X +inf Y -inf +0 1959// X -inf Y -inf +0 1960 1961// X +inf Y <0 +0 1962// X -inf Y <0, !odd integer +0 1963// X -inf Y <0, odd integer -0 1964 1965// X + inf Y=+0 +1 1966// X + inf Y=-0 +1 1967// X - inf Y=+0 +1 1968// X - inf Y=-0 +1 1969 1970// p13 == Y negative 1971// p14 == Y positive 1972 1973// p6 == Y is a floating point number outside the integer. 1974// Hence it is an integer and is even. 1975// p13 == (Y negative) 1976// return +inf 1977// p14 == (Y positive) 1978// return +0 1979 1980// p7 == Y is a floating point number within the integer range. 1981// p9 == (int_Y = NORM_Y), Y is an integer, which may be odd or even. 1982// p11 odd 1983// p13 == (Y negative) 1984// return (sign_of_x)inf 1985// p14 == (Y positive) 1986// return (sign_of_x)0 1987// pxx even 1988// p13 == (Y negative) 1989// return +inf 1990// p14 == (Y positive) 1991// return +0 1992 1993// pxx == Y is not an integer 1994// p13 == (Y negative) 1995// return +inf 1996// p14 == (Y positive) 1997// return +0 1998// 1999 2000// If x=inf, test y and flag denormal 2001{ .mfi 2002 nop.m 999 2003 fcmp.eq.s0 p10,p11 = f9,f0 2004 nop.i 999 2005} 2006;; 2007 2008{ .mfi 2009 nop.m 999 2010 fcmp.lt.s0 p13,p14 = POW_NORM_Y,f0 2011 cmp.gt p6,p7 = pow_GR_exp_Y, pow_GR_10033 2012} 2013{ .mfi 2014 nop.m 999 2015 fclass.m p12,p0 = f9, 0x23 //@inf 2016 nop.i 999 2017} 2018;; 2019 2020{ .mfi 2021 nop.m 999 2022 fclass.m p15,p0 = f9, 0x07 //@zero 2023 nop.i 999 2024} 2025;; 2026 2027{ .mfb 2028 nop.m 999 2029(p15) fmerge.s f8 = f1,f1 // Return +1.0 if x=inf, y=0 2030(p15) br.ret.spnt b0 // Exit if x=inf, y=0 2031} 2032;; 2033 2034{ .mfi 2035 nop.m 999 2036(p14) frcpa.s1 f8,p10 = f1,f0 // If x=inf, y>0, assume result +inf 2037 nop.i 999 2038} 2039{ .mfb 2040 nop.m 999 2041(p13) fma.d.s0 f8 = f0,f0,f0 // If x=inf, y<0, assume result +0.0 2042(p12) br.ret.spnt b0 // Exit if x=inf, y=inf 2043} 2044;; 2045 2046// Here if x=inf, and 0 < |y| < inf. Need to correct results if y odd integer. 2047{ .mfi 2048 nop.m 999 2049(p7) fcmp.eq.unc.s1 p9,p0 = POW_float_int_Y, POW_NORM_Y // Is y integer? 2050 nop.i 999 2051} 2052;; 2053 2054{ .mfi 2055 nop.m 999 2056 nop.f 999 2057(p9) tbit.nz.unc p11,p0 = pow_GR_sig_int_Y,0 // Test for y odd integer 2058} 2059;; 2060 2061{ .mfb 2062 nop.m 999 2063(p11) fmerge.s f8 = POW_NORM_X,f8 // If y odd integer use sign of x 2064 br.ret.sptk b0 // Exit for x=inf, 0 < |y| < inf 2065} 2066;; 2067 2068 2069POW_X_NEG_Y_NONINT: 2070// When X is negative and Y is a non-integer, IEEE 2071// returns a qnan indefinite. 2072// We call error support with this value 2073 2074{ .mfb 2075 mov pow_GR_tag = 28 2076 frcpa.s0 f8,p6 = f0,f0 2077 br.cond.sptk __libm_error_region 2078} 2079;; 2080 2081POW_X_NAN: 2082// Here if x=nan, y not nan 2083{ .mfi 2084 nop.m 999 2085 fclass.m p9,p13 = f9, 0x07 // Test y=zero 2086 nop.i 999 2087} 2088;; 2089 2090{ .mfb 2091 nop.m 999 2092(p13) fma.d.s0 f8 = f8,f1,f0 2093(p13) br.ret.sptk b0 // Exit if x nan, y anything but zero or nan 2094} 2095;; 2096 2097POW_X_NAN_Y_0: 2098// When X is a NAN and Y is zero, IEEE returns 1. 2099// We call error support with this value. 2100{ .mfi 2101 nop.m 999 2102 fcmp.eq.s0 p6,p0 = f8,f0 // Dummy op to set invalid on snan 2103 nop.i 999 2104} 2105{ .mfb 2106 mov pow_GR_tag = 29 2107 fma.d.s0 f8 = f0,f0,f1 2108 br.cond.sptk __libm_error_region 2109} 2110;; 2111 2112 2113POW_OVER_UNDER_X_NOT_INF: 2114 2115// p8 is TRUE for overflow 2116// p9 is TRUE for underflow 2117 2118// if y is infinity, we should not over/underflow 2119 2120{ .mfi 2121 nop.m 999 2122 fcmp.eq.s1 p14, p13 = POW_xsq,f1 // Test |x|=1 2123 cmp.eq p8,p9 = pow_GR_sign_Y_Gpr, r0 2124} 2125;; 2126 2127{ .mfi 2128 nop.m 999 2129(p14) fclass.m.unc p15, p0 = f9, 0x23 // If |x|=1, test y=inf 2130 nop.i 999 2131} 2132{ .mfi 2133 nop.m 999 2134(p13) fclass.m.unc p11,p0 = f9, 0x23 // If |x| not 1, test y=inf 2135 nop.i 999 2136} 2137;; 2138 2139// p15 = TRUE if |x|=1, y=inf, return +1 2140{ .mfb 2141 nop.m 999 2142(p15) fma.d.s0 f8 = f1,f1,f0 // If |x|=1, y=inf, result +1 2143(p15) br.ret.spnt b0 // Exit if |x|=1, y=inf 2144} 2145;; 2146 2147.pred.rel "mutex",p8,p9 2148{ .mfb 2149(p8) setf.exp f8 = pow_GR_17ones // If exp(+big), result inf 2150(p9) fmerge.s f8 = f0,f0 // If exp(-big), result 0 2151(p11) br.ret.sptk b0 // Exit if |x| not 1, y=inf 2152} 2153;; 2154 2155{ .mfb 2156 nop.m 999 2157 nop.f 999 2158 br.cond.sptk POW_OVER_UNDER_ERROR // Branch if y not inf 2159} 2160;; 2161 2162 2163POW_Y_NAN: 2164// Here if y=nan, x anything 2165// If x = +1 then result is +1, else result is quiet Y 2166{ .mfi 2167 nop.m 999 2168 fcmp.eq.s1 p10,p9 = POW_NORM_X, f1 2169 nop.i 999 2170} 2171;; 2172 2173{ .mfi 2174 nop.m 999 2175(p10) fcmp.eq.s0 p6,p0 = f9,f1 // Set invalid, even if x=+1 2176 nop.i 999 2177} 2178;; 2179 2180{ .mfi 2181 nop.m 999 2182(p10) fma.d.s0 f8 = f1,f1,f0 2183 nop.i 999 2184} 2185{ .mfb 2186 nop.m 999 2187(p9) fma.d.s0 f8 = f9,f8,f0 2188 br.ret.sptk b0 // Exit y=nan 2189} 2190;; 2191 2192 2193POW_OVER_UNDER_ERROR: 2194// Here if we have overflow or underflow. 2195// Enter with p12 true if x negative and y odd int to force -0 or -inf 2196 2197{ .mfi 2198 sub pow_GR_17ones_m1 = pow_GR_17ones, r0, 1 2199 nop.f 999 2200 mov pow_GR_one = 0x1 2201} 2202;; 2203 2204// overflow, force inf with O flag 2205{ .mmb 2206(p8) mov pow_GR_tag = 24 2207(p8) setf.exp POW_tmp = pow_GR_17ones_m1 2208 nop.b 999 2209} 2210;; 2211 2212// underflow, force zero with I, U flags 2213{ .mmi 2214(p9) mov pow_GR_tag = 25 2215(p9) setf.exp POW_tmp = pow_GR_one 2216 nop.i 999 2217} 2218;; 2219 2220{ .mfi 2221 nop.m 999 2222 fma.d.s0 f8 = POW_tmp, POW_tmp, f0 2223 nop.i 999 2224} 2225;; 2226 2227// p12 x is negative and y is an odd integer, change sign of result 2228{ .mfi 2229 nop.m 999 2230(p12) fnma.d.s0 f8 = POW_tmp, POW_tmp, f0 2231 nop.i 999 2232} 2233;; 2234 2235WEAK_LIBM_END(pow) 2236libm_alias_double_other (__pow, pow) 2237#ifdef SHARED 2238.symver pow,pow@@GLIBC_2.29 2239.weak __pow_compat 2240.set __pow_compat,__pow 2241.symver __pow_compat,pow@GLIBC_2.2 2242#endif 2243 2244 2245LOCAL_LIBM_ENTRY(__libm_error_region) 2246 2247.prologue 2248{ .mfi 2249 add GR_Parameter_Y=-32,sp // Parameter 2 value 2250 nop.f 0 2251.save ar.pfs,GR_SAVE_PFS 2252 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs 2253} 2254{ .mfi 2255.fframe 64 2256 add sp=-64,sp // Create new stack 2257 nop.f 0 2258 mov GR_SAVE_GP=gp // Save gp 2259};; 2260 2261{ .mmi 2262 stfd [GR_Parameter_Y] = POW_NORM_Y,16 // STORE Parameter 2 on stack 2263 add GR_Parameter_X = 16,sp // Parameter 1 address 2264.save b0, GR_SAVE_B0 2265 mov GR_SAVE_B0=b0 // Save b0 2266};; 2267 2268.body 2269{ .mib 2270 stfd [GR_Parameter_X] = POW_NORM_X // STORE Parameter 1 on stack 2271 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address 2272 nop.b 0 2273} 2274{ .mib 2275 stfd [GR_Parameter_Y] = f8 // STORE Parameter 3 on stack 2276 add GR_Parameter_Y = -16,GR_Parameter_Y 2277 br.call.sptk b0=__libm_error_support# // Call error handling function 2278};; 2279 2280{ .mmi 2281 add GR_Parameter_RESULT = 48,sp 2282 nop.m 0 2283 nop.i 0 2284};; 2285 2286{ .mmi 2287 ldfd f8 = [GR_Parameter_RESULT] // Get return result off stack 2288.restore sp 2289 add sp = 64,sp // Restore stack pointer 2290 mov b0 = GR_SAVE_B0 // Restore return address 2291};; 2292 2293{ .mib 2294 mov gp = GR_SAVE_GP // Restore gp 2295 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs 2296 br.ret.sptk b0 // Return 2297};; 2298 2299LOCAL_LIBM_END(__libm_error_region) 2300 2301.type __libm_error_support#,@function 2302.global __libm_error_support# 2303