sysdeps/ia64/fpu/e_scalbl.S

   1 .file "scalbl.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 // 01/26/01 Scalb completely reworked and now standalone version
  44 // 05/20/02 Cleaned up namespace and sf0 syntax
  45 // 02/10/03 Reordered header: .section, .global, .proc, .align
  46 // 08/06/03 Improved performance
  47 //
  48 // API
  49 //==============================================================
  50 // long double = scalbl  (long double x, long double n)
  51 // input  floating point f8 and floating point f9
  52 // output floating point f8
  53 //
  54 // int_type = 0 if int is 32 bits
  55 // int_type = 1 if int is 64 bits
  56 //
  57 // Returns x* 2**n using an fma and detects overflow
  58 // and underflow.
  59 //
  60 //
  61 // Strategy:
  62 //  Compute biased exponent of result exp_Result = N + exp_X
  63 //  Break into ranges:
  64 //   exp_Result > 0x13ffe                 -> Certain overflow
  65 //   exp_Result = 0x13ffe                 -> Possible overflow
  66 //   0x0c001 <= exp_Result < 0x13ffe      -> No over/underflow (main path)
  67 //   0x0c001 - 63 <= exp_Result < 0x0c001 -> Possible underflow
  68 //   exp_Result < 0x0c001 - 63            -> Certain underflow
  69
  70 FR_Big         = f6
  71 FR_NBig        = f7
  72 FR_Floating_X  = f8
  73 FR_Result      = f8
  74 FR_Floating_N  = f9
  75 FR_Result2     = f9
  76 FR_Result3     = f10
  77 FR_Norm_X      = f11
  78 FR_Two_N       = f12
  79 FR_N_float_int = f13
  80 FR_Norm_N      = f14
  81
  82 GR_neg_ov_limit= r14
  83 GR_big_exp     = r14
  84 GR_N_Biased    = r15
  85 GR_Big         = r16
  86 GR_exp_Result  = r18
  87 GR_pos_ov_limit= r19
  88 GR_exp_sure_ou = r19
  89 GR_Bias        = r20
  90 GR_N_as_int    = r21
  91 GR_signexp_X   = r22
  92 GR_exp_X       = r23
  93 GR_exp_mask    = r24
  94 GR_max_exp     = r25
  95 GR_min_exp     = r26
  96 GR_min_den_exp = r27
  97 GR_Scratch     = r28
  98 GR_signexp_N   = r29
  99 GR_exp_N       = r30
 100
 101 GR_SAVE_B0          = r32
 102 GR_SAVE_GP          = r33
 103 GR_SAVE_PFS         = r34
 104 GR_Parameter_X      = r35
 105 GR_Parameter_Y      = r36
 106 GR_Parameter_RESULT = r37
 107 GR_Tag              = r38
 108
 109 .section .text
 110 GLOBAL_IEEE754_ENTRY(scalbl)
 111
 112 //
 113 //   Is x NAN, INF, ZERO, +-?
 114 //   Build the exponent Bias
 115 //
 116 {    .mfi
 117      getf.exp      GR_signexp_N = FR_Floating_N // Get signexp of n
 118      fclass.m      p6,p0 = FR_Floating_X, 0xe7  // @snan | @qnan | @inf | @zero
 119      mov           GR_Bias = 0x0ffff
 120 }
 121 {    .mfi
 122      mov           GR_Big = 35000      // If N this big then certain overflow
 123      fcvt.fx.trunc.s1   FR_N_float_int = FR_Floating_N // Get N in significand
 124      nop.i         0
 125 }
 126 ;;
 127
 128 {    .mfi
 129      getf.exp      GR_signexp_X = FR_Floating_X // Get signexp of x
 130      fclass.m      p7,p0 = FR_Floating_N, 0x0b  // Test for n=unorm
 131      nop.i         0
 132 }
 133 //
 134 //   Normalize n
 135 //
 136 {    .mfi
 137      mov           GR_exp_mask = 0x1ffff     // Exponent mask
 138      fnorm.s1      FR_Norm_N = FR_Floating_N
 139      nop.i         0
 140 }
 141 ;;
 142
 143 //
 144 //   Is n NAN, INF, ZERO, +-?
 145 //
 146 {    .mfi
 147      mov           GR_big_exp = 0x1003e      // Exponent at which n is integer
 148      fclass.m      p9,p0 = FR_Floating_N, 0xe7  // @snan | @qnan | @inf | @zero
 149      mov           GR_max_exp = 0x13ffe      // Exponent of maximum long double
 150 }
 151 //
 152 //   Normalize x
 153 //
 154 { .mfb
 155      nop.m         0
 156      fnorm.s1      FR_Norm_X = FR_Floating_X
 157 (p7) br.cond.spnt  SCALBL_N_UNORM             // Branch if n=unorm
 158 }
 159 ;;
 160
 161 SCALBL_COMMON1:
 162 // Main path continues.  Also return here from u=unorm path.
 163 //   Handle special cases if x = Nan, Inf, Zero
 164 { .mfb
 165      nop.m         0
 166      fcmp.lt.s1    p7,p0 = FR_Floating_N, f0  // Test N negative
 167 (p6) br.cond.spnt  SCALBL_NAN_INF_ZERO
 168 }
 169 ;;
 170
 171 //   Handle special cases if n = Nan, Inf, Zero
 172 {    .mfi
 173      getf.sig      GR_N_as_int = FR_N_float_int // Get n from significand
 174      fclass.m      p8,p0 = FR_Floating_X, 0x0b // Test for x=unorm
 175      mov           GR_exp_sure_ou = 0x1000e // Exp_N where x*2^N sure over/under
 176 }
 177 {    .mfb
 178      mov           GR_min_exp = 0x0c001      // Exponent of minimum long double
 179      fcvt.xf       FR_N_float_int = FR_N_float_int // Convert N to FP integer
 180 (p9) br.cond.spnt  SCALBL_NAN_INF_ZERO
 181 }
 182 ;;
 183
 184 {    .mmi
 185      and           GR_exp_N = GR_exp_mask, GR_signexp_N // Get exponent of N
 186 (p7) sub           GR_Big = r0, GR_Big          // Limit for N
 187      nop.i         0
 188 }
 189 ;;
 190
 191 {    .mib
 192      cmp.lt        p9,p0 = GR_exp_N, GR_big_exp // N possible non-integer?
 193      cmp.ge        p6,p0 = GR_exp_N, GR_exp_sure_ou // N certain over/under?
 194 (p8) br.cond.spnt  SCALBL_X_UNORM             // Branch if x=unorm
 195 }
 196 ;;
 197
 198 SCALBL_COMMON2:
 199 // Main path continues.  Also return here from x=unorm path.
 200 //   Create biased exponent for 2**N
 201 {    .mmi
 202 (p6) mov           GR_N_as_int = GR_Big      // Limit N
 203 ;;
 204      add           GR_N_Biased = GR_Bias,GR_N_as_int
 205      nop.i         0
 206 }
 207 ;;
 208
 209 {    .mfi
 210      setf.exp      FR_Two_N = GR_N_Biased               // Form 2**N
 211 (p9) fcmp.neq.unc.s1 p9,p0 = FR_Norm_N, FR_N_float_int  // Test if N an integer
 212      and           GR_exp_X = GR_exp_mask, GR_signexp_X // Get exponent of X
 213 }
 214 ;;
 215
 216 //
 217 //   Compute biased result exponent
 218 //   Branch if N is not an integer
 219 //
 220 {    .mib
 221      add           GR_exp_Result = GR_exp_X, GR_N_as_int
 222      mov           GR_min_den_exp = 0x0c001 - 63 // Exp of min denorm long dble
 223 (p9) br.cond.spnt  SCALBL_N_NOT_INT
 224 }
 225 ;;
 226
 227 //
 228 //   Raise Denormal operand flag with compare
 229 //   Do final operation
 230 //
 231 {    .mfi
 232      cmp.lt        p7,p6 = GR_exp_Result, GR_max_exp  // Test no overflow
 233      fcmp.ge.s0    p0,p11 = FR_Floating_X,FR_Floating_N  // Dummy to set denorm
 234      cmp.lt        p9,p0 = GR_exp_Result, GR_min_den_exp // Test sure underflow
 235 }
 236 {    .mfb
 237      nop.m         0
 238      fma.s0        FR_Result = FR_Two_N,FR_Norm_X,f0
 239 (p9) br.cond.spnt  SCALBL_UNDERFLOW           // Branch if certain underflow
 240 }
 241 ;;
 242
 243 {    .mib
 244 (p6) cmp.gt.unc    p6,p8 = GR_exp_Result, GR_max_exp  // Test sure overflow
 245 (p7) cmp.ge.unc    p7,p9 = GR_exp_Result, GR_min_exp  // Test no over/underflow
 246 (p7) br.ret.sptk   b0                         // Return from main path
 247 }
 248 ;;
 249
 250 {    .bbb
 251 (p6) br.cond.spnt  SCALBL_OVERFLOW            // Branch if certain overflow
 252 (p8) br.cond.spnt  SCALBL_POSSIBLE_OVERFLOW   // Branch if possible overflow
 253 (p9) br.cond.spnt  SCALBL_POSSIBLE_UNDERFLOW  // Branch if possible underflow
 254 }
 255 ;;
 256
 257 // Here if possible underflow.
 258 // Resulting exponent: 0x0c001-63 <= exp_Result < 0x0c001
 259 SCALBL_POSSIBLE_UNDERFLOW:
 260 //
 261 // Here if possible overflow.
 262 // Resulting exponent: 0x13ffe = exp_Result
 263 SCALBL_POSSIBLE_OVERFLOW:
 264
 265 //   Set up necessary status fields
 266 //
 267 //   S0 user supplied status
 268 //   S2 user supplied status + WRE + TD  (Overflows)
 269 //   S3 user supplied status + FZ + TD   (Underflows)
 270 //
 271 {    .mfi
 272      mov           GR_pos_ov_limit = 0x13fff // Exponent for positive overflow
 273      fsetc.s3      0x7F,0x41
 274      nop.i         0
 275 }
 276 {    .mfi
 277      mov           GR_neg_ov_limit = 0x33fff // Exponent for negative overflow
 278      fsetc.s2      0x7F,0x42
 279      nop.i         0
 280 }
 281 ;;
 282
 283 //
 284 //   Do final operation with s2 and s3
 285 //
 286 {    .mfi
 287      setf.exp      FR_NBig = GR_neg_ov_limit
 288      fma.s3        FR_Result3 = FR_Two_N,FR_Norm_X,f0
 289      nop.i         0
 290 }
 291 {    .mfi
 292      setf.exp      FR_Big = GR_pos_ov_limit
 293      fma.s2        FR_Result2 = FR_Two_N,FR_Norm_X,f0
 294      nop.i         0
 295 }
 296 ;;
 297
 298 //   Check for overflow or underflow.
 299 //   Restore s3
 300 //   Restore s2
 301 //
 302 {    .mfi
 303      nop.m         0
 304      fsetc.s3      0x7F,0x40
 305      nop.i         0
 306 }
 307 {    .mfi
 308      nop.m         0
 309      fsetc.s2      0x7F,0x40
 310      nop.i         0
 311 }
 312 ;;
 313
 314 //
 315 //   Is the result zero?
 316 //
 317 {    .mfi
 318      nop.m         0
 319      fclass.m      p6, p0 =  FR_Result3, 0x007
 320      nop.i         0
 321 }
 322 {    .mfi
 323      nop.m         0
 324      fcmp.ge.s1    p7, p8 = FR_Result2 , FR_Big
 325      nop.i         0
 326 }
 327 ;;
 328
 329 //
 330 //   Detect masked underflow - Tiny + Inexact Only
 331 //
 332 {    .mfi
 333      nop.m         0
 334 (p6) fcmp.neq.unc.s1 p6, p0 = FR_Result , FR_Result2
 335      nop.i         0
 336 }
 337 ;;
 338
 339 //
 340 //   Is result bigger the allowed range?
 341 //   Branch out for underflow
 342 //
 343 {    .mfb
 344      nop.m          0
 345 (p8) fcmp.le.unc.s1 p9, p10 = FR_Result2 , FR_NBig
 346 (p6) br.cond.spnt   SCALBL_UNDERFLOW
 347 }
 348 ;;
 349
 350 //
 351 //   Branch out for overflow
 352 //
 353 { .bbb
 354 (p7) br.cond.spnt   SCALBL_OVERFLOW
 355 (p9) br.cond.spnt   SCALBL_OVERFLOW
 356      br.ret.sptk    b0             //   Return from main path.
 357 }
 358 ;;
 359
 360 // Here if result overflows
 361 SCALBL_OVERFLOW:
 362 { .mib
 363      alloc         r32=ar.pfs,3,0,4,0
 364      addl          GR_Tag = 51, r0     // Set error tag for overflow
 365      br.cond.sptk  __libm_error_region // Call error support for overflow
 366 }
 367 ;;
 368
 369 // Here if result underflows
 370 SCALBL_UNDERFLOW:
 371 { .mib
 372      alloc         r32=ar.pfs,3,0,4,0
 373      addl          GR_Tag = 52, r0     // Set error tag for underflow
 374      br.cond.sptk  __libm_error_region // Call error support for underflow
 375 }
 376 ;;
 377
 378 SCALBL_NAN_INF_ZERO:
 379
 380 //
 381 //   Before entry, N has been converted to a fp integer in significand of
 382 //     FR_N_float_int
 383 //
 384 //   Convert  N_float_int to floating point value
 385 //
 386 {    .mfi
 387      getf.sig     GR_N_as_int = FR_N_float_int
 388      fclass.m     p6,p0 = FR_Floating_N, 0xc3 //@snan | @qnan
 389      nop.i        0
 390 }
 391 {    .mfi
 392      addl         GR_Scratch = 1,r0
 393      fcvt.xf      FR_N_float_int = FR_N_float_int
 394      nop.i        0
 395 }
 396 ;;
 397
 398 {    .mfi
 399      nop.m        0
 400      fclass.m     p7,p0 = FR_Floating_X, 0xc3 //@snan | @qnan
 401      shl          GR_Scratch = GR_Scratch,63
 402 }
 403 ;;
 404
 405 {    .mfi
 406      nop.m        0
 407      fclass.m     p8,p0 = FR_Floating_N, 0x21 // @inf
 408      nop.i        0
 409 }
 410 {    .mfi
 411      nop.m        0
 412      fclass.m     p9,p0 = FR_Floating_N, 0x22 // @-inf
 413      nop.i        0
 414 }
 415 ;;
 416
 417 //
 418 //   Either X or N is a Nan, return result and possible raise invalid.
 419 //
 420 {    .mfb
 421      nop.m        0
 422 (p6) fma.s0       FR_Result = FR_Floating_N,FR_Floating_X,f0
 423 (p6) br.ret.spnt  b0
 424 }
 425 ;;
 426
 427 {    .mfb
 428      nop.m        0
 429 (p7) fma.s0       FR_Result = FR_Floating_N,FR_Floating_X,f0
 430 (p7) br.ret.spnt  b0
 431 }
 432 ;;
 433
 434 //
 435 //   If N + Inf do something special
 436 //   For N = -Inf, create Int
 437 //
 438 {    .mfb
 439      nop.m        0
 440 (p8) fma.s0       FR_Result = FR_Floating_X, FR_Floating_N,f0
 441 (p8) br.ret.spnt  b0
 442 }
 443 {    .mfi
 444      nop.m        0
 445 (p9) fnma.s0      FR_Floating_N = FR_Floating_N, f1, f0
 446      nop.i        0
 447 }
 448 ;;
 449
 450 //
 451 //   If N==-Inf,return x/(-N)
 452 //
 453 {    .mfb
 454      cmp.ne       p7,p0 = GR_N_as_int,GR_Scratch
 455 (p9) frcpa.s0     FR_Result,p0 = FR_Floating_X,FR_Floating_N
 456 (p9) br.ret.spnt  b0
 457 }
 458 ;;
 459
 460 //
 461 //   Is N an integer.
 462 //
 463 {    .mfi
 464      nop.m        0
 465 (p7) fcmp.neq.unc.s1 p7,p0 = FR_Norm_N, FR_N_float_int
 466      nop.i        0
 467 }
 468 ;;
 469
 470 //
 471 //   If N not an int, return NaN and raise invalid.
 472 //
 473 {    .mfb
 474      nop.m        0
 475 (p7) frcpa.s0     FR_Result,p0 = f0,f0
 476 (p7) br.ret.spnt  b0
 477 }
 478 ;;
 479
 480 //
 481 //   Always return x in other path.
 482 //
 483 {    .mfb
 484      nop.m        0
 485      fma.s0       FR_Result = FR_Floating_X,f1,f0
 486      br.ret.sptk  b0
 487 }
 488 ;;
 489
 490 // Here if n not int
 491 // Return NaN and raise invalid.
 492 SCALBL_N_NOT_INT:
 493 {    .mfb
 494      nop.m        0
 495      frcpa.s0     FR_Result,p0 = f0,f0
 496      br.ret.sptk  b0
 497 }
 498 ;;
 499
 500 // Here if n=unorm
 501 SCALBL_N_UNORM:
 502 { .mfb
 503      getf.exp      GR_signexp_N = FR_Norm_N // Get signexp of normalized n
 504      fcvt.fx.trunc.s1   FR_N_float_int = FR_Norm_N // Get N in significand
 505      br.cond.sptk  SCALBL_COMMON1            // Return to main path
 506 }
 507 ;;
 508
 509 // Here if x=unorm
 510 SCALBL_X_UNORM:
 511 { .mib
 512      getf.exp      GR_signexp_X = FR_Norm_X // Get signexp of normalized x
 513      nop.i         0
 514      br.cond.sptk  SCALBL_COMMON2            // Return to main path
 515 }
 516 ;;
 517
 518 GLOBAL_IEEE754_END(scalbl)
 519 LOCAL_LIBM_ENTRY(__libm_error_region)
 520
 521 //
 522 // Get stack address of N
 523 //
 524 .prologue
 525 { .mfi
 526     add   GR_Parameter_Y=-32,sp
 527     nop.f 0
 528 .save   ar.pfs,GR_SAVE_PFS
 529     mov  GR_SAVE_PFS=ar.pfs
 530 }
 531 //
 532 // Adjust sp
 533 //
 534 { .mfi
 535 .fframe 64
 536    add sp=-64,sp
 537    nop.f 0
 538    mov GR_SAVE_GP=gp
 539 };;
 540
 541 //
 542 //  Store N on stack in correct position
 543 //  Locate the address of x on stack
 544 //
 545 { .mmi
 546    stfe [GR_Parameter_Y] = FR_Norm_N,16
 547    add GR_Parameter_X = 16,sp
 548 .save   b0, GR_SAVE_B0
 549    mov GR_SAVE_B0=b0
 550 };;
 551
 552 //
 553 // Store x on the stack.
 554 // Get address for result on stack.
 555 //
 556 .body
 557 { .mib
 558    stfe [GR_Parameter_X] = FR_Norm_X
 559    add   GR_Parameter_RESULT = 0,GR_Parameter_Y
 560    nop.b 0
 561 }
 562 { .mib
 563    stfe [GR_Parameter_Y] = FR_Result
 564    add   GR_Parameter_Y = -16,GR_Parameter_Y
 565    br.call.sptk b0=__libm_error_support#
 566 };;
 567
 568 //
 569 //  Get location of result on stack
 570 //
 571 { .mmi
 572    add   GR_Parameter_RESULT = 48,sp
 573    nop.m 0
 574    nop.i 0
 575 };;
 576
 577 //
 578 //  Get the new result
 579 //
 580 { .mmi
 581    ldfe  FR_Result = [GR_Parameter_RESULT]
 582 .restore sp
 583    add   sp = 64,sp
 584    mov   b0 = GR_SAVE_B0
 585 };;
 586
 587 //
 588 //  Restore gp, ar.pfs and return
 589 //
 590 { .mib
 591    mov   gp = GR_SAVE_GP
 592    mov   ar.pfs = GR_SAVE_PFS
 593    br.ret.sptk     b0
 594 };;
 595
 596 LOCAL_LIBM_END(__libm_error_region)
 597
 598 .type   __libm_error_support#,@function
 599 .global __libm_error_support#