sysdeps/ia64/fpu/e_fmodl.S

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