1.file "powf.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^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 tabluate 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// 175 176// Table Generation 177//============================================================== 178 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). 242 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 332 333 334// Special Cases 335//============================================================== 336 337// double float 338// overflow error 24 30 339 340// underflow error 25 31 341 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 347 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 355 356// X zero Y positve 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 365 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 371 372// X - Y not integer QNAN error 28 34 invalid 373 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 378 379 380// X +inf Y >0 +inf 381// X -inf Y >0, !odd integer +inf 382// X -inf Y >0, odd integer -inf 383 384// X +inf Y <0 +0 385// X -inf Y <0, !odd integer +0 386// X -inf Y <0, odd integer -0 387 388// X +inf Y =0 +1 389// X -inf Y =0 +1 390 391// |X|<1 Y +inf +0 392// |X|<1 Y -inf +inf 393// |X|>1 Y +inf +inf 394// |X|>1 Y -inf +0 395 396// X any Y =0 +1 397 398// Assembly macros 399//============================================================== 400 401// integer registers used 402 403pow_GR_exp_half = r10 404pow_GR_signexp_Xm1 = r11 405pow_GR_tmp = r11 406 407pow_GR_signexp_X = r14 408pow_GR_17ones = r15 409pow_GR_Fpsr = r15 410pow_AD_P = r16 411pow_GR_rcs0_mask = r16 412pow_GR_exp_2tom8 = r17 413pow_GR_rcs0 = r17 414pow_GR_sig_X = r18 415pow_GR_10033 = r19 416pow_GR_16ones = r20 417 418pow_AD_Tt = r21 419pow_GR_exp_X = r22 420pow_AD_Q = r23 421pow_GR_true_exp_X = r24 422pow_GR_y_zero = r25 423 424pow_GR_exp_Y = r26 425pow_AD_tbl1 = r27 426pow_AD_tbl2 = r28 427pow_GR_offset = r29 428pow_GR_exp_Xm1 = r30 429pow_GR_xneg_yodd = r31 430 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_big_pos = r55 451 452pow_GR_big_neg = r56 453 454GR_SAVE_B0 = r50 455GR_SAVE_GP = r51 456GR_SAVE_PFS = r52 457 458GR_Parameter_X = r53 459GR_Parameter_Y = r54 460GR_Parameter_RESULT = r55 461pow_GR_tag = r56 462 463 464// floating point registers used 465 466POW_B = f32 467POW_NORM_X = f33 468POW_Xm1 = f34 469POW_r1 = f34 470 471POW_NORM_Y = f37 472POW_Q2 = f38 473POW_eps = f39 474POW_P2 = f40 475 476POW_P0 = f42 477POW_log2_lo = f43 478POW_r = f44 479POW_Q0_half = f45 480 481POW_tmp = f47 482POW_log2_hi = f48 483POW_Q1 = f49 484POW_P1 = f50 485 486POW_log2_by_128_hi = f51 487POW_inv_log2_by_128 = f52 488POW_rsq = f53 489POW_Yrcub = f54 490POW_log2_by_128_lo = f55 491 492POW_xsq = f57 493POW_v2 = f59 494POW_T = f60 495 496POW_RSHF = f62 497POW_v210 = f63 498POW_twoV = f65 499 500POW_U = f66 501POW_G = f67 502POW_delta = f68 503POW_V = f70 504 505POW_p = f71 506POW_Z = f72 507POW_e3 = f73 508POW_Z2 = f75 509 510POW_W1 = f77 511POW_Z3 = f80 512 513POW_Z3sq = f85 514 515POW_Nfloat = f87 516POW_f3 = f89 517POW_q = f90 518 519POW_T1 = f96 520POW_T2 = f97 521POW_2M = f98 522POW_s = f99 523POW_f12 = f100 524 525POW_ssq = f101 526POW_T1T2 = f102 527POW_1ps = f103 528POW_A = f104 529POW_es = f105 530 531POW_Xp1 = f106 532POW_int_K = f107 533POW_K = f108 534POW_f123 = f109 535POW_Gpr = f110 536 537POW_Y_Gpr = f111 538POW_int_Y = f112 539POW_2Mqp1 = f113 540 541POW_float_int_Y = f116 542POW_ftz_urm_f8 = f117 543POW_wre_urm_f8 = f118 544POW_big_neg = f119 545POW_big_pos = f120 546 547// Data tables 548//============================================================== 549 550RODATA 551 552.align 16 553 554LOCAL_OBJECT_START(pow_table_P) 555data8 0x80000000000018E5, 0x0000BFFD // P_1 556data8 0xb8aa3b295c17f0bc, 0x00004006 // inv_ln2_by_128 557// 558// 559data8 0x3FA5555555554A9E // Q_2 560data8 0x0000000000000000 // Pad 561data8 0x3FC5555555554733 // Q_1 562data8 0x43e8000000000000 // Right shift constant for exp 563data8 0xc9e3b39803f2f6af, 0x00003fb7 // ln2_by_128_lo 564LOCAL_OBJECT_END(pow_table_P) 565 566LOCAL_OBJECT_START(pow_table_Q) 567data8 0xCCCCCCCC4ED2BA7F, 0x00003FFC // P_2 568data8 0xAAAAAAAAAAAAB505, 0x00003FFD // P_0 569data8 0x3fe62e42fefa39e8, 0x3cccd5e4f1d9cc02 // log2 hi lo = +6.93147e-001 570data8 0xb17217f7d1cf79ab, 0x00003ff7 // ln2_by_128_hi 571LOCAL_OBJECT_END(pow_table_Q) 572 573 574LOCAL_OBJECT_START(pow_Tt) 575data8 0x3f60040155d58800 // log(1/frcpa(1+0/256))= +1.95503e-003 576data8 0x3f78121214586a00 // log(1/frcpa(1+1/256))= +5.87661e-003 577data8 0x3f841929f9683200 // log(1/frcpa(1+2/256))= +9.81362e-003 578data8 0x3f8c317384c75f00 // log(1/frcpa(1+3/256))= +1.37662e-002 579data8 0x3f91a6b91ac73380 // log(1/frcpa(1+4/256))= +1.72376e-002 580data8 0x3f95ba9a5d9ac000 // log(1/frcpa(1+5/256))= +2.12196e-002 581data8 0x3f99d2a807432580 // log(1/frcpa(1+6/256))= +2.52177e-002 582data8 0x3f9d6b2725979800 // log(1/frcpa(1+7/256))= +2.87291e-002 583data8 0x3fa0c58fa19dfa80 // log(1/frcpa(1+8/256))= +3.27573e-002 584data8 0x3fa2954c78cbce00 // log(1/frcpa(1+9/256))= +3.62953e-002 585data8 0x3fa4a94d2da96c40 // log(1/frcpa(1+10/256))= +4.03542e-002 586data8 0x3fa67c94f2d4bb40 // log(1/frcpa(1+11/256))= +4.39192e-002 587data8 0x3fa85188b630f040 // log(1/frcpa(1+12/256))= +4.74971e-002 588data8 0x3faa6b8abe73af40 // log(1/frcpa(1+13/256))= +5.16017e-002 589data8 0x3fac441e06f72a80 // log(1/frcpa(1+14/256))= +5.52072e-002 590data8 0x3fae1e6713606d00 // log(1/frcpa(1+15/256))= +5.88257e-002 591data8 0x3faffa6911ab9300 // log(1/frcpa(1+16/256))= +6.24574e-002 592data8 0x3fb0ec139c5da600 // log(1/frcpa(1+17/256))= +6.61022e-002 593data8 0x3fb1dbd2643d1900 // log(1/frcpa(1+18/256))= +6.97605e-002 594data8 0x3fb2cc7284fe5f00 // log(1/frcpa(1+19/256))= +7.34321e-002 595data8 0x3fb3bdf5a7d1ee60 // log(1/frcpa(1+20/256))= +7.71173e-002 596data8 0x3fb4b05d7aa012e0 // log(1/frcpa(1+21/256))= +8.08161e-002 597data8 0x3fb580db7ceb5700 // log(1/frcpa(1+22/256))= +8.39975e-002 598data8 0x3fb674f089365a60 // log(1/frcpa(1+23/256))= +8.77219e-002 599data8 0x3fb769ef2c6b5680 // log(1/frcpa(1+24/256))= +9.14602e-002 600data8 0x3fb85fd927506a40 // log(1/frcpa(1+25/256))= +9.52125e-002 601data8 0x3fb9335e5d594980 // log(1/frcpa(1+26/256))= +9.84401e-002 602data8 0x3fba2b0220c8e5e0 // log(1/frcpa(1+27/256))= +1.02219e-001 603data8 0x3fbb0004ac1a86a0 // log(1/frcpa(1+28/256))= +1.05469e-001 604data8 0x3fbbf968769fca00 // log(1/frcpa(1+29/256))= +1.09274e-001 605data8 0x3fbccfedbfee13a0 // log(1/frcpa(1+30/256))= +1.12548e-001 606data8 0x3fbda727638446a0 // log(1/frcpa(1+31/256))= +1.15832e-001 607data8 0x3fbea3257fe10f60 // log(1/frcpa(1+32/256))= +1.19677e-001 608data8 0x3fbf7be9fedbfde0 // log(1/frcpa(1+33/256))= +1.22985e-001 609data8 0x3fc02ab352ff25f0 // log(1/frcpa(1+34/256))= +1.26303e-001 610data8 0x3fc097ce579d2040 // log(1/frcpa(1+35/256))= +1.29633e-001 611data8 0x3fc1178e8227e470 // log(1/frcpa(1+36/256))= +1.33531e-001 612data8 0x3fc185747dbecf30 // log(1/frcpa(1+37/256))= +1.36885e-001 613data8 0x3fc1f3b925f25d40 // log(1/frcpa(1+38/256))= +1.40250e-001 614data8 0x3fc2625d1e6ddf50 // log(1/frcpa(1+39/256))= +1.43627e-001 615data8 0x3fc2d1610c868130 // log(1/frcpa(1+40/256))= +1.47015e-001 616data8 0x3fc340c597411420 // log(1/frcpa(1+41/256))= +1.50414e-001 617data8 0x3fc3b08b6757f2a0 // log(1/frcpa(1+42/256))= +1.53825e-001 618data8 0x3fc40dfb08378000 // log(1/frcpa(1+43/256))= +1.56677e-001 619data8 0x3fc47e74e8ca5f70 // log(1/frcpa(1+44/256))= +1.60109e-001 620data8 0x3fc4ef51f6466de0 // log(1/frcpa(1+45/256))= +1.63553e-001 621data8 0x3fc56092e02ba510 // log(1/frcpa(1+46/256))= +1.67010e-001 622data8 0x3fc5d23857cd74d0 // log(1/frcpa(1+47/256))= +1.70478e-001 623data8 0x3fc6313a37335d70 // log(1/frcpa(1+48/256))= +1.73377e-001 624data8 0x3fc6a399dabbd380 // log(1/frcpa(1+49/256))= +1.76868e-001 625data8 0x3fc70337dd3ce410 // log(1/frcpa(1+50/256))= +1.79786e-001 626data8 0x3fc77654128f6120 // log(1/frcpa(1+51/256))= +1.83299e-001 627data8 0x3fc7e9d82a0b0220 // log(1/frcpa(1+52/256))= +1.86824e-001 628data8 0x3fc84a6b759f5120 // log(1/frcpa(1+53/256))= +1.89771e-001 629data8 0x3fc8ab47d5f5a300 // log(1/frcpa(1+54/256))= +1.92727e-001 630data8 0x3fc91fe490965810 // log(1/frcpa(1+55/256))= +1.96286e-001 631data8 0x3fc981634011aa70 // log(1/frcpa(1+56/256))= +1.99261e-001 632data8 0x3fc9f6c407089660 // log(1/frcpa(1+57/256))= +2.02843e-001 633data8 0x3fca58e729348f40 // log(1/frcpa(1+58/256))= +2.05838e-001 634data8 0x3fcabb55c31693a0 // log(1/frcpa(1+59/256))= +2.08842e-001 635data8 0x3fcb1e104919efd0 // log(1/frcpa(1+60/256))= +2.11855e-001 636data8 0x3fcb94ee93e367c0 // log(1/frcpa(1+61/256))= +2.15483e-001 637data8 0x3fcbf851c0675550 // log(1/frcpa(1+62/256))= +2.18516e-001 638data8 0x3fcc5c0254bf23a0 // log(1/frcpa(1+63/256))= +2.21558e-001 639data8 0x3fccc000c9db3c50 // log(1/frcpa(1+64/256))= +2.24609e-001 640data8 0x3fcd244d99c85670 // log(1/frcpa(1+65/256))= +2.27670e-001 641data8 0x3fcd88e93fb2f450 // log(1/frcpa(1+66/256))= +2.30741e-001 642data8 0x3fcdedd437eaef00 // log(1/frcpa(1+67/256))= +2.33820e-001 643data8 0x3fce530effe71010 // log(1/frcpa(1+68/256))= +2.36910e-001 644data8 0x3fceb89a1648b970 // log(1/frcpa(1+69/256))= +2.40009e-001 645data8 0x3fcf1e75fadf9bd0 // log(1/frcpa(1+70/256))= +2.43117e-001 646data8 0x3fcf84a32ead7c30 // log(1/frcpa(1+71/256))= +2.46235e-001 647data8 0x3fcfeb2233ea07c0 // log(1/frcpa(1+72/256))= +2.49363e-001 648data8 0x3fd028f9c7035c18 // log(1/frcpa(1+73/256))= +2.52501e-001 649data8 0x3fd05c8be0d96358 // log(1/frcpa(1+74/256))= +2.55649e-001 650data8 0x3fd085eb8f8ae790 // log(1/frcpa(1+75/256))= +2.58174e-001 651data8 0x3fd0b9c8e32d1910 // log(1/frcpa(1+76/256))= +2.61339e-001 652data8 0x3fd0edd060b78080 // log(1/frcpa(1+77/256))= +2.64515e-001 653data8 0x3fd122024cf00638 // log(1/frcpa(1+78/256))= +2.67701e-001 654data8 0x3fd14be2927aecd0 // log(1/frcpa(1+79/256))= +2.70257e-001 655data8 0x3fd180618ef18ad8 // log(1/frcpa(1+80/256))= +2.73461e-001 656data8 0x3fd1b50bbe2fc638 // log(1/frcpa(1+81/256))= +2.76675e-001 657data8 0x3fd1df4cc7cf2428 // log(1/frcpa(1+82/256))= +2.79254e-001 658data8 0x3fd214456d0eb8d0 // log(1/frcpa(1+83/256))= +2.82487e-001 659data8 0x3fd23ec5991eba48 // log(1/frcpa(1+84/256))= +2.85081e-001 660data8 0x3fd2740d9f870af8 // log(1/frcpa(1+85/256))= +2.88333e-001 661data8 0x3fd29ecdabcdfa00 // log(1/frcpa(1+86/256))= +2.90943e-001 662data8 0x3fd2d46602adcce8 // log(1/frcpa(1+87/256))= +2.94214e-001 663data8 0x3fd2ff66b04ea9d0 // log(1/frcpa(1+88/256))= +2.96838e-001 664data8 0x3fd335504b355a30 // log(1/frcpa(1+89/256))= +3.00129e-001 665data8 0x3fd360925ec44f58 // log(1/frcpa(1+90/256))= +3.02769e-001 666data8 0x3fd38bf1c3337e70 // log(1/frcpa(1+91/256))= +3.05417e-001 667data8 0x3fd3c25277333180 // log(1/frcpa(1+92/256))= +3.08735e-001 668data8 0x3fd3edf463c16838 // log(1/frcpa(1+93/256))= +3.11399e-001 669data8 0x3fd419b423d5e8c0 // log(1/frcpa(1+94/256))= +3.14069e-001 670data8 0x3fd44591e0539f48 // log(1/frcpa(1+95/256))= +3.16746e-001 671data8 0x3fd47c9175b6f0a8 // log(1/frcpa(1+96/256))= +3.20103e-001 672data8 0x3fd4a8b341552b08 // log(1/frcpa(1+97/256))= +3.22797e-001 673data8 0x3fd4d4f390890198 // log(1/frcpa(1+98/256))= +3.25498e-001 674data8 0x3fd501528da1f960 // log(1/frcpa(1+99/256))= +3.28206e-001 675data8 0x3fd52dd06347d4f0 // log(1/frcpa(1+100/256))= +3.30921e-001 676data8 0x3fd55a6d3c7b8a88 // log(1/frcpa(1+101/256))= +3.33644e-001 677data8 0x3fd5925d2b112a58 // log(1/frcpa(1+102/256))= +3.37058e-001 678data8 0x3fd5bf406b543db0 // log(1/frcpa(1+103/256))= +3.39798e-001 679data8 0x3fd5ec433d5c35a8 // log(1/frcpa(1+104/256))= +3.42545e-001 680data8 0x3fd61965cdb02c18 // log(1/frcpa(1+105/256))= +3.45300e-001 681data8 0x3fd646a84935b2a0 // log(1/frcpa(1+106/256))= +3.48063e-001 682data8 0x3fd6740add31de90 // log(1/frcpa(1+107/256))= +3.50833e-001 683data8 0x3fd6a18db74a58c0 // log(1/frcpa(1+108/256))= +3.53610e-001 684data8 0x3fd6cf31058670e8 // log(1/frcpa(1+109/256))= +3.56396e-001 685data8 0x3fd6f180e852f0b8 // log(1/frcpa(1+110/256))= +3.58490e-001 686data8 0x3fd71f5d71b894e8 // log(1/frcpa(1+111/256))= +3.61289e-001 687data8 0x3fd74d5aefd66d58 // log(1/frcpa(1+112/256))= +3.64096e-001 688data8 0x3fd77b79922bd378 // log(1/frcpa(1+113/256))= +3.66911e-001 689data8 0x3fd7a9b9889f19e0 // log(1/frcpa(1+114/256))= +3.69734e-001 690data8 0x3fd7d81b037eb6a0 // log(1/frcpa(1+115/256))= +3.72565e-001 691data8 0x3fd8069e33827230 // log(1/frcpa(1+116/256))= +3.75404e-001 692data8 0x3fd82996d3ef8bc8 // log(1/frcpa(1+117/256))= +3.77538e-001 693data8 0x3fd85855776dcbf8 // log(1/frcpa(1+118/256))= +3.80391e-001 694data8 0x3fd8873658327cc8 // log(1/frcpa(1+119/256))= +3.83253e-001 695data8 0x3fd8aa75973ab8c8 // log(1/frcpa(1+120/256))= +3.85404e-001 696data8 0x3fd8d992dc8824e0 // log(1/frcpa(1+121/256))= +3.88280e-001 697data8 0x3fd908d2ea7d9510 // log(1/frcpa(1+122/256))= +3.91164e-001 698data8 0x3fd92c59e79c0e50 // log(1/frcpa(1+123/256))= +3.93332e-001 699data8 0x3fd95bd750ee3ed0 // log(1/frcpa(1+124/256))= +3.96231e-001 700data8 0x3fd98b7811a3ee58 // log(1/frcpa(1+125/256))= +3.99138e-001 701data8 0x3fd9af47f33d4068 // log(1/frcpa(1+126/256))= +4.01323e-001 702data8 0x3fd9df270c1914a0 // log(1/frcpa(1+127/256))= +4.04245e-001 703data8 0x3fda0325ed14fda0 // log(1/frcpa(1+128/256))= +4.06442e-001 704data8 0x3fda33440224fa78 // log(1/frcpa(1+129/256))= +4.09379e-001 705data8 0x3fda57725e80c380 // log(1/frcpa(1+130/256))= +4.11587e-001 706data8 0x3fda87d0165dd198 // log(1/frcpa(1+131/256))= +4.14539e-001 707data8 0x3fdaac2e6c03f890 // log(1/frcpa(1+132/256))= +4.16759e-001 708data8 0x3fdadccc6fdf6a80 // log(1/frcpa(1+133/256))= +4.19726e-001 709data8 0x3fdb015b3eb1e790 // log(1/frcpa(1+134/256))= +4.21958e-001 710data8 0x3fdb323a3a635948 // log(1/frcpa(1+135/256))= +4.24941e-001 711data8 0x3fdb56fa04462908 // log(1/frcpa(1+136/256))= +4.27184e-001 712data8 0x3fdb881aa659bc90 // log(1/frcpa(1+137/256))= +4.30182e-001 713data8 0x3fdbad0bef3db160 // log(1/frcpa(1+138/256))= +4.32437e-001 714data8 0x3fdbd21297781c28 // log(1/frcpa(1+139/256))= +4.34697e-001 715data8 0x3fdc039236f08818 // log(1/frcpa(1+140/256))= +4.37718e-001 716data8 0x3fdc28cb1e4d32f8 // log(1/frcpa(1+141/256))= +4.39990e-001 717data8 0x3fdc4e19b84723c0 // log(1/frcpa(1+142/256))= +4.42267e-001 718data8 0x3fdc7ff9c74554c8 // log(1/frcpa(1+143/256))= +4.45311e-001 719data8 0x3fdca57b64e9db00 // log(1/frcpa(1+144/256))= +4.47600e-001 720data8 0x3fdccb130a5ceba8 // log(1/frcpa(1+145/256))= +4.49895e-001 721data8 0x3fdcf0c0d18f3268 // log(1/frcpa(1+146/256))= +4.52194e-001 722data8 0x3fdd232075b5a200 // log(1/frcpa(1+147/256))= +4.55269e-001 723data8 0x3fdd490246defa68 // log(1/frcpa(1+148/256))= +4.57581e-001 724data8 0x3fdd6efa918d25c8 // log(1/frcpa(1+149/256))= +4.59899e-001 725data8 0x3fdd9509707ae528 // log(1/frcpa(1+150/256))= +4.62221e-001 726data8 0x3fddbb2efe92c550 // log(1/frcpa(1+151/256))= +4.64550e-001 727data8 0x3fddee2f3445e4a8 // log(1/frcpa(1+152/256))= +4.67663e-001 728data8 0x3fde148a1a2726c8 // log(1/frcpa(1+153/256))= +4.70004e-001 729data8 0x3fde3afc0a49ff38 // log(1/frcpa(1+154/256))= +4.72350e-001 730data8 0x3fde6185206d5168 // log(1/frcpa(1+155/256))= +4.74702e-001 731data8 0x3fde882578823d50 // log(1/frcpa(1+156/256))= +4.77060e-001 732data8 0x3fdeaedd2eac9908 // log(1/frcpa(1+157/256))= +4.79423e-001 733data8 0x3fded5ac5f436be0 // log(1/frcpa(1+158/256))= +4.81792e-001 734data8 0x3fdefc9326d16ab8 // log(1/frcpa(1+159/256))= +4.84166e-001 735data8 0x3fdf2391a21575f8 // log(1/frcpa(1+160/256))= +4.86546e-001 736data8 0x3fdf4aa7ee031928 // log(1/frcpa(1+161/256))= +4.88932e-001 737data8 0x3fdf71d627c30bb0 // log(1/frcpa(1+162/256))= +4.91323e-001 738data8 0x3fdf991c6cb3b378 // log(1/frcpa(1+163/256))= +4.93720e-001 739data8 0x3fdfc07ada69a908 // log(1/frcpa(1+164/256))= +4.96123e-001 740data8 0x3fdfe7f18eb03d38 // log(1/frcpa(1+165/256))= +4.98532e-001 741data8 0x3fe007c053c5002c // log(1/frcpa(1+166/256))= +5.00946e-001 742data8 0x3fe01b942198a5a0 // log(1/frcpa(1+167/256))= +5.03367e-001 743data8 0x3fe02f74400c64e8 // log(1/frcpa(1+168/256))= +5.05793e-001 744data8 0x3fe04360be7603ac // log(1/frcpa(1+169/256))= +5.08225e-001 745data8 0x3fe05759ac47fe30 // log(1/frcpa(1+170/256))= +5.10663e-001 746data8 0x3fe06b5f1911cf50 // log(1/frcpa(1+171/256))= +5.13107e-001 747data8 0x3fe078bf0533c568 // log(1/frcpa(1+172/256))= +5.14740e-001 748data8 0x3fe08cd9687e7b0c // log(1/frcpa(1+173/256))= +5.17194e-001 749data8 0x3fe0a10074cf9018 // log(1/frcpa(1+174/256))= +5.19654e-001 750data8 0x3fe0b5343a234474 // log(1/frcpa(1+175/256))= +5.22120e-001 751data8 0x3fe0c974c89431cc // log(1/frcpa(1+176/256))= +5.24592e-001 752data8 0x3fe0ddc2305b9884 // log(1/frcpa(1+177/256))= +5.27070e-001 753data8 0x3fe0eb524bafc918 // log(1/frcpa(1+178/256))= +5.28726e-001 754data8 0x3fe0ffb54213a474 // log(1/frcpa(1+179/256))= +5.31214e-001 755data8 0x3fe114253da97d9c // log(1/frcpa(1+180/256))= +5.33709e-001 756data8 0x3fe128a24f1d9afc // log(1/frcpa(1+181/256))= +5.36210e-001 757data8 0x3fe1365252bf0864 // log(1/frcpa(1+182/256))= +5.37881e-001 758data8 0x3fe14ae558b4a92c // log(1/frcpa(1+183/256))= +5.40393e-001 759data8 0x3fe15f85a19c7658 // log(1/frcpa(1+184/256))= +5.42910e-001 760data8 0x3fe16d4d38c119f8 // log(1/frcpa(1+185/256))= +5.44592e-001 761data8 0x3fe18203c20dd130 // log(1/frcpa(1+186/256))= +5.47121e-001 762data8 0x3fe196c7bc4b1f38 // log(1/frcpa(1+187/256))= +5.49656e-001 763data8 0x3fe1a4a738b7a33c // log(1/frcpa(1+188/256))= +5.51349e-001 764data8 0x3fe1b981c0c9653c // log(1/frcpa(1+189/256))= +5.53895e-001 765data8 0x3fe1ce69e8bb1068 // log(1/frcpa(1+190/256))= +5.56447e-001 766data8 0x3fe1dc619de06944 // log(1/frcpa(1+191/256))= +5.58152e-001 767data8 0x3fe1f160a2ad0da0 // log(1/frcpa(1+192/256))= +5.60715e-001 768data8 0x3fe2066d7740737c // log(1/frcpa(1+193/256))= +5.63285e-001 769data8 0x3fe2147dba47a390 // log(1/frcpa(1+194/256))= +5.65001e-001 770data8 0x3fe229a1bc5ebac0 // log(1/frcpa(1+195/256))= +5.67582e-001 771data8 0x3fe237c1841a502c // log(1/frcpa(1+196/256))= +5.69306e-001 772data8 0x3fe24cfce6f80d98 // log(1/frcpa(1+197/256))= +5.71898e-001 773data8 0x3fe25b2c55cd5760 // log(1/frcpa(1+198/256))= +5.73630e-001 774data8 0x3fe2707f4d5f7c40 // log(1/frcpa(1+199/256))= +5.76233e-001 775data8 0x3fe285e0842ca380 // log(1/frcpa(1+200/256))= +5.78842e-001 776data8 0x3fe294294708b770 // log(1/frcpa(1+201/256))= +5.80586e-001 777data8 0x3fe2a9a2670aff0c // log(1/frcpa(1+202/256))= +5.83207e-001 778data8 0x3fe2b7fb2c8d1cc0 // log(1/frcpa(1+203/256))= +5.84959e-001 779data8 0x3fe2c65a6395f5f4 // log(1/frcpa(1+204/256))= +5.86713e-001 780data8 0x3fe2dbf557b0df40 // log(1/frcpa(1+205/256))= +5.89350e-001 781data8 0x3fe2ea64c3f97654 // log(1/frcpa(1+206/256))= +5.91113e-001 782data8 0x3fe3001823684d70 // log(1/frcpa(1+207/256))= +5.93762e-001 783data8 0x3fe30e97e9a8b5cc // log(1/frcpa(1+208/256))= +5.95531e-001 784data8 0x3fe32463ebdd34e8 // log(1/frcpa(1+209/256))= +5.98192e-001 785data8 0x3fe332f4314ad794 // log(1/frcpa(1+210/256))= +5.99970e-001 786data8 0x3fe348d90e7464cc // log(1/frcpa(1+211/256))= +6.02643e-001 787data8 0x3fe35779f8c43d6c // log(1/frcpa(1+212/256))= +6.04428e-001 788data8 0x3fe36621961a6a98 // log(1/frcpa(1+213/256))= +6.06217e-001 789data8 0x3fe37c299f3c3668 // log(1/frcpa(1+214/256))= +6.08907e-001 790data8 0x3fe38ae2171976e4 // log(1/frcpa(1+215/256))= +6.10704e-001 791data8 0x3fe399a157a603e4 // log(1/frcpa(1+216/256))= +6.12504e-001 792data8 0x3fe3afccfe77b9d0 // log(1/frcpa(1+217/256))= +6.15210e-001 793data8 0x3fe3be9d503533b4 // log(1/frcpa(1+218/256))= +6.17018e-001 794data8 0x3fe3cd7480b4a8a0 // log(1/frcpa(1+219/256))= +6.18830e-001 795data8 0x3fe3e3c43918f76c // log(1/frcpa(1+220/256))= +6.21554e-001 796data8 0x3fe3f2acb27ed6c4 // log(1/frcpa(1+221/256))= +6.23373e-001 797data8 0x3fe4019c2125ca90 // log(1/frcpa(1+222/256))= +6.25197e-001 798data8 0x3fe4181061389720 // log(1/frcpa(1+223/256))= +6.27937e-001 799data8 0x3fe42711518df544 // log(1/frcpa(1+224/256))= +6.29769e-001 800data8 0x3fe436194e12b6bc // log(1/frcpa(1+225/256))= +6.31604e-001 801data8 0x3fe445285d68ea68 // log(1/frcpa(1+226/256))= +6.33442e-001 802data8 0x3fe45bcc464c8938 // log(1/frcpa(1+227/256))= +6.36206e-001 803data8 0x3fe46aed21f117fc // log(1/frcpa(1+228/256))= +6.38053e-001 804data8 0x3fe47a1527e8a2d0 // log(1/frcpa(1+229/256))= +6.39903e-001 805data8 0x3fe489445efffcc8 // log(1/frcpa(1+230/256))= +6.41756e-001 806data8 0x3fe4a018bcb69834 // log(1/frcpa(1+231/256))= +6.44543e-001 807data8 0x3fe4af5a0c9d65d4 // log(1/frcpa(1+232/256))= +6.46405e-001 808data8 0x3fe4bea2a5bdbe84 // log(1/frcpa(1+233/256))= +6.48271e-001 809data8 0x3fe4cdf28f10ac44 // log(1/frcpa(1+234/256))= +6.50140e-001 810data8 0x3fe4dd49cf994058 // log(1/frcpa(1+235/256))= +6.52013e-001 811data8 0x3fe4eca86e64a680 // log(1/frcpa(1+236/256))= +6.53889e-001 812data8 0x3fe503c43cd8eb68 // log(1/frcpa(1+237/256))= +6.56710e-001 813data8 0x3fe513356667fc54 // log(1/frcpa(1+238/256))= +6.58595e-001 814data8 0x3fe522ae0738a3d4 // log(1/frcpa(1+239/256))= +6.60483e-001 815data8 0x3fe5322e26867854 // log(1/frcpa(1+240/256))= +6.62376e-001 816data8 0x3fe541b5cb979808 // log(1/frcpa(1+241/256))= +6.64271e-001 817data8 0x3fe55144fdbcbd60 // log(1/frcpa(1+242/256))= +6.66171e-001 818data8 0x3fe560dbc45153c4 // log(1/frcpa(1+243/256))= +6.68074e-001 819data8 0x3fe5707a26bb8c64 // log(1/frcpa(1+244/256))= +6.69980e-001 820data8 0x3fe587f60ed5b8fc // log(1/frcpa(1+245/256))= +6.72847e-001 821data8 0x3fe597a7977c8f30 // log(1/frcpa(1+246/256))= +6.74763e-001 822data8 0x3fe5a760d634bb88 // log(1/frcpa(1+247/256))= +6.76682e-001 823data8 0x3fe5b721d295f10c // log(1/frcpa(1+248/256))= +6.78605e-001 824data8 0x3fe5c6ea94431ef8 // log(1/frcpa(1+249/256))= +6.80532e-001 825data8 0x3fe5d6bb22ea86f4 // log(1/frcpa(1+250/256))= +6.82462e-001 826data8 0x3fe5e6938645d38c // log(1/frcpa(1+251/256))= +6.84397e-001 827data8 0x3fe5f673c61a2ed0 // log(1/frcpa(1+252/256))= +6.86335e-001 828data8 0x3fe6065bea385924 // log(1/frcpa(1+253/256))= +6.88276e-001 829data8 0x3fe6164bfa7cc068 // log(1/frcpa(1+254/256))= +6.90222e-001 830data8 0x3fe62643fecf9740 // log(1/frcpa(1+255/256))= +6.92171e-001 831LOCAL_OBJECT_END(pow_Tt) 832 833 834// Table 1 is 2^(index_1/128) where 835// index_1 goes from 0 to 15 836LOCAL_OBJECT_START(pow_tbl1) 837data8 0x8000000000000000 , 0x00003FFF 838data8 0x80B1ED4FD999AB6C , 0x00003FFF 839data8 0x8164D1F3BC030773 , 0x00003FFF 840data8 0x8218AF4373FC25EC , 0x00003FFF 841data8 0x82CD8698AC2BA1D7 , 0x00003FFF 842data8 0x8383594EEFB6EE37 , 0x00003FFF 843data8 0x843A28C3ACDE4046 , 0x00003FFF 844data8 0x84F1F656379C1A29 , 0x00003FFF 845data8 0x85AAC367CC487B15 , 0x00003FFF 846data8 0x8664915B923FBA04 , 0x00003FFF 847data8 0x871F61969E8D1010 , 0x00003FFF 848data8 0x87DB357FF698D792 , 0x00003FFF 849data8 0x88980E8092DA8527 , 0x00003FFF 850data8 0x8955EE03618E5FDD , 0x00003FFF 851data8 0x8A14D575496EFD9A , 0x00003FFF 852data8 0x8AD4C6452C728924 , 0x00003FFF 853LOCAL_OBJECT_END(pow_tbl1) 854 855 856// Table 2 is 2^(index_1/8) where 857// index_2 goes from 0 to 7 858LOCAL_OBJECT_START(pow_tbl2) 859data8 0x8000000000000000 , 0x00003FFF 860data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF 861data8 0x9837F0518DB8A96F , 0x00003FFF 862data8 0xA5FED6A9B15138EA , 0x00003FFF 863data8 0xB504F333F9DE6484 , 0x00003FFF 864data8 0xC5672A115506DADD , 0x00003FFF 865data8 0xD744FCCAD69D6AF4 , 0x00003FFF 866data8 0xEAC0C6E7DD24392F , 0x00003FFF 867LOCAL_OBJECT_END(pow_tbl2) 868 869.section .text 870WEAK_LIBM_ENTRY(powf) 871 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;; 884 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;; 897 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;; 910 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;; 923 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 interchangably 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;; 937 938// Continue normal and denormal paths here 939POW_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;; 952 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;; 960 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;; 972 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;; 986 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;; 993 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;; 1006 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;; 1018 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;; 1031 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;; 1044 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;; 1056 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;; 1068 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;; 1080 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;; 1093 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;; 1105 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;; 1120 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;; 1132 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;; 1146 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;; 1161 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;; 1169 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;; 1183 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;; 1197 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;; 1210 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;; 1222 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;; 1229 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;; 1241 1242{ .mmi 1243 ldfe POW_T1 = [pow_AD_T1] 1244 ldfe POW_T2 = [pow_AD_T2] 1245 nop.i 999 1246} 1247;; 1248 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;; 1261 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;; 1268 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;; 1275 1276// p8 TRUE ==> |Y(G + r)| >= 7 1277 1278// single 1279// -2^7 -2^6 2^6 2^7 1280// -----+-----+----+ ... +-----+-----+----- 1281// p8 | p9 | p8 1282// | | p10 | | 1283 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;; 1296 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;; 1310 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;; 1322 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;; 1334 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;; 1346 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;; 1354 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;; 1361 1362// POSSIBLE_OVER_UNDER 1363// p6 = TRUE ==> Y_Gpr negative 1364// Result is already computed. We just need to know if over/underflow occurred. 1365 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;; 1372 1373// POSSIBLE_OVER 1374// We got an answer. 1375// overflow is a possibility, not a certainty 1376 1377 1378// We define an overflow when the answer with 1379// WRE set 1380// user-defined rounding mode 1381 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. 1388 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. 1395 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 1403 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;; 1411 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;; 1418 1419// Return s2 to default 1420{ .mfi 1421 nop.m 999 1422 fsetc.s2 0x7F,0x40 1423 nop.i 999 1424} 1425;; 1426 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;; 1434 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;; 1441 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;; 1448 1449 1450POW_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 1454 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 1472 1473// Form small constant (2^-170) to correct underflow result near region of 1474// smallest denormal in round-nearest. 1475 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;; 1489 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;; 1496 1497// Return s2 to default 1498{ .mfi 1499 nop.m 999 1500 fsetc.s2 0x7F,0x40 1501 nop.i 999 1502} 1503;; 1504 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;; 1512 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;; 1520 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;; 1534 1535POW_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;; 1544 1545{ .mib 1546 getf.sig pow_GR_sig_X = POW_NORM_X 1547 nop.i 999 1548 br.cond.sptk POW_COMMON 1549} 1550;; 1551 1552POW_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;; 1578 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;; 1586 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;; 1599 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;; 1606 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;; 1614 1615POW_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 1618 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;; 1625 1626POW_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 1631 1632// +0 !-odd int +inf 1633// -0 !-odd int +inf 1634 1635// p6 == Y is a floating point number outside the integer. 1636// Hence it is an integer and is even. 1637// return +inf 1638 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// 1648 1649{ .mfi 1650 nop.m 999 1651 nop.f 999 1652 cmp.gt p6,p7 = pow_GR_exp_Y, pow_GR_10033 1653} 1654;; 1655 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;; 1662 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;; 1669 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;; 1676 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;; 1684 1685{ .mfi 1686 nop.m 999 1687(p12) frcpa.s0 f8,p13 = f1,f0 1688 nop.i 999 1689} 1690;; 1691 1692{ .mfb 1693 nop.m 999 1694(p11) frcpa.s0 f8,p13 = f1,f8 1695 br.cond.sptk __libm_error_region 1696} 1697;; 1698 1699 1700POW_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;; 1715 1716 1717POW_X_INF: 1718// Here when X is +-inf 1719 1720// X +inf Y +inf +inf 1721// X -inf Y +inf +inf 1722 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 1726 1727// X +inf Y -inf +0 1728// X -inf Y -inf +0 1729 1730// X +inf Y <0 +0 1731// X -inf Y <0, !odd integer +0 1732// X -inf Y <0, odd integer -0 1733 1734// X + inf Y=+0 +1 1735// X + inf Y=-0 +1 1736// X - inf Y=+0 +1 1737// X - inf Y=-0 +1 1738 1739// p13 == Y negative 1740// p14 == Y positive 1741 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 1748 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 1761 1762// pxx == Y is not an integer 1763// p13 == (Y negative) 1764// return +inf 1765// p14 == (Y positive) 1766// return +0 1767// 1768 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;; 1776 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;; 1788 1789{ .mfi 1790 nop.m 999 1791 fclass.m p15,p0 = f9, 0x07 //@zero 1792 nop.i 999 1793} 1794;; 1795 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;; 1802 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;; 1814 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;; 1822 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;; 1829 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;; 1836 1837 1838POW_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 1842 1843{ .mfb 1844 mov pow_GR_tag = 34 1845 frcpa.s0 f8,p6 = f0,f0 1846 br.cond.sptk __libm_error_region 1847} 1848;; 1849 1850POW_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;; 1858 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;; 1865 1866POW_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;; 1880 1881 1882POW_OVER_UNDER_X_NOT_INF: 1883 1884// p8 is TRUE for overflow 1885// p9 is TRUE for underflow 1886 1887// if y is infinity, we should not over/underflow 1888 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;; 1895 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;; 1907 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;; 1915 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;; 1923 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;; 1930 1931 1932POW_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;; 1941 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;; 1948 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;; 1960 1961 1962POW_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 1965 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;; 1972 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;; 1980 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;; 1988 1989{ .mfi 1990 nop.m 999 1991 fma.s.s0 f8 = POW_tmp, POW_tmp, f0 1992 nop.i 999 1993} 1994;; 1995 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;; 2003 2004WEAK_LIBM_END(powf) 2005libm_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 2012 2013 2014LOCAL_LIBM_ENTRY(__libm_error_region) 2015 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};; 2029 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};; 2036 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};; 2048 2049{ .mmi 2050 add GR_Parameter_RESULT = 48,sp 2051 nop.m 0 2052 nop.i 0 2053};; 2054 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};; 2061 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};; 2067 2068LOCAL_LIBM_END(__libm_error_region) 2069 2070.type __libm_error_support#,@function 2071.global __libm_error_support# 2072