sysdeps/ia64/fpu/e_fmod.S

   1 .file "fmod.s"
   2
   3 // Copyright (C) 2000, 2001, Intel Corporation
   4 // All rights reserved.
   5 //
   6 // Contributed 2/2/2000 by John Harrison, Cristina Iordache, Ted Kubaska,
   7 // Bob Norin, Shane Story, and Ping Tak Peter Tang of the Computational
   8 // Software Lab, Intel Corporation.
   9 //
  10 // Redistribution and use in source and binary forms, with or without
  11 // modification, are permitted provided that the following conditions are
  12 // met:
  13 //
  14 // * Redistributions of source code must retain the above copyright
  15 // notice, this list of conditions and the following disclaimer.
  16 //
  17 // * Redistributions in binary form must reproduce the above copyright
  18 // notice, this list of conditions and the following disclaimer in the
  19 // documentation and/or other materials provided with the distribution.
  20 //
  21 // * The name of Intel Corporation may not be used to endorse or promote
  22 // products derived from this software without specific prior written
  23 // permission.
  24 //
  25 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  26 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  27 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  28 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
  29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  33 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
  34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  36 //
  37 // Intel Corporation is the author of this code, and requests that all
  38 // problem reports or change requests be submitted to it directly at
  39 // http://developer.intel.com/opensource.
  40 //
  41 // History
  42 //====================================================================
  43 // 2/02/00  Initial version
  44 // 3/02/00  New Algorithm
  45 // 4/04/00  Unwind support added
  46 // 8/15/00  Bundle added after call to __libm_error_support to properly
  47 //          set [the previously overwritten] GR_Parameter_RESULT.
  48 //11/28/00  Set FR_Y to f9
  49 //
  50 // API
  51 //====================================================================
  52 // double fmod(double,double);
  53 //
  54 // Overview of operation
  55 //====================================================================
  56 //  fmod(a,b)=a-i*b,
  57 //  where i is an integer such that, if b!=0,
  58 //  |i|<|a/b| and |a/b-i|<1
  59 //
  60 // Algorithm
  61 //====================================================================
  62 // a). if |a|<|b|, return a
  63 // b). get quotient and reciprocal overestimates accurate to
  64 //     33 bits (q2,y2)
  65 // c). if the exponent difference (exponent(a)-exponent(b))
  66 //     is less than 32, truncate quotient to integer and
  67 //     finish in one iteration
  68 // d). if exponent(a)-exponent(b)>=32 (q2>=2^32)
  69 //     round quotient estimate to single precision (k=RN(q2)),
  70 //     calculate partial remainder (a'=a-k*b),
  71 //     get quotient estimate (a'*y2), and repeat from c).
  72 //
  73 // Special cases
  74 //====================================================================
  75 // b=+/-0: return NaN, call libm_error_support
  76 // a=+/-Inf, a=NaN or b=NaN: return NaN
  77 //
  78 // Registers used
  79 //====================================================================
  80 // Predicate registers: p6-p11
  81 // General registers:   r2,r29,r32 (ar.pfs), r33-r39
  82 // Floating point registers: f6-f15
  83
  84 #include "libm_support.h"
  85
  86 .section .text
  87
  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 .proc fmod#
 105 .align 32
 106 .global fmod#
 107 .align 32
 108
 109 fmod:
 110 #ifdef _LIBC
 111 .global __ieee754_fmod
 112 .type __ieee754_fmod,@function
 113 __ieee754_fmod:
 114 #endif
 115 // inputs in f8, f9
 116 // result in f8
 117
 118 { .mfi
 119   alloc r32=ar.pfs,1,4,4,0
 120   // f6=|a|
 121   fmerge.s f6=f0,f8
 122   mov r2 = 0x0ffdd
 123 }
 124   {.mfi
 125   nop.m 0
 126   // f7=|b|
 127   fmerge.s f7=f0,f9
 128   nop.i 0;;
 129 }
 130
 131 { .mfi
 132   setf.exp f11 = r2
 133   // (1) y0
 134   frcpa.s1 f10,p6=f6,f7
 135   nop.i 0
 136 }
 137
 138 // Y +-NAN, +-inf, +-0?     p7
 139 { .mfi
 140       nop.m 999
 141 (p0)  fclass.m.unc  p7,p0 = f9, 0xe7
 142       nop.i 999;;
 143 }
 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
 150 { .mfi
 151       nop.m 999
 152 (p0)  fclass.m.unc  p9,p0 = f8, 0xe3
 153       nop.i 999
 154 }
 155
 156 // |x| < |y|? Return x p8
 157 { .mfi
 158       nop.m 999
 159 (p0)  fcmp.lt.unc.s1 p8,p0 = f6,f7
 160       nop.i 999 ;;
 161 }
 162
 163 { .mfi
 164   nop.m 0
 165   // normalize y (if |x|<|y|)
 166   (p8) fma.s0 f9=f9,f1,f0
 167   nop.i 0;;
 168 }
 169
 170   { .mfi
 171   mov r2=0x1001f
 172   // (2) q0=a*y0
 173   (p6) fma.s1 f13=f6,f10,f0
 174   nop.i 0
 175 }
 176 { .mfi
 177   nop.m 0
 178   // (3) e0 = 1 - b * y0
 179   (p6) fnma.s1 f12=f7,f10,f1
 180   nop.i 0;;
 181 }
 182
 183   {.mfi
 184   nop.m 0
 185   // normalize x (if |x|<|y|)
 186   (p8) fma.d.s0 f8=f8,f1,f0
 187   nop.i 0
 188 }
 189 {.bbb
 190   (p9) br.cond.spnt L(FMOD_X_NAN_INF)
 191   (p7) br.cond.spnt L(FMOD_Y_NAN_INF_ZERO)
 192   // if |x|<|y|, return
 193   (p8) br.ret.spnt    b0;;
 194 }
 195
 196   {.mfi
 197   nop.m 0
 198   // normalize x
 199   fma.s0 f6=f6,f1,f0
 200   nop.i 0
 201 }
 202 {.mfi
 203   nop.m 0
 204   // normalize y
 205   fma.s0 f7=f7,f1,f0
 206   nop.i 0;;
 207 }
 208
 209   {.mfi
 210   // f15=2^32
 211   setf.exp f15=r2
 212   // (4) q1=q0+e0*q0
 213   (p6) fma.s1 f13=f12,f13,f13
 214   nop.i 0
 215 }
 216 { .mfi
 217   nop.m 0
 218   // (5) e1 = e0 * e0 + 2^-34
 219   (p6) fma.s1 f14=f12,f12,f11
 220   nop.i 0;;
 221 }
 222 {.mlx
 223   nop.m 0
 224   movl r2=0x33a00000;;
 225 }
 226 { .mfi
 227   nop.m 0
 228   // (6) y1 = y0 + e0 * y0
 229   (p6) fma.s1 f10=f12,f10,f10
 230   nop.i 0;;
 231 }
 232 {.mfi
 233   // set f12=1.25*2^{-24}
 234   setf.s f12=r2
 235   // (7) q2=q1+e1*q1
 236   (p6) fma.s1 f13=f13,f14,f13
 237   nop.i 0;;
 238 }
 239 {.mfi
 240   nop.m 0
 241   fmerge.s f9=f8,f9
 242   nop.i 0
 243 }
 244 { .mfi
 245   nop.m 0
 246   // (8) y2 = y1 + e1 * y1
 247   (p6) fma.s1 f10=f14,f10,f10
 248   // set p6=0, p10=0
 249   cmp.ne.and p6,p10=r0,r0;;
 250 }
 251
 252 .align 32
 253 L(loop53):
 254   {.mfi
 255   nop.m 0
 256   // compare q2, 2^32
 257   fcmp.lt.unc.s1 p8,p7=f13,f15
 258   nop.i 0
 259 }
 260   {.mfi
 261   nop.m 0
 262   // will truncate quotient to integer, if exponent<32 (in advance)
 263   fcvt.fx.trunc.s1 f11=f13
 264   nop.i 0;;
 265 }
 266   {.mfi
 267   nop.m 0
 268   // if exponent>32, round quotient to single precision (perform in advance)
 269   fma.s.s1 f13=f13,f1,f0
 270   nop.i 0;;
 271 }
 272   {.mfi
 273   nop.m 0
 274   // set f12=sgn(a)
 275   (p8) fmerge.s f12=f8,f1
 276   nop.i 0
 277 }
 278   {.mfi
 279   nop.m 0
 280   // normalize truncated quotient
 281   (p8) fcvt.xf f13=f11
 282   nop.i 0;;
 283 }
 284   { .mfi
 285   nop.m 0
 286   // calculate remainder (assuming f13=RZ(Q))
 287   (p7) fnma.s1 f14=f13,f7,f6
 288   nop.i 0
 289 }
 290   {.mfi
 291   nop.m 0
 292   // also if exponent>32, round quotient to single precision
 293   // and subtract 1 ulp: q=q-q*(1.25*2^{-24})
 294   (p7) fnma.s.s1 f11=f13,f12,f13
 295   nop.i 0;;
 296 }
 297
 298   {.mfi
 299   nop.m 0
 300   // (p8) calculate remainder (82-bit format)
 301   (p8) fnma.s1 f11=f13,f7,f6
 302   nop.i 0
 303 }
 304   {.mfi
 305   nop.m 0
 306   // (p7) calculate remainder (assuming f11=RZ(Q))
 307   (p7) fnma.s1 f6=f11,f7,f6
 308   nop.i 0;;
 309 }
 310
 311
 312   {.mfi
 313   nop.m 0
 314   // Final iteration (p8): is f6 the correct remainder (quotient was not overestimated) ?
 315   (p8) fcmp.lt.unc.s1 p6,p10=f11,f0
 316   nop.i 0;;
 317 }
 318   {.mfi
 319   nop.m 0
 320   // get new quotient estimation: a'*y2
 321   (p7) fma.s1 f13=f14,f10,f0
 322   nop.i 0
 323 }
 324   {.mfb
 325   nop.m 0
 326   // was f14=RZ(Q) ? (then new remainder f14>=0)
 327   (p7) fcmp.lt.unc.s1 p7,p9=f14,f0
 328   nop.b 0;;
 329 }
 330
 331
 332 .pred.rel "mutex",p6,p10
 333   {.mfb
 334   nop.m 0
 335   // add b to estimated remainder (to cover the case when the quotient was overestimated)
 336   // also set correct sign by using f9=|b|*sgn(a), f12=sgn(a)
 337   (p6) fma.d.s0 f8=f11,f12,f9
 338   nop.b 0
 339 }
 340   {.mfb
 341   nop.m 0
 342   // calculate remainder (single precision)
 343   // set correct sign of result before returning
 344   (p10) fma.d.s0 f8=f11,f12,f0
 345   (p8) br.ret.sptk b0;;
 346 }
 347   {.mfi
 348   nop.m 0
 349   // if f13!=RZ(Q), get alternative quotient estimation: a''*y2
 350   (p7) fma.s1 f13=f6,f10,f0
 351   nop.i 0
 352 }
 353   {.mfb
 354   nop.m 0
 355   // if f14 was RZ(Q), set remainder to f14
 356   (p9) mov f6=f14
 357   br.cond.sptk L(loop53);;
 358 }
 359
 360
 361
 362 L(FMOD_X_NAN_INF):
 363
 364 // Y zero ?
 365 {.mfi
 366   nop.m 0
 367   fma.s1 f10=f9,f1,f0
 368   nop.i 0;;
 369 }
 370 {.mfi
 371  nop.m 0
 372  fcmp.eq.unc.s1 p11,p0=f10,f0
 373  nop.i 0;;
 374 }
 375 {.mib
 376   nop.m 0
 377   nop.i 0
 378   // if Y zero
 379   (p11) br.cond.spnt L(FMOD_Y_ZERO);;
 380 }
 381
 382 // X infinity? Return QNAN indefinite
 383 { .mfi
 384       nop.m 999
 385 (p0)  fclass.m.unc  p8,p9 = f8, 0x23
 386       nop.i 999;;
 387 }
 388 // Y NaN ?
 389 {.mfi
 390          nop.m 999
 391 (p8) fclass.m p9,p8=f9,0xc3
 392          nop.i 0;;
 393 }
 394 {.mfi
 395           nop.m 999
 396 (p8)  frcpa.s0 f8,p0 = f8,f8
 397       nop.i 0
 398 }
 399 { .mfi
 400       nop.m 999
 401         // also set Denormal flag if necessary
 402 (p8)  fma.s0 f9=f9,f1,f0
 403       nop.i 999 ;;
 404 }
 405
 406 { .mfb
 407       nop.m 999
 408 (p8)  fma.d f8=f8,f1,f0
 409           nop.b 999 ;;
 410 }
 411
 412 { .mfb
 413       nop.m 999
 414 (p9)  frcpa.s0 f8,p7=f8,f9
 415       br.ret.sptk   b0 ;;
 416 }
 417
 418
 419 L(FMOD_Y_NAN_INF_ZERO):
 420
 421 // Y INF
 422 { .mfi
 423       nop.m 999
 424 (p0)  fclass.m.unc  p7,p0 = f9, 0x23
 425       nop.i 999 ;;
 426 }
 427
 428 { .mfb
 429       nop.m 999
 430 (p7)  fma.d f8=f8,f1,f0
 431 (p7)  br.ret.spnt    b0 ;;
 432 }
 433
 434 // Y NAN?
 435 { .mfi
 436       nop.m 999
 437 (p0)  fclass.m.unc  p9,p0 = f9, 0xc3
 438       nop.i 999 ;;
 439 }
 440
 441 { .mfb
 442       nop.m 999
 443 (p9)  fma.d f8=f9,f1,f0
 444 (p9)  br.ret.spnt    b0 ;;
 445 }
 446
 447 L(FMOD_Y_ZERO):
 448 // Y zero? Must be zero at this point
 449 // because it is the only choice left.
 450 // Return QNAN indefinite
 451
 452 {.mfi
 453   nop.m 0
 454   // set Invalid
 455   frcpa f12,p0=f0,f0
 456   nop.i 0
 457 }
 458 // X NAN?
 459 { .mfi
 460       nop.m 999
 461 (p0)  fclass.m.unc  p9,p10 = f8, 0xc3
 462       nop.i 999 ;;
 463 }
 464 { .mfi
 465       nop.m 999
 466 (p10)  fclass.nm  p9,p10 = f8, 0xff
 467       nop.i 999 ;;
 468 }
 469
 470 {.mfi
 471  nop.m 999
 472  (p9) frcpa f11,p7=f8,f0
 473  nop.i 0;;
 474 }
 475
 476 { .mfi
 477       nop.m 999
 478 (p10)  frcpa         f11,p7 = f9,f9
 479 (p0)  mov        GR_Parameter_TAG = 121 ;;
 480 }
 481
 482 { .mfi
 483       nop.m 999
 484 (p0)  fmerge.s      f10 = f8, f8
 485       nop.i 999
 486 }
 487
 488 { .mfb
 489       nop.m 999
 490 (p0)  fma.d f8=f11,f1,f0
 491 (p0)  br.sptk __libm_error_region;;
 492 }
 493
 494 .endp fmod
 495 ASM_SIZE_DIRECTIVE(fmod)
 496 ASM_SIZE_DIRECTIVE(__ieee754_fmod)
 497
 498 .proc __libm_error_region
 499 __libm_error_region:
 500 .prologue
 501 { .mfi
 502         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 503         nop.f 0
 504 .save   ar.pfs,GR_SAVE_PFS
 505         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 506 }
 507 { .mfi
 508 .fframe 64
 509         add sp=-64,sp                           // Create new stack
 510         nop.f 0
 511         mov GR_SAVE_GP=gp                       // Save gp
 512 };;
 513 { .mmi
 514         stfd [GR_Parameter_Y] = FR_Y,16         // Save Parameter 2 on stack
 515         add GR_Parameter_X = 16,sp              // Parameter 1 address
 516 .save   b0, GR_SAVE_B0
 517         mov GR_SAVE_B0=b0                       // Save b0
 518 };;
 519 .body
 520 { .mib
 521         stfd [GR_Parameter_X] = FR_X            // Store Parameter 1 on stack
 522         add   GR_Parameter_RESULT = 0,GR_Parameter_Y
 523         nop.b 0                                 // Parameter 3 address
 524 }
 525 { .mib
 526         stfd [GR_Parameter_Y] = FR_RESULT      // Store Parameter 3 on stack
 527         add   GR_Parameter_Y = -16,GR_Parameter_Y
 528         br.call.sptk b0=__libm_error_support#  // Call error handling function
 529 };;
 530 { .mmi
 531         nop.m 0
 532         nop.m 0
 533         add   GR_Parameter_RESULT = 48,sp
 534 };;
 535 { .mmi
 536         ldfd  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 537 .restore sp
 538         add   sp = 64,sp                       // Restore stack pointer
 539         mov   b0 = GR_SAVE_B0                  // Restore return address
 540 };;
 541 { .mib
 542         mov   gp = GR_SAVE_GP                  // Restore gp
 543         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 544         br.ret.sptk     b0                     // Return
 545 };;
 546
 547 .endp __libm_error_region
 548 ASM_SIZE_DIRECTIVE(__libm_error_region)
 549
 550 .type   __libm_error_support#,@function
 551 .global __libm_error_support#