sysdeps/ia64/fpu/e_fmodl.S

   1 .file "fmodl.s"
   2
   3
   4 // Copyright (c) 2000 - 2004, Intel Corporation
   5 // All rights reserved.
   6 //
   7 // Contributed 2000 by the Intel Numerics Group, Intel Corporation
   8 //
   9 // Redistribution and use in source and binary forms, with or without
  10 // modification, are permitted provided that the following conditions are
  11 // met:
  12 //
  13 // * Redistributions of source code must retain the above copyright
  14 // notice, this list of conditions and the following disclaimer.
  15 //
  16 // * Redistributions in binary form must reproduce the above copyright
  17 // notice, this list of conditions and the following disclaimer in the
  18 // documentation and/or other materials provided with the distribution.
  19 //
  20 // * The name of Intel Corporation may not be used to endorse or promote
  21 // products derived from this software without specific prior written
  22 // permission.
  23
  24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
  28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  29 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  30 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  31 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
  33 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  34 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  35 //
  36 // Intel Corporation is the author of this code, and requests that all
  37 // problem reports or change requests be submitted to it directly at
  38 // http://www.intel.com/software/products/opensource/libraries/num.htm.
  39 //
  40 // History
  41 //====================================================================
  42 // 02/02/00 Initial version
  43 // 03/02/00 New Algorithm
  44 // 04/04/00 Unwind support added
  45 // 08/15/00 Bundle added after call to __libm_error_support to properly
  46 // set [ the previously overwritten ] GR_Parameter_RESULT.
  47 // 11/28/00 Set FR_Y to f9
  48 // 03/11/02 Fixed flags for fmodl(qnan, zero)
  49 // 05/20/02 Cleaned up namespace and sf0 syntax
  50 // 02/10/03 Reordered header:.section,.global,.proc,.align
  51 // 04/28/03 Fix: fmod(sNaN, 0) no longer sets errno
  52 // 11/23/04 Reformatted routine and improved speed
  53 //
  54 // API
  55 //====================================================================
  56 // long double fmodl(long double, long double);
  57 //
  58 // Overview of operation
  59 //====================================================================
  60 // fmod(a, b)= a-i*b,
  61 // where i is an integer such that, if b!= 0,
  62 // |i|<|a/b| and |a/b-i|<1
  63 //
  64 // Algorithm
  65 //====================================================================
  66 // a). if |a|<|b|, return a
  67 // b). get quotient and reciprocal overestimates accurate to
  68 // 33 bits (q2, y2)
  69 // c). if the exponent difference (exponent(a)-exponent(b))
  70 // is less than 32, truncate quotient to integer and
  71 // finish in one iteration
  72 // d). if exponent(a)-exponent(b)>= 32 (q2>= 2^32)
  73 // round quotient estimate to single precision (k= RN(q2)),
  74 // calculate partial remainder (a'= a-k*b),
  75 // get quotient estimate (a'*y2), and repeat from c).
  76 //
  77 // Registers used
  78 //====================================================================
  79
  80 GR_SMALLBIASEXP     = r2
  81 GR_2P32             = r3
  82 GR_SMALLBIASEXP     = r20
  83 GR_ROUNDCONST       = r21
  84 GR_SIG_B            = r22
  85 GR_ARPFS            = r23
  86 GR_TMP1             = r24
  87 GR_TMP2             = r25
  88 GR_TMP3             = r26
  89
  90 GR_SAVE_B0          = r33
  91 GR_SAVE_PFS         = r34
  92 GR_SAVE_GP          = r35
  93 GR_SAVE_SP          = r36
  94
  95 GR_Parameter_X      = r37
  96 GR_Parameter_Y      = r38
  97 GR_Parameter_RESULT = r39
  98 GR_Parameter_TAG    = r40
  99
 100 FR_X                = f10
 101 FR_Y                = f9
 102 FR_RESULT           = f8
 103
 104 FR_ABS_A            = f6
 105 FR_ABS_B            = f7
 106 FR_Y_INV            = f10
 107 FR_SMALLBIAS        = f11
 108 FR_E0               = f12
 109 FR_Q                = f13
 110 FR_E1               = f14
 111 FR_2P32             = f15
 112 FR_TMPX             = f32
 113 FR_TMPY             = f33
 114 FR_ROUNDCONST       = f34
 115 FR_QINT             = f35
 116 FR_QRND24           = f36
 117 FR_NORM_B           = f37
 118 FR_TMP              = f38
 119 FR_TMP2             = f39
 120 FR_DFLAG            = f40
 121 FR_Y_INV0           = f41
 122 FR_Y_INV1           = f42
 123 FR_Q0               = f43
 124 FR_Q1               = f44
 125 FR_QINT_Z           = f45
 126 FR_QREM             = f46
 127 FR_B_SGN_A          = f47
 128
 129 .section .text
 130 GLOBAL_IEEE754_ENTRY(fmodl)
 131
 132 // inputs in f8, f9
 133 // result in f8
 134
 135 { .mfi
 136        getf.sig GR_SIG_B = f9
 137        // FR_ABS_A = |a|
 138        fmerge.s FR_ABS_A = f0, f8
 139        mov GR_SMALLBIASEXP = 0x0ffdd
 140 }
 141 { .mfi
 142        nop.m 0
 143        // FR_ABS_B = |b|
 144        fmerge.s FR_ABS_B = f0, f9
 145        nop.i 0
 146 }
 147 ;;
 148
 149 { .mfi
 150        setf.exp FR_SMALLBIAS = GR_SMALLBIASEXP
 151        // (1) y0
 152        frcpa.s1 FR_Y_INV0, p6 = FR_ABS_A, FR_ABS_B
 153        nop.i 0
 154 }
 155 ;;
 156
 157 { .mlx
 158        nop.m 0
 159        movl GR_ROUNDCONST = 0x33a00000
 160 }
 161 ;;
 162
 163 // eliminate special cases
 164 { .mmi
 165        nop.m 0
 166        nop.m 0
 167        // y pseudo-zero ?
 168        cmp.eq p7, p10 = GR_SIG_B, r0
 169 }
 170 ;;
 171
 172 // set p7 if b +/-NAN, +/-inf, +/-0
 173 { .mfi
 174        nop.m 0
 175  (p10) fclass.m p7, p10 = f9, 0xe7
 176        nop.i 0
 177 }
 178 ;;
 179
 180 { .mfi
 181        mov GR_2P32 = 0x1001f
 182        // (2) q0 = a*y0
 183  (p6)  fma.s1 FR_Q0 = FR_ABS_A, FR_Y_INV0, f0
 184        nop.i 0
 185 }
 186 { .mfi
 187        nop.m 0
 188        // (3) e0 = 1 - b * y0
 189  (p6)  fnma.s1 FR_E0 = FR_ABS_B, FR_Y_INV0, f1
 190        nop.i 0
 191 }
 192 ;;
 193
 194 // set p9 if a +/-NAN, +/-inf
 195 { .mfi
 196        nop.m 0
 197        fclass.m.unc p9, p11 = f8, 0xe3
 198        nop.i 0
 199 }
 200        // |a| < |b|? Return a, p8=1
 201 { .mfi
 202        nop.m 0
 203  (p10) fcmp.lt.unc.s1 p8, p0 = FR_ABS_A, FR_ABS_B
 204        nop.i 0
 205 }
 206 ;;
 207
 208 // set p7 if b +/-NAN, +/-inf, +/-0
 209 { .mfi
 210        nop.m 0
 211        // pseudo-NaN ?
 212  (p10) fclass.nm p7, p0 = f9, 0xff
 213        nop.i 0
 214 }
 215 ;;
 216
 217 // set p9 if a is +/-NaN, +/-Inf
 218 { .mfi
 219        nop.m 0
 220  (p11) fclass.nm p9, p0 = f8, 0xff
 221        nop.i 0
 222 }
 223 { .mfi
 224        nop.m 0
 225        // b denormal ? set D flag (if |a|<|b|)
 226  (p8)  fnma.s0 FR_DFLAG = f9, f1, f9
 227        nop.i 0
 228 }
 229 ;;
 230
 231 { .mfi
 232        // FR_2P32 = 2^32
 233        setf.exp FR_2P32 = GR_2P32
 234        // (4) q1 = q0+e0*q0
 235  (p6)  fma.s1 FR_Q1 = FR_E0, FR_Q0, FR_Q0
 236        nop.i 0
 237 }
 238 { .mfi
 239        nop.m 0
 240        // (5) e1 = e0 * e0 + 2^-34
 241  (p6)  fma.s1 FR_E1 = FR_E0, FR_E0, FR_SMALLBIAS
 242        nop.i 0
 243 }
 244 ;;
 245
 246 { .mfi
 247        nop.m 0
 248        // normalize a (if |a|<|b|)
 249  (p8)  fma.s0 f8 = f8, f1, f0
 250        nop.i 0
 251 }
 252 { .bbb
 253  (p9) br.cond.spnt FMOD_A_NAN_INF
 254  (p7) br.cond.spnt FMOD_B_NAN_INF_ZERO
 255        // if |a|<|b|, return
 256  (p8) br.ret.spnt b0
 257 }
 258 ;;
 259
 260
 261 { .mfi
 262        nop.m 0
 263        // (6) y1 = y0 + e0 * y0
 264  (p6)  fma.s1 FR_Y_INV1 = FR_E0, FR_Y_INV0, FR_Y_INV0
 265        nop.i 0
 266 }
 267 ;;
 268
 269 { .mfi
 270        nop.m 0
 271        // a denormal ? set D flag
 272        // b denormal ? set D flag
 273        fcmp.eq.s0 p12,p0 = FR_ABS_A, FR_ABS_B
 274        nop.i 0
 275 }
 276 { .mfi
 277        // set FR_ROUNDCONST = 1.25*2^{-24}
 278        setf.s FR_ROUNDCONST = GR_ROUNDCONST
 279        // (7) q2 = q1+e1*q1
 280  (p6)  fma.s1 FR_Q = FR_Q1, FR_E1, FR_Q1
 281        nop.i 0
 282 }
 283 ;;
 284
 285 { .mfi
 286        nop.m 0
 287        fmerge.s FR_B_SGN_A = f8, f9
 288        nop.i 0
 289 }
 290 { .mfi
 291        nop.m 0
 292        // (8) y2 = y1 + e1 * y1
 293  (p6)  fma.s1 FR_Y_INV = FR_E1, FR_Y_INV1, FR_Y_INV1
 294        // set p6 = 0, p10 = 0
 295        cmp.ne.and p6, p10 = r0, r0
 296 }
 297 ;;
 298
 299 //   will compute integer quotient bits (24 bits per iteration)
 300 .align 32
 301 loop64:
 302 { .mfi
 303        nop.m 0
 304        // compare q2, 2^32
 305        fcmp.lt.unc.s1 p8, p7 = FR_Q, FR_2P32
 306        nop.i 0
 307 }
 308 { .mfi
 309        nop.m 0
 310        // will truncate quotient to integer, if exponent<32 (in advance)
 311        fcvt.fx.trunc.s1 FR_QINT = FR_Q
 312        nop.i 0
 313 }
 314 ;;
 315
 316 { .mfi
 317        nop.m 0
 318        // if exponent>32 round quotient to single precision (perform in advance)
 319        fma.s.s1 FR_QRND24 = FR_Q, f1, f0
 320        nop.i 0
 321 }
 322 ;;
 323
 324 { .mfi
 325        nop.m 0
 326        // set FR_ROUNDCONST = sgn(a)
 327  (p8)  fmerge.s FR_ROUNDCONST = f8, f1
 328        nop.i 0
 329 }
 330 { .mfi
 331        nop.m 0
 332        // normalize truncated quotient
 333  (p8)  fcvt.xf FR_QRND24 = FR_QINT
 334        nop.i 0
 335 }
 336 ;;
 337
 338 { .mfi
 339        nop.m 0
 340        // calculate remainder (assuming FR_QRND24 = RZ(Q))
 341  (p7)  fnma.s1 FR_E1 = FR_QRND24, FR_ABS_B, FR_ABS_A
 342        nop.i 0
 343 }
 344 { .mfi
 345        nop.m 0
 346        // also if exponent>32, round quotient to single precision
 347        // and subtract 1 ulp: q = q-q*(1.25*2^{-24})
 348  (p7)  fnma.s.s1 FR_QINT_Z = FR_QRND24, FR_ROUNDCONST, FR_QRND24
 349        nop.i 0
 350 }
 351 ;;
 352
 353 { .mfi
 354        nop.m 0
 355        // (p8) calculate remainder (82-bit format)
 356  (p8)  fnma.s1 FR_QREM = FR_QRND24, FR_ABS_B, FR_ABS_A
 357        nop.i 0
 358 }
 359 { .mfi
 360        nop.m 0
 361        // (p7) calculate remainder (assuming FR_QINT_Z = RZ(Q))
 362  (p7)  fnma.s1 FR_ABS_A = FR_QINT_Z, FR_ABS_B, FR_ABS_A
 363        nop.i 0
 364 }
 365 ;;
 366
 367 { .mfi
 368        nop.m 0
 369        // Final iteration (p8): is FR_ABS_A the correct remainder
 370        // (quotient was not overestimated) ?
 371  (p8)  fcmp.lt.unc.s1 p6, p10 = FR_QREM, f0
 372        nop.i 0
 373 }
 374 ;;
 375
 376 { .mfi
 377        nop.m 0
 378        // get new quotient estimation: a'*y2
 379  (p7)  fma.s1 FR_Q = FR_E1, FR_Y_INV, f0
 380        nop.i 0
 381 }
 382 { .mfb
 383        nop.m 0
 384        // was FR_Q = RZ(Q) ? (then new remainder FR_E1> = 0)
 385  (p7)  fcmp.lt.unc.s1 p7, p9 = FR_E1, f0
 386        nop.b 0
 387 }
 388 ;;
 389
 390 .pred.rel "mutex", p6, p10
 391 { .mfb
 392        nop.m 0
 393        // add b to estimated remainder (to cover the case when the quotient was
 394        // overestimated)
 395        // also set correct sign by using
 396        // FR_B_SGN_A = |b|*sgn(a), FR_ROUNDCONST = sgn(a)
 397  (p6)  fma.s0 f8 = FR_QREM, FR_ROUNDCONST, FR_B_SGN_A
 398        nop.b 0
 399 }
 400 { .mfb
 401        nop.m 0
 402        // set correct sign of result before returning: FR_ROUNDCONST = sgn(a)
 403  (p10) fma.s0 f8 = FR_QREM, FR_ROUNDCONST, f0
 404  (p8)  br.ret.sptk b0
 405 }
 406 ;;
 407
 408 { .mfi
 409        nop.m 0
 410        // if f13! = RZ(Q), get alternative quotient estimation: a''*y2
 411  (p7)  fma.s1 FR_Q = FR_ABS_A, FR_Y_INV, f0
 412        nop.i 0
 413 }
 414 { .mfb
 415        nop.m 0
 416        // if FR_E1 was RZ(Q), set remainder to FR_E1
 417  (p9)  fma.s1 FR_ABS_A = FR_E1, f1, f0
 418        br.cond.sptk loop64
 419 }
 420 ;;
 421
 422 FMOD_A_NAN_INF:
 423
 424 // b zero ?
 425 { .mfi
 426        nop.m 0
 427        fclass.m p10, p0 = f8, 0xc3 // Test a = nan
 428        nop.i 0
 429 }
 430 { .mfi
 431        nop.m 0
 432        fma.s1 FR_NORM_B = f9, f1, f0
 433        nop.i 0
 434 }
 435 ;;
 436
 437 { .mfi
 438        nop.m 0
 439        fma.s0 f8 = f8, f1, f0
 440        nop.i 0
 441 }
 442 { .mfi
 443        nop.m 0
 444  (p10) fclass.m p10, p0 = f9, 0x07 // Test x = nan, and y = zero
 445        nop.i 0
 446 }
 447 ;;
 448
 449 { .mfb
 450        nop.m 0
 451        fcmp.eq.unc.s1 p11, p0 = FR_NORM_B, f0
 452  (p10) br.ret.spnt b0 // Exit with result = a if a = nan and b = zero
 453 }
 454 ;;
 455
 456 { .mib
 457        nop.m 0
 458        nop.i 0
 459        // if Y zero
 460  (p11) br.cond.spnt FMOD_B_ZERO
 461 }
 462 ;;
 463
 464 // a= infinity? Return QNAN indefinite
 465 { .mfi
 466        // set p7 t0 0
 467        cmp.ne p7, p0 = r0, r0
 468        fclass.m.unc p8, p9 = f8, 0x23
 469        nop.i 0
 470 }
 471 ;;
 472
 473 // b NaN ?
 474 { .mfi
 475        nop.m 0
 476  (p8)  fclass.m p9, p8 = f9, 0xc3
 477        nop.i 0
 478 }
 479 ;;
 480
 481 // b not pseudo-zero ? (GR_SIG_B holds significand)
 482 { .mii
 483        nop.m 0
 484  (p8)  cmp.ne p7, p0 = GR_SIG_B, r0
 485        nop.i 0
 486 }
 487 ;;
 488
 489 { .mfi
 490        nop.m 0
 491  (p8)  frcpa.s0 f8, p0 = f8, f8
 492        nop.i 0
 493 }
 494 { .mfi
 495        nop.m 0
 496        // also set Denormal flag if necessary
 497  (p7)  fnma.s0 f9 = f9, f1, f9
 498        nop.i 0
 499 }
 500 ;;
 501
 502 { .mfb
 503        nop.m 0
 504  (p8)  fma.s0 f8 = f8, f1, f0
 505        nop.b 0
 506 }
 507 ;;
 508
 509 { .mfb
 510        nop.m 0
 511  (p9)  frcpa.s0 f8, p7 = f8, f9
 512        br.ret.sptk b0
 513 }
 514 ;;
 515
 516 FMOD_B_NAN_INF_ZERO:
 517 // b INF
 518 { .mfi
 519        nop.m 0
 520        fclass.m.unc p7, p0 = f9, 0x23
 521        nop.i 0
 522 }
 523 ;;
 524
 525 { .mfb
 526        nop.m 0
 527  (p7)  fma.s0 f8 = f8, f1, f0
 528  (p7)  br.ret.spnt b0
 529 }
 530 ;;
 531
 532 // b NAN?
 533 { .mfi
 534        nop.m 0
 535        fclass.m.unc p9, p10 = f9, 0xc3
 536        nop.i 0
 537 }
 538 ;;
 539
 540 { .mfi
 541        nop.m 0
 542  (p10) fclass.nm p9, p0 = f9, 0xff
 543        nop.i 0
 544 }
 545 ;;
 546
 547 { .mfb
 548        nop.m 0
 549  (p9)  fma.s0 f8 = f9, f1, f0
 550  (p9)  br.ret.spnt b0
 551 }
 552 ;;
 553
 554 FMOD_B_ZERO:
 555 // Y zero? Must be zero at this point
 556 // because it is the only choice left.
 557 // Return QNAN indefinite
 558
 559 { .mfi
 560        nop.m 0
 561        // set Invalid
 562        frcpa.s0 FR_TMP, p0 = f0, f0
 563        nop.i 0
 564 }
 565 ;;
 566
 567 // a NAN?
 568 { .mfi
 569        nop.m 0
 570        fclass.m.unc p9, p10 = f8, 0xc3
 571        nop.i 0
 572 }
 573 ;;
 574
 575 { .mfi
 576        alloc GR_ARPFS = ar.pfs, 1, 4, 4, 0
 577  (p10) fclass.nm p9, p10 = f8, 0xff
 578        nop.i 0
 579 }
 580 ;;
 581
 582 { .mfi
 583        nop.m 0
 584  (p9)  frcpa.s0 FR_TMP2, p7 = f8, f0
 585        nop.i 0
 586 }
 587 ;;
 588
 589 { .mfi
 590        nop.m 0
 591  (p10) frcpa.s0 FR_TMP2, p7 = f9, f9
 592        mov GR_Parameter_TAG = 120
 593 }
 594 ;;
 595
 596 { .mfi
 597        nop.m 0
 598        fmerge.s FR_X = f8, f8
 599        nop.i 0
 600 }
 601 { .mfb
 602        nop.m 0
 603        fma.s0 f8 = FR_TMP2, f1, f0
 604        br.sptk __libm_error_region
 605 }
 606 ;;
 607
 608 GLOBAL_IEEE754_END(fmodl)
 609
 610 LOCAL_LIBM_ENTRY(__libm_error_region)
 611 .prologue
 612 { .mfi
 613        add GR_Parameter_Y = -32, sp // Parameter 2 value
 614        nop.f 0
 615 .save ar.pfs, GR_SAVE_PFS
 616        mov GR_SAVE_PFS = ar.pfs     // Save ar.pfs
 617 }
 618 { .mfi
 619 .fframe 64
 620        add sp = -64, sp             // Create new stack
 621        nop.f 0
 622        mov GR_SAVE_GP = gp          // Save gp
 623 }
 624 ;;
 625
 626 { .mmi
 627        stfe [ GR_Parameter_Y ] = FR_Y, 16 // Save Parameter 2 on stack
 628        add GR_Parameter_X = 16, sp  // Parameter 1 address
 629 .save b0, GR_SAVE_B0
 630        mov GR_SAVE_B0 = b0          // Save b0
 631 }
 632 ;;
 633
 634 .body
 635 { .mib
 636        stfe [ GR_Parameter_X ] = FR_X // Store Parameter 1 on stack
 637        add GR_Parameter_RESULT = 0, GR_Parameter_Y
 638        nop.b 0                      // Parameter 3 address
 639 }
 640 { .mib
 641        stfe [ GR_Parameter_Y ] = FR_RESULT // Store Parameter 3 on stack
 642        add GR_Parameter_Y = -16, GR_Parameter_Y
 643        br.call.sptk b0 = __libm_error_support# // Call error handling function
 644 }
 645 ;;
 646
 647 { .mmi
 648        nop.m 0
 649        nop.m 0
 650        add GR_Parameter_RESULT = 48, sp
 651 }
 652 ;;
 653
 654 { .mmi
 655        ldfe f8 = [ GR_Parameter_RESULT ] // Get return result off stack
 656 .restore sp
 657        add sp = 64, sp                   // Restore stack pointer
 658        mov b0 = GR_SAVE_B0               // Restore return address
 659 }
 660 ;;
 661
 662 { .mib
 663        mov gp = GR_SAVE_GP               // Restore gp
 664        mov ar.pfs = GR_SAVE_PFS          // Restore ar.pfs
 665        br.ret.sptk b0                    // Return
 666 }
 667 ;;
 668
 669 LOCAL_LIBM_END(__libm_error_region)
 670
 671 .type __libm_error_support#, @function
 672 .global __libm_error_support#