1.file "exp_m1.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// 04/04/00 Unwind support added 43// 08/15/00 Bundle added after call to __libm_error_support to properly 44// set [the previously overwritten] GR_Parameter_RESULT. 45// 07/07/01 Improved speed of all paths 46// 05/20/02 Cleaned up namespace and sf0 syntax 47// 11/20/02 Improved speed, algorithm based on exp 48// 03/31/05 Reformatted delimiters between data tables 49 50// API 51//============================================================== 52// double expm1(double) 53 54// Overview of operation 55//============================================================== 56// 1. Inputs of Nan, Inf, Zero, NatVal handled with special paths 57// 58// 2. |x| < 2^-60 59// Result = x, computed by x + x*x to handle appropriate flags and rounding 60// 61// 3. 2^-60 <= |x| < 2^-2 62// Result determined by 13th order Taylor series polynomial 63// expm1f(x) = x + Q2*x^2 + ... + Q13*x^13 64// 65// 4. x < -48.0 66// Here we know result is essentially -1 + eps, where eps only affects 67// rounded result. Set I. 68// 69// 5. x >= 709.7827 70// Result overflows. Set I, O, and call error support 71// 72// 6. 2^-2 <= x < 709.7827 or -48.0 <= x < -2^-2 73// This is the main path. The algorithm is described below: 74 75// Take the input x. w is "how many log2/128 in x?" 76// w = x * 128/log2 77// n = int(w) 78// x = n log2/128 + r + delta 79 80// n = 128M + index_1 + 2^4 index_2 81// x = M log2 + (log2/128) index_1 + (log2/8) index_2 + r + delta 82 83// exp(x) = 2^M 2^(index_1/128) 2^(index_2/8) exp(r) exp(delta) 84// Construct 2^M 85// Get 2^(index_1/128) from table_1; 86// Get 2^(index_2/8) from table_2; 87// Calculate exp(r) by series by 5th order polynomial 88// r = x - n (log2/128)_high 89// delta = - n (log2/128)_low 90// Calculate exp(delta) as 1 + delta 91 92 93// Special values 94//============================================================== 95// expm1(+0) = +0.0 96// expm1(-0) = -0.0 97 98// expm1(+qnan) = +qnan 99// expm1(-qnan) = -qnan 100// expm1(+snan) = +qnan 101// expm1(-snan) = -qnan 102 103// expm1(-inf) = -1.0 104// expm1(+inf) = +inf 105 106// Overflow and Underflow 107//======================= 108// expm1(x) = largest double normal when 109// x = 709.7827 = 40862e42fefa39ef 110// 111// Underflow is handled as described in case 2 above. 112 113 114// Registers used 115//============================================================== 116// Floating Point registers used: 117// f8, input 118// f9 -> f15, f32 -> f75 119 120// General registers used: 121// r14 -> r40 122 123// Predicate registers used: 124// p6 -> p15 125 126// Assembly macros 127//============================================================== 128 129rRshf = r14 130rAD_TB1 = r15 131rAD_T1 = r15 132rAD_TB2 = r16 133rAD_T2 = r16 134rAD_Ln2_lo = r17 135rAD_P = r17 136 137rN = r18 138rIndex_1 = r19 139rIndex_2_16 = r20 140 141rM = r21 142rBiased_M = r21 143rIndex_1_16 = r22 144rSignexp_x = r23 145rExp_x = r24 146rSig_inv_ln2 = r25 147 148rAD_Q1 = r26 149rAD_Q2 = r27 150rTmp = r27 151rExp_bias = r28 152rExp_mask = r29 153rRshf_2to56 = r30 154 155rGt_ln = r31 156rExp_2tom56 = r31 157 158 159GR_SAVE_B0 = r33 160GR_SAVE_PFS = r34 161GR_SAVE_GP = r35 162GR_SAVE_SP = r36 163 164GR_Parameter_X = r37 165GR_Parameter_Y = r38 166GR_Parameter_RESULT = r39 167GR_Parameter_TAG = r40 168 169 170FR_X = f10 171FR_Y = f1 172FR_RESULT = f8 173 174fRSHF_2TO56 = f6 175fINV_LN2_2TO63 = f7 176fW_2TO56_RSH = f9 177f2TOM56 = f11 178fP5 = f12 179fP54 = f50 180fP5432 = f50 181fP4 = f13 182fP3 = f14 183fP32 = f14 184fP2 = f15 185 186fLn2_by_128_hi = f33 187fLn2_by_128_lo = f34 188 189fRSHF = f35 190fNfloat = f36 191fW = f37 192fR = f38 193fF = f39 194 195fRsq = f40 196fRcube = f41 197 198f2M = f42 199fS1 = f43 200fT1 = f44 201 202fMIN_DBL_OFLOW_ARG = f45 203fMAX_DBL_MINUS_1_ARG = f46 204fMAX_DBL_NORM_ARG = f47 205fP_lo = f51 206fP_hi = f52 207fP = f53 208fS = f54 209 210fNormX = f56 211 212fWre_urm_f8 = f57 213 214fGt_pln = f58 215fTmp = f58 216 217fS2 = f59 218fT2 = f60 219fSm1 = f61 220 221fXsq = f62 222fX6 = f63 223fX4 = f63 224fQ7 = f64 225fQ76 = f64 226fQ7654 = f64 227fQ765432 = f64 228fQ6 = f65 229fQ5 = f66 230fQ54 = f66 231fQ4 = f67 232fQ3 = f68 233fQ32 = f68 234fQ2 = f69 235fQD = f70 236fQDC = f70 237fQDCBA = f70 238fQDCBA98 = f70 239fQDCBA98765432 = f70 240fQC = f71 241fQB = f72 242fQBA = f72 243fQA = f73 244fQ9 = f74 245fQ98 = f74 246fQ8 = f75 247 248// Data tables 249//============================================================== 250 251RODATA 252.align 16 253 254// ************* DO NOT CHANGE ORDER OF THESE TABLES ******************** 255 256// double-extended 1/ln(2) 257// 3fff b8aa 3b29 5c17 f0bb be87fed0691d3e88 258// 3fff b8aa 3b29 5c17 f0bc 259// For speed the significand will be loaded directly with a movl and setf.sig 260// and the exponent will be bias+63 instead of bias+0. Thus subsequent 261// computations need to scale appropriately. 262// The constant 128/ln(2) is needed for the computation of w. This is also 263// obtained by scaling the computations. 264// 265// Two shifting constants are loaded directly with movl and setf.d. 266// 1. fRSHF_2TO56 = 1.1000..00 * 2^(63-7) 267// This constant is added to x*1/ln2 to shift the integer part of 268// x*128/ln2 into the rightmost bits of the significand. 269// The result of this fma is fW_2TO56_RSH. 270// 2. fRSHF = 1.1000..00 * 2^(63) 271// This constant is subtracted from fW_2TO56_RSH * 2^(-56) to give 272// the integer part of w, n, as a floating-point number. 273// The result of this fms is fNfloat. 274 275 276LOCAL_OBJECT_START(exp_Table_1) 277data8 0x40862e42fefa39f0 // smallest dbl overflow arg 278data8 0xc048000000000000 // approx largest arg for minus one result 279data8 0x40862e42fefa39ef // largest dbl arg to give normal dbl result 280data8 0x0 // pad 281data8 0xb17217f7d1cf79ab , 0x00003ff7 // ln2/128 hi 282data8 0xc9e3b39803f2f6af , 0x00003fb7 // ln2/128 lo 283// 284// Table 1 is 2^(index_1/128) where 285// index_1 goes from 0 to 15 286// 287data8 0x8000000000000000 , 0x00003FFF 288data8 0x80B1ED4FD999AB6C , 0x00003FFF 289data8 0x8164D1F3BC030773 , 0x00003FFF 290data8 0x8218AF4373FC25EC , 0x00003FFF 291data8 0x82CD8698AC2BA1D7 , 0x00003FFF 292data8 0x8383594EEFB6EE37 , 0x00003FFF 293data8 0x843A28C3ACDE4046 , 0x00003FFF 294data8 0x84F1F656379C1A29 , 0x00003FFF 295data8 0x85AAC367CC487B15 , 0x00003FFF 296data8 0x8664915B923FBA04 , 0x00003FFF 297data8 0x871F61969E8D1010 , 0x00003FFF 298data8 0x87DB357FF698D792 , 0x00003FFF 299data8 0x88980E8092DA8527 , 0x00003FFF 300data8 0x8955EE03618E5FDD , 0x00003FFF 301data8 0x8A14D575496EFD9A , 0x00003FFF 302data8 0x8AD4C6452C728924 , 0x00003FFF 303LOCAL_OBJECT_END(exp_Table_1) 304 305// Table 2 is 2^(index_1/8) where 306// index_2 goes from 0 to 7 307LOCAL_OBJECT_START(exp_Table_2) 308data8 0x8000000000000000 , 0x00003FFF 309data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF 310data8 0x9837F0518DB8A96F , 0x00003FFF 311data8 0xA5FED6A9B15138EA , 0x00003FFF 312data8 0xB504F333F9DE6484 , 0x00003FFF 313data8 0xC5672A115506DADD , 0x00003FFF 314data8 0xD744FCCAD69D6AF4 , 0x00003FFF 315data8 0xEAC0C6E7DD24392F , 0x00003FFF 316LOCAL_OBJECT_END(exp_Table_2) 317 318 319LOCAL_OBJECT_START(exp_p_table) 320data8 0x3f8111116da21757 //P5 321data8 0x3fa55555d787761c //P4 322data8 0x3fc5555555555414 //P3 323data8 0x3fdffffffffffd6a //P2 324LOCAL_OBJECT_END(exp_p_table) 325 326LOCAL_OBJECT_START(exp_Q1_table) 327data8 0x3de6124613a86d09 // QD = 1/13! 328data8 0x3e21eed8eff8d898 // QC = 1/12! 329data8 0x3ec71de3a556c734 // Q9 = 1/9! 330data8 0x3efa01a01a01a01a // Q8 = 1/8! 331data8 0x8888888888888889,0x3ff8 // Q5 = 1/5! 332data8 0xaaaaaaaaaaaaaaab,0x3ffc // Q3 = 1/3! 333data8 0x0,0x0 // Pad to avoid bank conflicts 334LOCAL_OBJECT_END(exp_Q1_table) 335 336LOCAL_OBJECT_START(exp_Q2_table) 337data8 0x3e5ae64567f544e4 // QB = 1/11! 338data8 0x3e927e4fb7789f5c // QA = 1/10! 339data8 0x3f2a01a01a01a01a // Q7 = 1/7! 340data8 0x3f56c16c16c16c17 // Q6 = 1/6! 341data8 0xaaaaaaaaaaaaaaab,0x3ffa // Q4 = 1/4! 342data8 0x8000000000000000,0x3ffe // Q2 = 1/2! 343LOCAL_OBJECT_END(exp_Q2_table) 344 345 346.section .text 347GLOBAL_IEEE754_ENTRY(expm1) 348 349{ .mlx 350 getf.exp rSignexp_x = f8 // Must recompute if x unorm 351 movl rSig_inv_ln2 = 0xb8aa3b295c17f0bc // signif of 1/ln2 352} 353{ .mlx 354 addl rAD_TB1 = @ltoff(exp_Table_1), gp 355 movl rRshf_2to56 = 0x4768000000000000 // 1.10000 2^(63+56) 356} 357;; 358 359// We do this fnorm right at the beginning to normalize 360// any input unnormals so that SWA is not taken. 361{ .mfi 362 ld8 rAD_TB1 = [rAD_TB1] 363 fclass.m p6,p0 = f8,0x0b // Test for x=unorm 364 mov rExp_mask = 0x1ffff 365} 366{ .mfi 367 mov rExp_bias = 0xffff 368 fnorm.s1 fNormX = f8 369 mov rExp_2tom56 = 0xffff-56 370} 371;; 372 373// Form two constants we need 374// 1/ln2 * 2^63 to compute w = x * 1/ln2 * 128 375// 1.1000..000 * 2^(63+63-7) to right shift int(w) into the significand 376 377{ .mfi 378 setf.sig fINV_LN2_2TO63 = rSig_inv_ln2 // form 1/ln2 * 2^63 379 fclass.m p8,p0 = f8,0x07 // Test for x=0 380 nop.i 0 381} 382{ .mlx 383 setf.d fRSHF_2TO56 = rRshf_2to56 // Form 1.100 * 2^(63+56) 384 movl rRshf = 0x43e8000000000000 // 1.10000 2^63 for rshift 385} 386;; 387 388{ .mfi 389 setf.exp f2TOM56 = rExp_2tom56 // form 2^-56 for scaling Nfloat 390 fclass.m p9,p0 = f8,0x22 // Test for x=-inf 391 add rAD_TB2 = 0x140, rAD_TB1 // Point to Table 2 392} 393{ .mib 394 add rAD_Q1 = 0x1e0, rAD_TB1 // Point to Q table for small path 395 add rAD_Ln2_lo = 0x30, rAD_TB1 // Point to ln2_by_128_lo 396(p6) br.cond.spnt EXPM1_UNORM // Branch if x unorm 397} 398;; 399 400EXPM1_COMMON: 401{ .mfi 402 ldfpd fMIN_DBL_OFLOW_ARG, fMAX_DBL_MINUS_1_ARG = [rAD_TB1],16 403 fclass.m p10,p0 = f8,0x1e1 // Test for x=+inf, NaN, NaT 404 add rAD_Q2 = 0x50, rAD_Q1 // Point to Q table for small path 405} 406{ .mfb 407 nop.m 0 408 nop.f 0 409(p8) br.ret.spnt b0 // Exit for x=0, return x 410} 411;; 412 413{ .mfi 414 ldfd fMAX_DBL_NORM_ARG = [rAD_TB1],16 415 nop.f 0 416 and rExp_x = rExp_mask, rSignexp_x // Biased exponent of x 417} 418{ .mfb 419 setf.d fRSHF = rRshf // Form right shift const 1.100 * 2^63 420(p9) fms.d.s0 f8 = f0,f0,f1 // quick exit for x=-inf 421(p9) br.ret.spnt b0 422} 423;; 424 425{ .mfi 426 ldfpd fQD, fQC = [rAD_Q1], 16 // Load coeff for small path 427 nop.f 0 428 sub rExp_x = rExp_x, rExp_bias // True exponent of x 429} 430{ .mfb 431 ldfpd fQB, fQA = [rAD_Q2], 16 // Load coeff for small path 432(p10) fma.d.s0 f8 = f8, f1, f0 // For x=+inf, NaN, NaT 433(p10) br.ret.spnt b0 // Exit for x=+inf, NaN, NaT 434} 435;; 436 437{ .mfi 438 ldfpd fQ9, fQ8 = [rAD_Q1], 16 // Load coeff for small path 439 fma.s1 fXsq = fNormX, fNormX, f0 // x*x for small path 440 cmp.gt p7, p8 = -2, rExp_x // Test |x| < 2^(-2) 441} 442{ .mfi 443 ldfpd fQ7, fQ6 = [rAD_Q2], 16 // Load coeff for small path 444 nop.f 0 445 nop.i 0 446} 447;; 448 449{ .mfi 450 ldfe fQ5 = [rAD_Q1], 16 // Load coeff for small path 451 nop.f 0 452 nop.i 0 453} 454{ .mib 455 ldfe fQ4 = [rAD_Q2], 16 // Load coeff for small path 456(p7) cmp.gt.unc p6, p7 = -60, rExp_x // Test |x| < 2^(-60) 457(p7) br.cond.spnt EXPM1_SMALL // Branch if 2^-60 <= |x| < 2^-2 458} 459;; 460 461// W = X * Inv_log2_by_128 462// By adding 1.10...0*2^63 we shift and get round_int(W) in significand. 463// We actually add 1.10...0*2^56 to X * Inv_log2 to do the same thing. 464 465{ .mfi 466 ldfe fLn2_by_128_hi = [rAD_TB1],32 467 fma.s1 fW_2TO56_RSH = fNormX, fINV_LN2_2TO63, fRSHF_2TO56 468 nop.i 0 469} 470{ .mfb 471 ldfe fLn2_by_128_lo = [rAD_Ln2_lo] 472(p6) fma.d.s0 f8 = f8, f8, f8 // If x < 2^-60, result=x+x*x 473(p6) br.ret.spnt b0 // Exit if x < 2^-60 474} 475;; 476 477// Divide arguments into the following categories: 478// Certain minus one p11 - -inf < x <= MAX_DBL_MINUS_1_ARG 479// Possible Overflow p14 - MAX_DBL_NORM_ARG < x < MIN_DBL_OFLOW_ARG 480// Certain Overflow p15 - MIN_DBL_OFLOW_ARG <= x < +inf 481// 482// If the input is really a double arg, then there will never be "Possible 483// Overflow" arguments. 484// 485 486// After that last load, rAD_TB1 points to the beginning of table 1 487 488{ .mfi 489 nop.m 0 490 fcmp.ge.s1 p15,p14 = fNormX,fMIN_DBL_OFLOW_ARG 491 nop.i 0 492} 493;; 494 495{ .mfi 496 add rAD_P = 0x80, rAD_TB2 497 fcmp.le.s1 p11,p0 = fNormX,fMAX_DBL_MINUS_1_ARG 498 nop.i 0 499} 500;; 501 502{ .mfb 503 ldfpd fP5, fP4 = [rAD_P] ,16 504(p14) fcmp.gt.unc.s1 p14,p0 = fNormX,fMAX_DBL_NORM_ARG 505(p15) br.cond.spnt EXPM1_CERTAIN_OVERFLOW 506} 507;; 508 509// Nfloat = round_int(W) 510// The signficand of fW_2TO56_RSH contains the rounded integer part of W, 511// as a twos complement number in the lower bits (that is, it may be negative). 512// That twos complement number (called N) is put into rN. 513 514// Since fW_2TO56_RSH is scaled by 2^56, it must be multiplied by 2^-56 515// before the shift constant 1.10000 * 2^63 is subtracted to yield fNfloat. 516// Thus, fNfloat contains the floating point version of N 517 518{ .mfb 519 ldfpd fP3, fP2 = [rAD_P] 520 fms.s1 fNfloat = fW_2TO56_RSH, f2TOM56, fRSHF 521(p11) br.cond.spnt EXPM1_CERTAIN_MINUS_ONE 522} 523;; 524 525{ .mfi 526 getf.sig rN = fW_2TO56_RSH 527 nop.f 0 528 nop.i 0 529} 530;; 531 532// rIndex_1 has index_1 533// rIndex_2_16 has index_2 * 16 534// rBiased_M has M 535// rIndex_1_16 has index_1 * 16 536 537// r = x - Nfloat * ln2_by_128_hi 538// f = 1 - Nfloat * ln2_by_128_lo 539{ .mfi 540 and rIndex_1 = 0x0f, rN 541 fnma.s1 fR = fNfloat, fLn2_by_128_hi, fNormX 542 shr rM = rN, 0x7 543} 544{ .mfi 545 and rIndex_2_16 = 0x70, rN 546 fnma.s1 fF = fNfloat, fLn2_by_128_lo, f1 547 nop.i 0 548} 549;; 550 551// rAD_T1 has address of T1 552// rAD_T2 has address if T2 553 554{ .mmi 555 add rBiased_M = rExp_bias, rM 556 add rAD_T2 = rAD_TB2, rIndex_2_16 557 shladd rAD_T1 = rIndex_1, 4, rAD_TB1 558} 559;; 560 561// Create Scale = 2^M 562// Load T1 and T2 563{ .mmi 564 setf.exp f2M = rBiased_M 565 ldfe fT2 = [rAD_T2] 566 nop.i 0 567} 568;; 569 570{ .mfi 571 ldfe fT1 = [rAD_T1] 572 fmpy.s0 fTmp = fLn2_by_128_lo, fLn2_by_128_lo // Force inexact 573 nop.i 0 574} 575;; 576 577{ .mfi 578 nop.m 0 579 fma.s1 fP54 = fR, fP5, fP4 580 nop.i 0 581} 582{ .mfi 583 nop.m 0 584 fma.s1 fP32 = fR, fP3, fP2 585 nop.i 0 586} 587;; 588 589{ .mfi 590 nop.m 0 591 fma.s1 fRsq = fR, fR, f0 592 nop.i 0 593} 594;; 595 596{ .mfi 597 nop.m 0 598 fma.s1 fP5432 = fRsq, fP54, fP32 599 nop.i 0 600} 601;; 602 603{ .mfi 604 nop.m 0 605 fma.s1 fS2 = fF,fT2,f0 606 nop.i 0 607} 608{ .mfi 609 nop.m 0 610 fma.s1 fS1 = f2M,fT1,f0 611 nop.i 0 612} 613;; 614 615{ .mfi 616 nop.m 0 617 fma.s1 fP = fRsq, fP5432, fR 618 nop.i 0 619} 620;; 621 622{ .mfi 623 nop.m 0 624 fms.s1 fSm1 = fS1,fS2,f1 // S - 1.0 625 nop.i 0 626} 627{ .mfb 628 nop.m 0 629 fma.s1 fS = fS1,fS2,f0 630(p14) br.cond.spnt EXPM1_POSSIBLE_OVERFLOW 631} 632;; 633 634{ .mfb 635 nop.m 0 636 fma.d.s0 f8 = fS, fP, fSm1 637 br.ret.sptk b0 // Normal path exit 638} 639;; 640 641// Here if 2^-60 <= |x| <2^-2 642// Compute 13th order polynomial 643EXPM1_SMALL: 644{ .mmf 645 ldfe fQ3 = [rAD_Q1], 16 646 ldfe fQ2 = [rAD_Q2], 16 647 fma.s1 fX4 = fXsq, fXsq, f0 648} 649;; 650 651{ .mfi 652 nop.m 0 653 fma.s1 fQDC = fQD, fNormX, fQC 654 nop.i 0 655} 656{ .mfi 657 nop.m 0 658 fma.s1 fQBA = fQB, fNormX, fQA 659 nop.i 0 660} 661;; 662 663{ .mfi 664 nop.m 0 665 fma.s1 fQ98 = fQ9, fNormX, fQ8 666 nop.i 0 667} 668{ .mfi 669 nop.m 0 670 fma.s1 fQ76= fQ7, fNormX, fQ6 671 nop.i 0 672} 673;; 674 675{ .mfi 676 nop.m 0 677 fma.s1 fQ54 = fQ5, fNormX, fQ4 678 nop.i 0 679} 680;; 681 682{ .mfi 683 nop.m 0 684 fma.s1 fX6 = fX4, fXsq, f0 685 nop.i 0 686} 687{ .mfi 688 nop.m 0 689 fma.s1 fQ32= fQ3, fNormX, fQ2 690 nop.i 0 691} 692;; 693 694{ .mfi 695 nop.m 0 696 fma.s1 fQDCBA = fQDC, fXsq, fQBA 697 nop.i 0 698} 699{ .mfi 700 nop.m 0 701 fma.s1 fQ7654 = fQ76, fXsq, fQ54 702 nop.i 0 703} 704;; 705 706{ .mfi 707 nop.m 0 708 fma.s1 fQDCBA98 = fQDCBA, fXsq, fQ98 709 nop.i 0 710} 711{ .mfi 712 nop.m 0 713 fma.s1 fQ765432 = fQ7654, fXsq, fQ32 714 nop.i 0 715} 716;; 717 718{ .mfi 719 nop.m 0 720 fma.s1 fQDCBA98765432 = fQDCBA98, fX6, fQ765432 721 nop.i 0 722} 723;; 724 725{ .mfb 726 nop.m 0 727 fma.d.s0 f8 = fQDCBA98765432, fXsq, fNormX 728 br.ret.sptk b0 // Exit small branch 729} 730;; 731 732 733EXPM1_POSSIBLE_OVERFLOW: 734 735// Here if fMAX_DBL_NORM_ARG < x < fMIN_DBL_OFLOW_ARG 736// This cannot happen if input is a double, only if input higher precision. 737// Overflow is a possibility, not a certainty. 738 739// Recompute result using status field 2 with user's rounding mode, 740// and wre set. If result is larger than largest double, then we have 741// overflow 742 743{ .mfi 744 mov rGt_ln = 0x103ff // Exponent for largest dbl + 1 ulp 745 fsetc.s2 0x7F,0x42 // Get user's round mode, set wre 746 nop.i 0 747} 748;; 749 750{ .mfi 751 setf.exp fGt_pln = rGt_ln // Create largest double + 1 ulp 752 fma.d.s2 fWre_urm_f8 = fS, fP, fSm1 // Result with wre set 753 nop.i 0 754} 755;; 756 757{ .mfi 758 nop.m 0 759 fsetc.s2 0x7F,0x40 // Turn off wre in sf2 760 nop.i 0 761} 762;; 763 764{ .mfi 765 nop.m 0 766 fcmp.ge.s1 p6, p0 = fWre_urm_f8, fGt_pln // Test for overflow 767 nop.i 0 768} 769;; 770 771{ .mfb 772 nop.m 0 773 nop.f 0 774(p6) br.cond.spnt EXPM1_CERTAIN_OVERFLOW // Branch if overflow 775} 776;; 777 778{ .mfb 779 nop.m 0 780 fma.d.s0 f8 = fS, fP, fSm1 781 br.ret.sptk b0 // Exit if really no overflow 782} 783;; 784 785EXPM1_CERTAIN_OVERFLOW: 786{ .mmi 787 sub rTmp = rExp_mask, r0, 1 788;; 789 setf.exp fTmp = rTmp 790 nop.i 0 791} 792;; 793 794{ .mfi 795 alloc r32=ar.pfs,1,4,4,0 796 fmerge.s FR_X = f8,f8 797 nop.i 0 798} 799{ .mfb 800 mov GR_Parameter_TAG = 41 801 fma.d.s0 FR_RESULT = fTmp, fTmp, f0 // Set I,O and +INF result 802 br.cond.sptk __libm_error_region 803} 804;; 805 806// Here if x unorm 807EXPM1_UNORM: 808{ .mfb 809 getf.exp rSignexp_x = fNormX // Must recompute if x unorm 810 fcmp.eq.s0 p6, p0 = f8, f0 // Set D flag 811 br.cond.sptk EXPM1_COMMON 812} 813;; 814 815// here if result will be -1 and inexact, x <= -48.0 816EXPM1_CERTAIN_MINUS_ONE: 817{ .mmi 818 mov rTmp = 1 819;; 820 setf.exp fTmp = rTmp 821 nop.i 0 822} 823;; 824 825{ .mfb 826 nop.m 0 827 fms.d.s0 FR_RESULT = fTmp, fTmp, f1 // Set I, rounded -1+eps result 828 br.ret.sptk b0 829} 830;; 831 832GLOBAL_IEEE754_END(expm1) 833libm_alias_double_other (__expm1, expm1) 834 835 836LOCAL_LIBM_ENTRY(__libm_error_region) 837.prologue 838{ .mfi 839 add GR_Parameter_Y=-32,sp // Parameter 2 value 840 nop.f 0 841.save ar.pfs,GR_SAVE_PFS 842 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs 843} 844{ .mfi 845.fframe 64 846 add sp=-64,sp // Create new stack 847 nop.f 0 848 mov GR_SAVE_GP=gp // Save gp 849};; 850{ .mmi 851 stfd [GR_Parameter_Y] = FR_Y,16 // STORE Parameter 2 on stack 852 add GR_Parameter_X = 16,sp // Parameter 1 address 853.save b0, GR_SAVE_B0 854 mov GR_SAVE_B0=b0 // Save b0 855};; 856.body 857{ .mib 858 stfd [GR_Parameter_X] = FR_X // STORE Parameter 1 on stack 859 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address 860 nop.b 0 861} 862{ .mib 863 stfd [GR_Parameter_Y] = FR_RESULT // STORE Parameter 3 on stack 864 add GR_Parameter_Y = -16,GR_Parameter_Y 865 br.call.sptk b0=__libm_error_support# // Call error handling function 866};; 867{ .mmi 868 add GR_Parameter_RESULT = 48,sp 869 nop.m 0 870 nop.i 0 871};; 872{ .mmi 873 ldfd f8 = [GR_Parameter_RESULT] // Get return result off stack 874.restore sp 875 add sp = 64,sp // Restore stack pointer 876 mov b0 = GR_SAVE_B0 // Restore return address 877};; 878{ .mib 879 mov gp = GR_SAVE_GP // Restore gp 880 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs 881 br.ret.sptk b0 // Return 882};; 883 884LOCAL_LIBM_END(__libm_error_region) 885.type __libm_error_support#,@function 886.global __libm_error_support# 887