sysdeps/ia64/fpu/e_fmodf.S

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