1.file "erfcf.s" 2 3 4// Copyright (c) 2002 - 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// 01/17/02 Initial version 42// 05/20/02 Cleaned up namespace and sf0 syntax 43// 02/06/03 Reordered header: .section, .global, .proc, .align 44// 03/31/05 Reformatted delimiters between data tables 45// 46// API 47//============================================================== 48// float erfcf(float) 49// 50// Overview of operation 51//============================================================== 52// 1. 0 <= x <= 10.06 53// 54// erfcf(x) = P15(x) * exp( -x^2 ) 55// 56// Comment: 57// 58// Let x(0)=0, x(i) = 2^(i), i=1,...3, x(4)= 10.06 59// 60// Let x(i)<= x < x(i+1). 61// We can find i as exponent of argument x (let i = 0 for 0<= x < 2 ) 62// 63// Let P15(x) - polynomial approximation of degree 15 for function 64// erfcf(x) * exp( x^2) and x(i) <= x <= x(i+1), i = 0,1,2,3 65// Polynomial coefficients we have in the table erfc_p_table. 66// 67// So we can find result for erfcf(x) as above. 68// Algorithm description for exp function see below. 69// 70// 2. -4.4 <= x < 0 71// 72// erfcf(x) = 2.0 - erfcf(-x) 73// 74// 3. x > 10.06 75// 76// erfcf(x) ~=~ 0.0 77// 78// 4. x < -4.4 79// 80// erfcf(x) ~=~ 2.0 81 82// Special values 83//============================================================== 84// erfcf(+0) = 1.0 85// erfcf(-0) = 1.0 86 87// erfcf(+qnan) = +qnan 88// erfcf(-qnan) = -qnan 89// erfcf(+snan) = +qnan 90// erfcf(-snan) = -qnan 91 92// erfcf(-inf) = 2.0 93// erfcf(+inf) = +0 94 95//============================================================== 96// Take double exp(double) from libm_64. 97// 98// Overview of operation 99//============================================================== 100// Take the input x. w is "how many log2/128 in x?" 101// w = x * 128/log2 102// n = int(w) 103// x = n log2/128 + r + delta 104 105// n = 128M + index_1 + 2^4 index_2 106// x = M log2 + (log2/128) index_1 + (log2/8) index_2 + r + delta 107 108// exp(x) = 2^M 2^(index_1/128) 2^(index_2/8) exp(r) exp(delta) 109// Construct 2^M 110// Get 2^(index_1/128) from table_1; 111// Get 2^(index_2/8) from table_2; 112// Calculate exp(r) by series 113// r = x - n (log2/128)_high 114// delta = - n (log2/128)_low 115// Calculate exp(delta) as 1 + delta 116// 117// Comment for erfcf: 118// 119// Let exp(r) = 1 + x + 0.5*x^2 + (1/6)*x^3 120// Let delta = 0. 121//============================================================== 122// 123// Registers used 124//============================================================== 125// Floating Point registers used: 126// f8, input 127// f6,f7,f9 -> f11, f32 -> f92 128 129// General registers used: 130// r14 -> r22,r32 -> r50 131 132// Predicate registers used: 133// p6 -> p15 134 135// Assembly macros 136//============================================================== 137EXP_AD_TB1 = r14 138exp_GR_sig_inv_ln2 = r15 139exp_TB1_size = r16 140exp_GR_rshf_2to56 = r17 141exp_GR_exp_2tom56 = r18 142 143exp_GR_rshf = r33 144EXP_AD_TB2 = r34 145EXP_AD_P = r35 146exp_GR_N = r36 147exp_GR_index_1 = r37 148exp_GR_index_2_16 = r38 149exp_GR_biased_M = r39 150EXP_AD_T1 = r40 151EXP_AD_T2 = r41 152exp_TB2_size = r42 153 154// GR for erfcf(x) 155//============================================================== 156GR_IndxPlusBias = r19 157GR_ExpMask = r20 158GR_BIAS = r21 159GR_ShftPi_bias = r22 160 161GR_P_POINT_1 = r43 162GR_P_POINT_2 = r44 163GR_P_POINT_3 = r45 164GR_P_POINT_4 = r46 165 166GR_ShftPi = r47 167GR_EpsNorm = r48 168 169GR_05 = r49 170GR_1_by_6 = r50 171 172// GR for __libm_support call 173//============================================================== 174 175GR_SAVE_B0 = r43 176GR_SAVE_PFS = r44 177GR_SAVE_GP = r45 178GR_SAVE_SP = r46 179 180GR_Parameter_X = r47 181GR_Parameter_Y = r48 182GR_Parameter_RESULT = r49 183GR_Parameter_TAG = r50 184 185 186// FR for exp(-x^2) 187//============================================================== 188FR_X = f10 189FR_Y = f1 190FR_RESULT = f8 191 192EXP_2TOM56 = f6 193EXP_INV_LN2_2TO63 = f7 194EXP_W_2TO56_RSH = f9 195exp_ln2_by_128_hi = f11 196 197EXP_RSHF_2TO56 = f32 198exp_ln2_by_128_lo = f33 199EXP_RSHF = f34 200EXP_Nfloat = f35 201exp_r = f36 202exp_rsq = f37 203EXP_2M = f38 204exp_S1 = f39 205exp_T1 = f40 206exp_P = f41 207exp_S = f42 208EXP_NORM_f8 = f43 209exp_S2 = f44 210exp_T2 = f45 211 212// FR for erfcf(x) 213//============================================================== 214FR_AbsArg = f46 215FR_Tmp = f47 216FR_Tmp1 = f48 217FR_Tmpf = f49 218FR_NormX = f50 219 220FR_A15 = f51 221FR_A14 = f52 222 223FR_A13 = f53 224FR_A12 = f54 225 226FR_A11 = f55 227FR_A10 = f56 228 229FR_A9 = f57 230FR_A8 = f58 231 232FR_A7 = f59 233FR_A6 = f60 234 235FR_A5 = f61 236FR_A4 = f62 237 238FR_A3 = f63 239FR_A2 = f64 240 241FR_A1 = f65 242FR_A0 = f66 243 244FR_P15_0_1 = f67 245FR_P15_1_1 = f68 246FR_P15_1_2 = f69 247FR_P15_2_1 = f70 248FR_P15_2_2 = f71 249FR_P15_3_1 = f72 250FR_P15_3_2 = f73 251FR_P15_4_1 = f74 252FR_P15_4_2 = f75 253FR_P15_7_1 = f76 254FR_P15_7_2 = f77 255FR_P15_8_1 = f78 256FR_P15_9_1 = f79 257FR_P15_9_2 = f80 258FR_P15_13_1 = f81 259FR_P15_14_1 = f82 260FR_P15_14_2 = f83 261 262FR_2 = f84 263FR_05 = f85 264FR_1_by_6 = f86 265FR_Pol = f87 266FR_Exp = f88 267 268FR_POS_ARG_ASYMP = f89 269FR_NEG_ARG_ASYMP = f90 270 271FR_UnfBound = f91 272FR_EpsNorm = f92 273 274// Data tables 275//============================================================== 276RODATA 277.align 16 278 279// ************* DO NOT CHANGE ORDER OF THESE TABLES ******************** 280 281// double-extended 1/ln(2) 282// 3fff b8aa 3b29 5c17 f0bb be87fed0691d3e88 283// 3fff b8aa 3b29 5c17 f0bc 284// For speed the significand will be loaded directly with a movl and setf.sig 285// and the exponent will be bias+63 instead of bias+0. Thus subsequent 286// computations need to scale appropriately. 287// The constant 128/ln(2) is needed for the computation of w. This is also 288// obtained by scaling the computations. 289// 290// Two shifting constants are loaded directly with movl and setf.d. 291// 1. EXP_RSHF_2TO56 = 1.1000..00 * 2^(63-7) 292// This constant is added to x*1/ln2 to shift the integer part of 293// x*128/ln2 into the rightmost bits of the significand. 294// The result of this fma is EXP_W_2TO56_RSH. 295// 2. EXP_RSHF = 1.1000..00 * 2^(63) 296// This constant is subtracted from EXP_W_2TO56_RSH * 2^(-56) to give 297// the integer part of w, n, as a floating-point number. 298// The result of this fms is EXP_Nfloat. 299 300 301LOCAL_OBJECT_START(exp_table_1) 302 303data4 0x4120f5c3, 0x408ccccd //POS_ARG_ASYMP = 10.06, NEG_ARG_ASYMP = 4.4 304data4 0x41131Cdf, 0x00800000 //UnfBound ~=~ 9.1, EpsNorm ~=~ 1.1754944e-38 305// 306data8 0xb17217f7d1cf79ab , 0x00003ff7 // ln2/128 hi 307data8 0xc9e3b39803f2f6af , 0x00003fb7 // ln2/128 lo 308// 309// Table 1 is 2^(index_1/128) where 310// index_1 goes from 0 to 15 311// 312data8 0x8000000000000000 , 0x00003FFF 313data8 0x80B1ED4FD999AB6C , 0x00003FFF 314data8 0x8164D1F3BC030773 , 0x00003FFF 315data8 0x8218AF4373FC25EC , 0x00003FFF 316data8 0x82CD8698AC2BA1D7 , 0x00003FFF 317data8 0x8383594EEFB6EE37 , 0x00003FFF 318data8 0x843A28C3ACDE4046 , 0x00003FFF 319data8 0x84F1F656379C1A29 , 0x00003FFF 320data8 0x85AAC367CC487B15 , 0x00003FFF 321data8 0x8664915B923FBA04 , 0x00003FFF 322data8 0x871F61969E8D1010 , 0x00003FFF 323data8 0x87DB357FF698D792 , 0x00003FFF 324data8 0x88980E8092DA8527 , 0x00003FFF 325data8 0x8955EE03618E5FDD , 0x00003FFF 326data8 0x8A14D575496EFD9A , 0x00003FFF 327data8 0x8AD4C6452C728924 , 0x00003FFF 328LOCAL_OBJECT_END(exp_table_1) 329 330// Table 2 is 2^(index_1/8) where 331// index_2 goes from 0 to 7 332 333LOCAL_OBJECT_START(exp_table_2) 334 335data8 0x8000000000000000 , 0x00003FFF 336data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF 337data8 0x9837F0518DB8A96F , 0x00003FFF 338data8 0xA5FED6A9B15138EA , 0x00003FFF 339data8 0xB504F333F9DE6484 , 0x00003FFF 340data8 0xC5672A115506DADD , 0x00003FFF 341data8 0xD744FCCAD69D6AF4 , 0x00003FFF 342data8 0xEAC0C6E7DD24392F , 0x00003FFF 343LOCAL_OBJECT_END(exp_table_2) 344 345LOCAL_OBJECT_START(erfc_p_table) 346 347// Pol_0 348data8 0xBEA3260C63CB0446 //A15 = -5.70673541831883454676e-07 349data8 0x3EE63D6178077654 //A14 = +1.06047480138940182343e-05 350data8 0xBF18646BC5FC70A7 //A13 = -9.30491237309283694347e-05 351data8 0x3F40F92F909117FE //A12 = +5.17986512144075019133e-04 352data8 0xBF611344289DE1E6 //A11 = -2.08438217390159994419e-03 353data8 0x3F7AF9FE6AD16DC0 //A10 = +6.58606893292862351928e-03 354data8 0xBF91D219E196CBA7 //A9 = -1.74030345858217321001e-02 355data8 0x3FA4AFDDA355854C //A8 = +4.04042493708041968315e-02 356data8 0xBFB5D465BB7025AE //A7 = -8.52721769916999425445e-02 357data8 0x3FC54C15A95B717D //A6 = +1.66384418195672549029e-01 358data8 0xBFD340A75B4B1AB5 //A5 = -3.00821150926292166899e-01 359data8 0x3FDFFFC0BFCD247F //A4 = +4.99984919839853542841e-01 360data8 0xBFE81270C361852B //A3 = -7.52251035312075583309e-01 361data8 0x3FEFFFFFC67295FC //A2 = +9.99999892800303301771e-01 362data8 0xBFF20DD74F8CD2BF //A1 = -1.12837916445020868099e+00 363data8 0x3FEFFFFFFFFE7C1D //A0 = +9.99999999988975570714e-01 364// Pol_1 365data8 0xBDE8EC4BDD953B56 //A15 = -1.81338928934942767144e-10 366data8 0x3E43607F269E2A1C //A14 = +9.02309090272196442358e-09 367data8 0xBE8C4D9E69C10E02 //A13 = -2.10875261143659275328e-07 368data8 0x3EC9CF2F84566725 //A12 = +3.07671055805877356583e-06 369data8 0xBF007980B1B46A4D //A11 = -3.14228438702169818945e-05 370data8 0x3F2F4C3AD6DEF24A //A10 = +2.38783056770846320260e-04 371data8 0xBF56F5129F8D30FA //A9 = -1.40120333363130546426e-03 372data8 0x3F7AA6C7ABFC38EE //A8 = +6.50671002200751820429e-03 373data8 0xBF98E7522CB84BEF //A7 = -2.43199195666185511109e-02 374data8 0x3FB2F68EB1C3D073 //A6 = +7.40746673580490638637e-02 375data8 0xBFC7C16055AC6385 //A5 = -1.85588876564704611769e-01 376data8 0x3FD8A707AEF5A440 //A4 = +3.85194702967570635211e-01 377data8 0xBFE547BFE39AE2EA //A3 = -6.65008492032112467310e-01 378data8 0x3FEE7C91BDF13578 //A2 = +9.52706213932898128515e-01 379data8 0xBFF1CB5B61F8C589 //A1 = -1.11214769621105541214e+00 380data8 0x3FEFEA56BC81FD37 //A0 = +9.97355812243688815239e-01 381// Pol_2 382data8 0xBD302724A12F46E0 //A15 = -5.73866382814058809406e-14 383data8 0x3D98889B75D3102E //A14 = +5.57829983681360947356e-12 384data8 0xBDF16EA15074A1E9 //A13 = -2.53671153922423457844e-10 385data8 0x3E3EC6E688CFEE5F //A12 = +7.16581828336436419561e-09 386data8 0xBE82E5ED44C52609 //A11 = -1.40802202239825487803e-07 387data8 0x3EC120BE5CE42353 //A10 = +2.04180535157522081699e-06 388data8 0xBEF7B8B0311A1911 //A9 = -2.26225266204633600888e-05 389data8 0x3F29A281F43FC238 //A8 = +1.95577968156184077632e-04 390data8 0xBF55E19858B3B7A4 //A7 = -1.33552434527526534043e-03 391data8 0x3F7DAC8C3D12E5FD //A6 = +7.24463253680473816303e-03 392data8 0xBF9FF9C04613FB47 //A5 = -3.12261622211693854028e-02 393data8 0x3FBB3D5DBF9D9366 //A4 = +1.06405123978743883370e-01 394data8 0xBFD224DE9F62C258 //A3 = -2.83500342989133623476e-01 395data8 0x3FE28A95CB8C6D3E //A2 = +5.79417131000276437708e-01 396data8 0xBFEC21205D358672 //A1 = -8.79043752717008257224e-01 397data8 0x3FEDAE44D5EDFE5B //A0 = +9.27523057776805771830e-01 398// Pol_3 399data8 0xBCA3BCA734AC82F1 //A15 = -1.36952437983096410260e-16 400data8 0x3D16740DC3990612 //A14 = +1.99425676175410093285e-14 401data8 0xBD77F4353812C46A //A13 = -1.36162367755616790260e-12 402data8 0x3DCFD0BE13C73DB4 //A12 = +5.78718761040355136007e-11 403data8 0xBE1D728DF71189B4 //A11 = -1.71406885583934105120e-09 404data8 0x3E64252C8CB710B5 //A10 = +3.75233795940731111303e-08 405data8 0xBEA514B93180F33D //A9 = -6.28261292774310809962e-07 406data8 0x3EE1381118CC7151 //A8 = +8.21066421390821904504e-06 407data8 0xBF1634404FB0FA72 //A7 = -8.47019436358372148764e-05 408data8 0x3F46B2CBBCF0EB32 //A6 = +6.92700845213200923490e-04 409data8 0xBF725C2B445E6D81 //A5 = -4.48243046949004063741e-03 410data8 0x3F974E7CFA4D89D9 //A4 = +2.27603462002522228717e-02 411data8 0xBFB6D7BAC2E342D1 //A3 = -8.92292714882032736443e-02 412data8 0x3FD0D156AD9CE2A6 //A2 = +2.62777013343603696631e-01 413data8 0xBFE1C228572AADB0 //A1 = -5.54950876471982857725e-01 414data8 0x3FE8A739F48B9A3B //A0 = +7.70413377406675619766e-01 415LOCAL_OBJECT_END(erfc_p_table) 416 417 418.section .text 419GLOBAL_LIBM_ENTRY(erfcf) 420 421// Form index i for table erfc_p_table as exponent of x 422// We use i + bias in real calculations 423{ .mlx 424 getf.exp GR_IndxPlusBias = f8 // (sign + exp + bias) of x 425 movl exp_GR_sig_inv_ln2 = 0xb8aa3b295c17f0bc //signif.of 1/ln2 426} 427{ .mlx 428 addl EXP_AD_TB1 = @ltoff(exp_table_1), gp 429 movl exp_GR_rshf_2to56 = 0x4768000000000000 // 1.100 2^(63+56) 430} 431;; 432 433// Form argument EXP_NORM_f8 for exp(-x^2) 434{ .mfi 435 ld8 EXP_AD_TB1 = [EXP_AD_TB1] 436 fcmp.ge.s1 p6,p7 = f8, f0 // p6: x >= 0 ,p7: x<0 437 mov GR_BIAS = 0x0FFFF 438} 439{ .mfi 440 mov exp_GR_exp_2tom56 = 0xffff-56 441 fnma.s1 EXP_NORM_f8 = f8, f8, f0 // -x^2 442 mov GR_ExpMask = 0x1ffff 443} 444;; 445 446// Form two constants we need 447// 1/ln2 * 2^63 to compute w = x * 1/ln2 * 128 448// 1.1000..000 * 2^(63+63-7) to right shift int(w) into the significand 449 450// p9: x = 0,+inf,-inf,nan,unnorm. 451// p10: x!= 0,+inf,-inf,nan,unnorm. 452{ .mfi 453 setf.sig EXP_INV_LN2_2TO63 = exp_GR_sig_inv_ln2 // Form 1/ln2*2^63 454 fclass.m p9,p10 = f8,0xef 455 shl GR_ShftPi_bias = GR_BIAS, 7 456} 457{ .mfi 458 setf.d EXP_RSHF_2TO56 = exp_GR_rshf_2to56 //Const 1.10*2^(63+56) 459 nop.f 0 460 and GR_IndxPlusBias = GR_IndxPlusBias, GR_ExpMask // i + bias 461} 462;; 463 464{ .mfi 465 alloc r32 = ar.pfs, 0, 15, 4, 0 466(p6) fma.s1 FR_AbsArg = f1, f0, f8 // |x| if x >= 0 467 cmp.lt p15,p0 = GR_IndxPlusBias, GR_BIAS//p15: i < 0 (for |x|<1) 468} 469{ .mlx 470 setf.exp EXP_2TOM56 = exp_GR_exp_2tom56 //2^-56 for scaling Nfloat 471 movl exp_GR_rshf = 0x43e8000000000000 //1.10 2^63,right shift. 472} 473;; 474 475{ .mfi 476 ldfps FR_POS_ARG_ASYMP, FR_NEG_ARG_ASYMP = [EXP_AD_TB1],8 477 nop.f 0 478(p15) mov GR_IndxPlusBias = GR_BIAS //Let i = 0 if i < 0 479} 480{ .mlx 481 mov GR_P_POINT_3 = 0x1A0 482 movl GR_05 = 0x3fe0000000000000 483} 484;; 485 486// Form shift GR_ShftPi from the beginning of erfc_p_table 487// to the polynomial with number i 488{ .mfi 489 ldfps FR_UnfBound, FR_EpsNorm = [EXP_AD_TB1],8 490 nop.f 0 491 shl GR_ShftPi = GR_IndxPlusBias, 7 492} 493{ .mfi 494 setf.d EXP_RSHF = exp_GR_rshf // Form right shift 1.100 * 2^63 495(p7) fms.s1 FR_AbsArg = f1, f0, f8 // |x| if x < 0 496 mov exp_TB1_size = 0x100 497} 498;; 499 500// Form pointer GR_P_POINT_3 to the beginning of erfc_p_table 501{ .mfi 502 setf.d FR_05 = GR_05 503 nop.f 0 504 sub GR_ShftPi = GR_ShftPi,GR_ShftPi_bias 505} 506{ .mfb 507 add GR_P_POINT_3 = GR_P_POINT_3, EXP_AD_TB1 508 nop.f 0 509(p9) br.cond.spnt SPECIAL // For x = 0,+inf,-inf,nan,unnorm 510} 511;; 512 513{ .mfi 514 add GR_P_POINT_1 = GR_P_POINT_3, GR_ShftPi 515 nop.f 0 516 add GR_P_POINT_2 = GR_P_POINT_3, GR_ShftPi 517} 518{ .mfi 519 ldfe exp_ln2_by_128_hi = [EXP_AD_TB1],16 520 fma.s1 FR_NormX = f8,f1,f0 521 add GR_P_POINT_3 = GR_P_POINT_3, GR_ShftPi 522} 523;; 524 525// Load coefficients for polynomial P15(x) 526{ .mfi 527 ldfpd FR_A15, FR_A14 = [GR_P_POINT_1], 16 528 nop.f 0 529 add GR_P_POINT_3 = 0x30, GR_P_POINT_3 530} 531{ .mfi 532 ldfe exp_ln2_by_128_lo = [EXP_AD_TB1], 16 533 nop.f 0 534 add GR_P_POINT_2 = 0x20, GR_P_POINT_2 535} 536;; 537 538// Now EXP_AD_TB1 points to the beginning of table 1 539{ .mlx 540 ldfpd FR_A13, FR_A12 = [GR_P_POINT_1] 541 movl GR_1_by_6 = 0x3FC5555555555555 542} 543{ .mfi 544 add GR_P_POINT_4 = 0x30, GR_P_POINT_2 545 nop.f 0 546 nop.i 0 547} 548;; 549 550{ .mfi 551 ldfpd FR_A11, FR_A10 = [GR_P_POINT_2] 552 fma.s1 FR_2 = f1, f1, f1 553 mov exp_TB2_size = 0x80 554} 555{ .mfi 556 ldfpd FR_A9, FR_A8 = [GR_P_POINT_3],16 557 nop.f 0 558 add GR_P_POINT_1 = 0x60 ,GR_P_POINT_1 559} 560;; 561 562// W = X * Inv_log2_by_128 563// By adding 1.10...0*2^63 we shift and get round_int(W) in significand. 564// We actually add 1.10...0*2^56 to X * Inv_log2 to do the same thing. 565{ .mfi 566 ldfpd FR_A7, FR_A6 = [GR_P_POINT_3] 567 fma.s1 EXP_W_2TO56_RSH = EXP_NORM_f8,EXP_INV_LN2_2TO63,EXP_RSHF_2TO56 568 add EXP_AD_TB2 = exp_TB1_size, EXP_AD_TB1 569 570} 571{ .mfi 572 ldfpd FR_A5, FR_A4 = [GR_P_POINT_4], 16 573 nop.f 0 574 nop.i 0 575} 576;; 577 578{ .mfi 579 ldfpd FR_A3, FR_A2 = [GR_P_POINT_4] 580 fmerge.s FR_X = f8,f8 581 nop.i 0 582} 583{ .mfi 584 ldfpd FR_A1, FR_A0 = [GR_P_POINT_1] 585 nop.f 0 586 nop.i 0 587} 588;; 589 590//p14: x < - NEG_ARG_ASYMP = -4.4 -> erfcf(x) ~=~ 2.0 591{ .mfi 592 setf.d FR_1_by_6 = GR_1_by_6 593(p7) fcmp.gt.unc.s1 p14,p0 = FR_AbsArg, FR_NEG_ARG_ASYMP //p7: x < 0 594 nop.i 0 595} 596;; 597 598//p15: x > POS_ARG_ASYMP = 10.06 -> erfcf(x) ~=~ 0.0 599{ .mfi 600 nop.m 0 601(p6) fcmp.gt.unc.s1 p15,p0 = FR_AbsArg, FR_POS_ARG_ASYMP //p6: x > 0 602 nop.i 0 603} 604;; 605 606{ .mfi 607 nop.m 0 608 fcmp.le.s1 p8,p0 = FR_NormX, FR_UnfBound // p8: x <= UnfBound 609 nop.i 0 610} 611{ .mfb 612 nop.m 0 613(p14) fnma.s.s0 FR_RESULT = FR_EpsNorm, FR_EpsNorm, FR_2//y = 2 if x <-4.4 614(p14) br.ret.spnt b0 615} 616;; 617 618// Nfloat = round_int(W) 619// The signficand of EXP_W_2TO56_RSH contains the rounded integer part of W, 620// as a twos complement number in the lower bits (that is, it may be negative). 621// That twos complement number (called N) is put into exp_GR_N. 622 623// Since EXP_W_2TO56_RSH is scaled by 2^56, it must be multiplied by 2^-56 624// before the shift constant 1.10000 * 2^63 is subtracted to yield EXP_Nfloat. 625// Thus, EXP_Nfloat contains the floating point version of N 626 627{ .mfi 628 nop.m 0 629 fms.s1 EXP_Nfloat = EXP_W_2TO56_RSH, EXP_2TOM56, EXP_RSHF 630 nop.i 0 631} 632{ .mfb 633(p15) mov GR_Parameter_TAG = 209 634(p15) fma.s.s0 FR_RESULT = FR_EpsNorm,FR_EpsNorm,f0 //Result.for x>10.06 635(p15) br.cond.spnt __libm_error_region 636} 637;; 638 639// Now we can calculate polynomial P15(x) 640{ .mfi 641 nop.m 0 642 fma.s1 FR_P15_1_1 = FR_AbsArg, FR_AbsArg, f0 // x ^2 643 nop.i 0 644} 645{ .mfi 646 nop.m 0 647 fma.s1 FR_P15_0_1 = FR_A15, FR_AbsArg, FR_A14 648 nop.i 0 649} 650;; 651 652{ .mfi 653 nop.m 0 654 fma.s1 FR_P15_1_2 = FR_A13, FR_AbsArg, FR_A12 655 nop.i 0 656} 657;; 658 659{ .mfi 660 getf.sig exp_GR_N = EXP_W_2TO56_RSH 661 fma.s1 FR_P15_2_1 = FR_A9, FR_AbsArg, FR_A8 662 nop.i 0 663} 664{ .mfi 665 nop.m 0 666 fma.s1 FR_P15_2_2 = FR_A11, FR_AbsArg, FR_A10 667 nop.i 0 668} 669;; 670 671{ .mfi 672 nop.m 0 673 fma.s1 FR_P15_3_1 = FR_A5, FR_AbsArg, FR_A4 674 nop.i 0 675} 676{ .mfi 677 nop.m 0 678 fma.s1 FR_P15_3_2 = FR_A7, FR_AbsArg, FR_A6 679 nop.i 0 680} 681;; 682 683// exp_GR_index_1 has index_1 684// exp_GR_index_2_16 has index_2 * 16 685// exp_GR_biased_M has M 686// exp_GR_index_1_16 has index_1 * 16 687 688// r2 has true M 689{ .mfi 690 and exp_GR_index_1 = 0x0f, exp_GR_N 691 fma.s1 FR_P15_4_1 = FR_A1, FR_AbsArg, FR_A0 692 shr r2 = exp_GR_N, 0x7 693 694} 695{ .mfi 696 and exp_GR_index_2_16 = 0x70, exp_GR_N 697 fma.s1 FR_P15_4_2 = FR_A3, FR_AbsArg, FR_A2 698 nop.i 0 699} 700;; 701 702// EXP_AD_T1 has address of T1 703// EXP_AD_T2 has address if T2 704 705{ .mfi 706 add EXP_AD_T2 = EXP_AD_TB2, exp_GR_index_2_16 707 nop.f 0 708 shladd EXP_AD_T1 = exp_GR_index_1, 4, EXP_AD_TB1 709} 710{ .mfi 711 addl exp_GR_biased_M = 0xffff, r2 712 fnma.s1 exp_r = EXP_Nfloat, exp_ln2_by_128_hi, EXP_NORM_f8 713 nop.i 0 714} 715;; 716 717// Create Scale = 2^M 718// r = x - Nfloat * ln2_by_128_hi 719 720{ .mfi 721 setf.exp EXP_2M = exp_GR_biased_M 722 fma.s1 FR_P15_7_1 = FR_P15_0_1, FR_P15_1_1, FR_P15_1_2 723 nop.i 0 724} 725{ .mfi 726 ldfe exp_T2 = [EXP_AD_T2] 727 nop.f 0 728 nop.i 0 729} 730;; 731 732// Load T1 and T2 733 734{ .mfi 735 ldfe exp_T1 = [EXP_AD_T1] 736 fma.s1 FR_P15_7_2 = FR_P15_1_1, FR_P15_1_1, f0 // x^4 737 nop.i 0 738} 739{ .mfi 740 nop.m 0 741 fma.s1 FR_P15_8_1 = FR_P15_1_1, FR_P15_2_2, FR_P15_2_1 742 nop.i 0 743} 744;; 745 746{ .mfi 747 nop.m 0 748 fma.s1 FR_P15_9_1 = FR_P15_1_1, FR_P15_4_2, FR_P15_4_1 749 nop.i 0 750} 751{ .mfi 752 nop.m 0 753 fma.s1 FR_P15_9_2 = FR_P15_1_1, FR_P15_3_2, FR_P15_3_1 754 nop.i 0 755} 756;; 757 758{ .mfi 759 nop.m 0 760 fma.s1 exp_P = FR_1_by_6, exp_r, FR_05 761 nop.i 0 762} 763{ .mfi 764 nop.m 0 765 fma.s1 exp_rsq = exp_r, exp_r, f0 766 nop.i 0 767} 768;; 769 770{ .mfi 771 nop.m 0 772 fma.s1 FR_P15_13_1 = FR_P15_7_2, FR_P15_7_1, FR_P15_8_1 773 nop.i 0 774} 775;; 776 777{ .mfi 778 nop.m 0 779 fma.s1 FR_P15_14_1 = FR_P15_7_2, FR_P15_9_2, FR_P15_9_1 780 nop.i 0 781} 782{ .mfi 783 nop.m 0 784 fma.s1 FR_P15_14_2 = FR_P15_7_2, FR_P15_7_2, f0 // x^8 785 nop.i 0 786} 787;; 788 789{ .mfi 790 nop.m 0 791 fma.s1 exp_P = exp_P, exp_rsq, exp_r 792 nop.i 0 793} 794{ .mfi 795 nop.m 0 796 fma.s1 exp_S1 = EXP_2M, exp_T2, f0 797 nop.i 0 798} 799;; 800 801{ .mfi 802 nop.m 0 803 fma.s1 FR_Pol = FR_P15_14_2, FR_P15_13_1, FR_P15_14_1 // P15(x) 804 nop.i 0 805} 806;; 807 808{ .mfi 809 nop.m 0 810 fma.s1 exp_S = exp_S1, exp_T1, f0 811 nop.i 0 812} 813;; 814 815{ .mfi 816 nop.m 0 817 fma.s1 FR_Exp = exp_S, exp_P, exp_S // exp(-x^2) 818 nop.i 0 819} 820;; 821 822{ .mfi 823 nop.m 0 824 fma.s.s0 FR_Tmpf = f8, f1, f0 // Flag d 825 nop.i 0 826} 827;; 828 829//p6: result for 0 < x < = POS_ARG_ASYMP 830//p7: result for - NEG_ARG_ASYMP <= x < 0 831//p8: exit for - NEG_ARG_ASYMP <= x <= UnfBound, x!=0 832.pred.rel "mutex",p6,p7 833{ .mfi 834 nop.m 0 835(p6) fma.s.s0 f8 = FR_Exp, FR_Pol, f0 836 nop.i 0 837} 838{ .mfb 839 mov GR_Parameter_TAG = 209 840(p7) fnma.s.s0 f8 = FR_Exp, FR_Pol, FR_2 841(p8) br.ret.sptk b0 842} 843;; 844 845//p10: branch for UnfBound < x < = POS_ARG_ASYMP 846{ .mfb 847 nop.m 0 848 nop.f 0 849(p10) br.cond.spnt __libm_error_region 850} 851;; 852 853//Only via (p9) br.cond.spnt SPECIAL for x = 0,+inf,-inf,nan,unnorm 854SPECIAL: 855 856{ .mfi 857 nop.m 0 858 fclass.m.unc p10,p0 = f8,0x07 // p10: x = 0 859 nop.i 0 860} 861;; 862 863{ .mfi 864 nop.m 0 865 fclass.m.unc p11,p0 = f8,0x21 // p11: x = +inf 866 nop.i 0 867} 868;; 869 870{ .mfi 871 nop.m 0 872 fclass.m.unc p12,p0 = f8,0x22 // p12 x = -inf 873 nop.i 0 874} 875{ .mfb 876 nop.m 0 877(p10) fma.s.s0 f8 = f1, f1, f0 878(p10) br.ret.sptk b0 // Quick exit for x = 0 879} 880;; 881 882{ .mfi 883 nop.m 0 884 fclass.m.unc p13,p0 = f8,0xc3 // p13: x = nan 885 nop.i 0 886} 887{ .mfb 888 nop.m 0 889(p11) fma.s.s0 f8 = f0, f1, f0 890(p11) br.ret.spnt b0 // Quick exit for x = +inf 891} 892;; 893{ .mfi 894 nop.m 0 895 fclass.m.unc p14,p0 = f8,0x0b // P14: x = unnormalized 896 nop.i 0 897} 898{ .mfb 899 nop.m 0 900(p12) fma.s.s0 f8 = f1, f1, f1 901(p12) br.ret.spnt b0 // Quick exit for x = -inf 902} 903;; 904 905{ .mfb 906 nop.m 0 907(p13) fma.s.s0 f8 = f8, f1, f0 908(p13) br.ret.sptk b0 // Quick exit for x = nan 909} 910;; 911 912{ .mfb 913 nop.m 0 914(p14) fnma.s.s0 f8 = f8, f1, f1 915(p14) br.ret.sptk b0 // Quick exit for x = unnormalized 916} 917;; 918 919GLOBAL_LIBM_END(erfcf) 920libm_alias_float_other (erfc, erfc) 921 922 923// Call via (p10) br.cond.spnt __libm_error_region 924// for UnfBound < x < = POS_ARG_ASYMP 925// and 926// 927// call via (p15) br.cond.spnt __libm_error_region 928// for x > POS_ARG_ASYMP 929 930LOCAL_LIBM_ENTRY(__libm_error_region) 931.prologue 932{ .mfi 933 add GR_Parameter_Y=-32,sp // Parameter 2 value 934 nop.f 0 935.save ar.pfs,GR_SAVE_PFS 936 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs 937} 938{ .mfi 939.fframe 64 940 add sp=-64,sp // Create new stack 941 nop.f 0 942 mov GR_SAVE_GP=gp // Save gp 943};; 944{ .mmi 945 stfs [GR_Parameter_Y] = FR_Y,16 // STORE Parameter 2 on stack 946 add GR_Parameter_X = 16,sp // Parameter 1 address 947.save b0, GR_SAVE_B0 948 mov GR_SAVE_B0=b0 // Save b0 949};; 950.body 951{ .mib 952 stfs [GR_Parameter_X] = FR_X // STORE Parameter 1 on stack 953 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address 954 nop.b 0 955} 956{ .mib 957 stfs [GR_Parameter_Y] = FR_RESULT // STORE Parameter 3 on stack 958 add GR_Parameter_Y = -16,GR_Parameter_Y 959 br.call.sptk b0=__libm_error_support# // Call error handling function 960};; 961{ .mmi 962 nop.m 0 963 nop.m 0 964 add GR_Parameter_RESULT = 48,sp 965};; 966{ .mmi 967 ldfs f8 = [GR_Parameter_RESULT] // Get return result off stack 968.restore sp 969 add sp = 64,sp // Restore stack pointer 970 mov b0 = GR_SAVE_B0 // Restore return address 971};; 972{ .mib 973 mov gp = GR_SAVE_GP // Restore gp 974 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs 975 br.ret.sptk b0 // Return 976};; 977 978LOCAL_LIBM_END(__libm_error_region) 979.type __libm_error_support#,@function 980.global __libm_error_support# 981