sysdeps/ia64/fpu/e_fmodf.S

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