sysdeps/ia64/fpu/e_remainder.S

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