sysdeps/ia64/fpu/e_remainder.S

   1   .file "remainder.asm"
   2 // Copyright (C) 2000, 2001, Intel Corporation
   3 // All rights reserved.
   4 //
   5 // Contributed 2/2/2000 by John Harrison, Cristina Iordache, Ted Kubaska, Bob Norin,
   6 // Shane Story, and Ping Tak Peter Tang of the Computational Software Lab,
   7 // 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://developer.intel.com/opensource.
  39 //
  40 // History
  41 //====================================================================
  42 // 2/02/00  Initial version
  43 // 3/02/00  New Algorithm
  44 // 4/04/00  Unwind support added
  45 // 7/21/00  Fixed quotient=2^{24*m+23}*1.q1...q23 1 bug
  46 // 8/15/00  Bundle added after call to __libm_error_support to properly
  47 //          set [the previously overwritten] GR_Parameter_RESULT.
  48 //11/29/00  Set FR_Y to f9
  49 //
  50 // API
  51 //====================================================================
  52 // double remainder(double,double);
  53 //
  54 // Overview of operation
  55 //====================================================================
  56 //  remainder(a,b)=a-i*b,
  57 //  where i is an integer such that, if b!=0 and a is finite,
  58 //  |a/b-i|<=1/2. If |a/b-i|=1/2, i is even.
  59 //
  60 // Algorithm
  61 //====================================================================
  62 // a). eliminate special cases
  63 // b). if |a/b|<0.25 (first quotient estimate), return a
  64 // c). use single precision divide algorithm to get quotient q
  65 //     rounded to 24 bits of precision
  66 // d). calculate partial remainders (using both q and q-ulp);
  67 //     select one and RZ(a/b) based on the sign of |a|-|b|*q
  68 // e). if the exponent difference (exponent(a)-exponent(b))
  69 //     is less than 24 (quotient estimate<2^{24}-2), use RZ(a/b)
  70 //     and sticky bits to round to integer; exit loop and
  71 //     calculate final remainder
  72 // f). if exponent(a)-exponent(b)>=24, select new value of a as
  73 //     the partial remainder calculated using RZ(a/b);
  74 //     repeat from c).
  75 //
  76 // Special cases
  77 //====================================================================
  78 // a=+/- Inf, or b=+/-0: return NaN, call libm_error_support
  79 // a=NaN or b=NaN: return NaN
  80
  81 #include "libm_support.h"
  82
  83 // Registers used
  84 //====================================================================
  85 // Predicate registers: p6-p14
  86 // General registers:   r2,r3,r28,r29,r32 (ar.pfs), r33-r39
  87 // Floating point registers: f6-f15,f32
  88
  89   .section .text
  90
  91 GR_SAVE_B0                    = r33
  92 GR_SAVE_PFS                   = r34
  93 GR_SAVE_GP                    = r35
  94 GR_SAVE_SP                    = r36
  95
  96 GR_Parameter_X                = r37
  97 GR_Parameter_Y                = r38
  98 GR_Parameter_RESULT           = r39
  99 GR_Parameter_TAG              = r40
 100
 101 FR_X             = f10
 102 FR_Y             = f9
 103 FR_RESULT        = f8
 104
 105
 106
 107   .proc  remainder#
 108   .align 32
 109   .global remainder#
 110   .align 32
 111
 112 remainder:
 113 #ifdef _LIBC
 114 .global __remainder
 115 .type __remainder,@function
 116 __remainder:
 117 #endif
 118 // inputs in f8, f9
 119 // result in f8
 120
 121 { .mfi
 122   alloc r32=ar.pfs,1,4,4,0
 123   // f13=|a|
 124   fmerge.s f13=f0,f8
 125   nop.i 0
 126 }
 127   {.mfi
 128   nop.m 0
 129   // f14=|b|
 130   fmerge.s f14=f0,f9
 131   nop.i 0;;
 132 }
 133  {.mlx
 134   mov r28=0x2ffdd
 135   // r2=2^{23}
 136   movl r3=0x4b000000;;
 137 }
 138
 139 // Y +-NAN, +-inf, +-0?     p11
 140 { .mfi
 141           setf.exp f32=r28
 142 (p0)  fclass.m.unc  p11,p0 = f9, 0xe7
 143       nop.i 999
 144 }
 145 // qnan snan inf norm     unorm 0 -+
 146 // 1    1    1   0        0     0 11
 147 // e                      3
 148 // X +-NAN, +-inf, ?        p9
 149 { .mfi
 150       nop.m 999
 151 (p0)  fclass.m.unc  p9,p0 = f8, 0xe3
 152       nop.i 999;;
 153 }
 154
 155 {.mfi
 156   nop.m 0
 157   mov f12=f0
 158   nop.i 0
 159 }
 160 { .mfi
 161   // set p7=1
 162   cmp.eq.unc p7,p0=r0,r0
 163   // Step (1)
 164   // y0 = 1 / b in f10
 165   frcpa.s1 f10,p6=f13,f14
 166   nop.i 0;;
 167 }
 168
 169 {.bbb
 170   (p9) br.cond.spnt L(FREM_X_NAN_INF)
 171   (p11) br.cond.spnt L(FREM_Y_NAN_INF_ZERO)
 172   nop.b 0
 173 }  {.mfi
 174    nop.m 0
 175    // set D flag if a (f8) is denormal
 176    fnma.s0 f6=f8,f1,f8
 177    nop.i 0;;
 178 }
 179
 180
 181 L(remloop24):
 182   { .mfi
 183   nop.m 0
 184   // Step (2)
 185   // q0 = a * y0 in f12
 186   (p6) fma.s1 f12=f13,f10,f0
 187   nop.i 0
 188 } { .mfi
 189   nop.m 0
 190   // Step (3)
 191   // e0 = 1 - b * y0 in f7
 192   (p6) fnma.s1 f7=f14,f10,f1
 193   nop.i 0;;
 194 }  {.mlx
 195   nop.m 0
 196   // r2=1.25*2^{-24}
 197   movl r2=0x33a00000;;
 198 }
 199
 200 {.mfi
 201   nop.m 0
 202   // q1=q0*(1+e0)
 203   fma.s1 f15=f12,f7,f12
 204   nop.i 0
 205 }
 206 { .mfi
 207   nop.m 0
 208   // Step (4)
 209   // e1 = e0 * e0 + E in f7
 210   (p6) fma.s1 f7=f7,f7,f32
 211   nop.i 0;;
 212 }
 213  {.mii
 214   (p7) getf.exp r29=f12
 215   (p7) mov r28=0xfffd
 216   nop.i 0;;
 217 }
 218  { .mfi
 219   // f12=2^{23}
 220   setf.s f12=r3
 221   // Step (5)
 222   // q2 = q1 + e1 * q1 in f11
 223   (p6) fma.s.s1 f11=f7,f15,f15
 224   nop.i 0
 225 } { .mfi
 226    nop.m 0
 227   // Step (6)
 228   // q2 = q1 + e1 * q1 in f6
 229   (p6) fma.s1 f6=f7,f15,f15
 230   nop.i 0;;
 231 }
 232
 233  {.mmi
 234   // f15=1.25*2^{-24}
 235   setf.s f15=r2
 236   // q<1/4 ? (i.e. expon< -2)
 237   (p7) cmp.gt p7,p0=r28,r29
 238   nop.i 0;;
 239 }
 240
 241 {.mfb
 242   // r29= -32+bias
 243   mov r29=0xffdf
 244  // if |a/b|<1/4, set D flag before returning
 245  (p7) fma.d.s0 f9=f9,f0,f8
 246   nop.b 0;;
 247 }
 248  {.mfb
 249  nop.m 0
 250  // can be combined with bundle above if sign of 0 or
 251  // FTZ enabled are not important
 252  (p7) fmerge.s f8=f8,f9
 253  // return if |a|<4*|b| (estimated quotient < 1/4)
 254  (p7) br.ret.spnt b0;;
 255 }
 256   {.mfi
 257   // f7=2^{-32}
 258   setf.exp f7=r29
 259   // set f8 to current a value | sign
 260   fmerge.s f8=f8,f13
 261   nop.i 0;;
 262 }
 263
 264
 265   {.mfi
 266   getf.exp r28=f6
 267   // last step ? (q<2^{23})
 268   fcmp.lt.unc.s1 p0,p12=f6,f12
 269   nop.i 0;;
 270 }
 271   {.mfi
 272   nop.m 0
 273   // r=a-b*q
 274   fnma.s1 f6=f14,f11,f13
 275   nop.i 0
 276 } {.mfi
 277   // r2=23+bias
 278   mov r2=0xffff+23
 279   // q'=q-q*(1.25*2^{-24})   (q'=q-ulp)
 280   fnma.s.s1 f15=f11,f15,f11
 281   nop.i 0;;
 282 }
 283   {.mmi
 284   nop.m 0
 285   cmp.eq p11,p14=r2,r28
 286   nop.i 0;;
 287 }
 288
 289 .pred.rel "mutex",p11,p14
 290   {.mfi
 291   nop.m 0
 292   // if exp_q=2^23, then r=a-b*2^{23}
 293   (p11) fnma.s1 f13=f12,f14,f13
 294   nop.i 0
 295 }
 296 {.mfi
 297   nop.m 0
 298   // r2=a-b*q'
 299   (p14) fnma.s1 f13=f14,f15,f13
 300   nop.i 0;;
 301 }
 302   {.mfi
 303   nop.m 0
 304   // r>0 iff q=RZ(a/b) and inexact
 305   fcmp.gt.unc.s1 p8,p0=f6,f0
 306   nop.i 0
 307 } {.mfi
 308   nop.m 0
 309   // r<0 iff q'=RZ(a/b) and inexact
 310   (p14) fcmp.lt.unc.s1 p9,p10=f6,f0
 311   nop.i 0;;
 312 }
 313
 314 .pred.rel "mutex",p8,p9
 315   {.mfi
 316    nop.m 0
 317   // (p8) Q=q+(last iteration ? sticky bits:0)
 318   // i.e. Q=q+q*x  (x=2^{-32} or 0)
 319   (p8) fma.s1 f11=f11,f7,f11
 320   nop.i 0
 321 } {.mfi
 322   nop.m 0
 323   // (p9) Q=q'+(last iteration ? sticky bits:0)
 324   // i.e. Q=q'+q'*x  (x=2^{-32} or 0)
 325   (p9) fma.s1 f11=f15,f7,f15
 326   nop.i 0;;
 327 }
 328
 329   {.mfb
 330   nop.m 0
 331   //  (p9) set r=r2 (new a, if not last iteration)
 332   // (p10) new a =r
 333   (p10) mov f13=f6
 334   (p12) br.cond.sptk L(remloop24);;
 335 }
 336
 337 // last iteration
 338   {.mfi
 339   nop.m 0
 340   // set f9=|b|*sgn(a)
 341   fmerge.s f9=f8,f9
 342   nop.i 0
 343 }
 344   {.mfi
 345   nop.m 0
 346   // round to integer
 347   fcvt.fx.s1 f11=f11
 348   nop.i 0;;
 349 }
 350   {.mfi
 351   nop.m 0
 352   // save sign of a
 353   fmerge.s f7=f8,f8
 354   nop.i 0
 355 } {.mfi
 356   nop.m 0
 357   // normalize
 358   fcvt.xf f11=f11
 359   nop.i 0;;
 360 }
 361   {.mfi
 362   nop.m 0
 363   // This can be removed if sign of 0 is not important
 364   // get remainder using sf1
 365   fnma.d.s1 f12=f9,f11,f8
 366   nop.i 0
 367 }
 368   {.mfi
 369   nop.m 0
 370   // get remainder
 371   fnma.d.s0 f8=f9,f11,f8
 372   nop.i 0;;
 373 }
 374   {.mfi
 375   nop.m 0
 376   // f12=0?
 377   // This can be removed if sign of 0 is not important
 378   fcmp.eq.unc.s1 p8,p0=f12,f0
 379   nop.i 0;;
 380 }
 381   {.mfb
 382   nop.m 0
 383   // if f8=0, set sign correctly
 384   // This can be removed if sign of 0 is not important
 385   (p8) fmerge.s f8=f7,f8
 386   // return
 387   br.ret.sptk b0;;
 388 }
 389
 390
 391 L(FREM_X_NAN_INF):
 392
 393 // Y zero ?
 394 {.mfi
 395   nop.m 0
 396   fma.s1 f10=f9,f1,f0
 397   nop.i 0;;
 398 }
 399 {.mfi
 400  nop.m 0
 401  fcmp.eq.unc.s1 p11,p0=f10,f0
 402  nop.i 0;;
 403 }
 404 {.mib
 405   nop.m 0
 406   nop.i 0
 407   // if Y zero
 408   (p11) br.cond.spnt L(FREM_Y_ZERO);;
 409 }
 410
 411 // X infinity? Return QNAN indefinite
 412 { .mfi
 413       nop.m 999
 414 (p0)  fclass.m.unc  p8,p0 = f8, 0x23
 415       nop.i 999
 416 }
 417 // X infinity? Return QNAN indefinite
 418 { .mfi
 419       nop.m 999
 420 (p0)  fclass.m.unc  p11,p0 = f8, 0x23
 421       nop.i 999;;
 422 }
 423 // Y NaN ?
 424 {.mfi
 425          nop.m 999
 426 (p8) fclass.m.unc p0,p8=f9,0xc3
 427          nop.i 0;;
 428 }
 429 {.mfi
 430         nop.m 999
 431         // also set Denormal flag if necessary
 432 (p8) fma.s0 f9=f9,f1,f0
 433     nop.i 0
 434 }
 435 { .mfi
 436       nop.m 999
 437 (p8)  frcpa.s0 f8,p7 = f8,f8
 438       nop.i 999 ;;
 439 }
 440
 441 {.mfi
 442       nop.m 999
 443 (p11) mov f10=f8
 444           nop.i 0
 445 }
 446 { .mfi
 447       nop.m 999
 448 (p8) fma.d f8=f8,f1,f0
 449           nop.i 0 ;;
 450 }
 451
 452 { .mfb
 453       nop.m 999
 454       frcpa.s0 f8,p7=f8,f9
 455           (p11) br.cond.spnt L(EXP_ERROR_RETURN);;
 456 }
 457 { .mib
 458         nop.m 0
 459         nop.i 0
 460         br.ret.spnt    b0 ;;
 461 }
 462
 463
 464 L(FREM_Y_NAN_INF_ZERO):
 465
 466 // Y INF
 467 { .mfi
 468       nop.m 999
 469 (p0)  fclass.m.unc  p7,p0 = f9, 0x23
 470       nop.i 999 ;;
 471 }
 472
 473 { .mfb
 474       nop.m 999
 475 (p7)  fma.d f8=f8,f1,f0
 476 (p7)  br.ret.spnt    b0 ;;
 477 }
 478
 479 // Y NAN?
 480 { .mfi
 481       nop.m 999
 482 (p0)  fclass.m.unc  p9,p0 = f9, 0xc3
 483       nop.i 999 ;;
 484 }
 485
 486 { .mfb
 487       nop.m 999
 488 (p9)  fma.d f8=f9,f1,f0
 489 (p9)  br.ret.spnt    b0 ;;
 490 }
 491
 492 L(FREM_Y_ZERO):
 493 // Y zero? Must be zero at this point
 494 // because it is the only choice left.
 495 // Return QNAN indefinite
 496
 497 // X NAN?
 498 { .mfi
 499       nop.m 999
 500 (p0)  fclass.m.unc  p9,p10 = f8, 0xc3
 501       nop.i 999 ;;
 502 }
 503 { .mfi
 504       nop.m 999
 505 (p10)  fclass.nm  p9,p10 = f8, 0xff
 506       nop.i 999 ;;
 507 }
 508
 509 {.mfi
 510  nop.m 999
 511  (p9) frcpa f11,p7=f8,f0
 512  nop.i 0;;
 513 }
 514
 515 { .mfi
 516       nop.m 999
 517 (p10)  frcpa         f11,p7 = f0,f0
 518           nop.i 999;;
 519 }
 520
 521 { .mfi
 522       nop.m 999
 523 (p0)  fmerge.s      f10 = f8, f8
 524       nop.i 999
 525 }
 526
 527 { .mfi
 528       nop.m 999
 529 (p0)  fma.d f8=f11,f1,f0
 530       nop.i 999
 531 }
 532
 533
 534 L(EXP_ERROR_RETURN):
 535
 536 { .mib
 537 (p0)  mov   GR_Parameter_TAG = 124
 538           nop.i 999
 539 (p0)  br.sptk __libm_error_region;;
 540 }
 541
 542 .endp remainder
 543 ASM_SIZE_DIRECTIVE(remainder)
 544 #ifdef _LIBC
 545 ASM_SIZE_DIRECTIVE(__remainder)
 546 #endif
 547
 548
 549
 550 .proc __libm_error_region
 551 __libm_error_region:
 552 .prologue
 553 { .mfi
 554         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 555         nop.f 0
 556 .save   ar.pfs,GR_SAVE_PFS
 557         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 558 }
 559 { .mfi
 560 .fframe 64
 561         add sp=-64,sp                           // Create new stack
 562         nop.f 0
 563         mov GR_SAVE_GP=gp                       // Save gp
 564 };;
 565 { .mmi
 566         stfd [GR_Parameter_Y] = FR_Y,16         // Save Parameter 2 on stack
 567         add GR_Parameter_X = 16,sp              // Parameter 1 address
 568 .save   b0, GR_SAVE_B0
 569         mov GR_SAVE_B0=b0                       // Save b0
 570 };;
 571 .body
 572 { .mib
 573         stfd [GR_Parameter_X] = FR_X            // Store Parameter 1 on stack
 574         add   GR_Parameter_RESULT = 0,GR_Parameter_Y
 575         nop.b 0                                 // Parameter 3 address
 576 }
 577 { .mib
 578         stfd [GR_Parameter_Y] = FR_RESULT      // Store Parameter 3 on stack
 579         add   GR_Parameter_Y = -16,GR_Parameter_Y
 580         br.call.sptk b0=__libm_error_support#  // Call error handling function
 581 };;
 582 { .mmi
 583         nop.m 0
 584         nop.m 0
 585         add   GR_Parameter_RESULT = 48,sp
 586 };;
 587 { .mmi
 588         ldfd  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 589 .restore sp
 590         add   sp = 64,sp                       // Restore stack pointer
 591         mov   b0 = GR_SAVE_B0                  // Restore return address
 592 };;
 593 { .mib
 594         mov   gp = GR_SAVE_GP                  // Restore gp
 595         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 596         br.ret.sptk     b0                     // Return
 597 };;
 598
 599 .endp __libm_error_region
 600 ASM_SIZE_DIRECTIVE(__libm_error_region)
 601
 602
 603
 604 .type   __libm_error_support#,@function
 605 .global __libm_error_support#