sysdeps/ia64/fpu/e_scalb.S

   1 .file "scalb.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 //
  47 // API
  48 //==============================================================
  49 // double = scalb  (double x, double n)
  50 // input  floating point f8 and floating point f9
  51 // output floating point f8
  52 //
  53 // Returns x* 2**n using an fma and detects overflow
  54 // and underflow.
  55 //
  56 //
  57
  58 FR_Floating_X  = f8
  59 FR_Result      = f8
  60 FR_Floating_N  = f9
  61 FR_Result2     = f9
  62 FR_Norm_N      = f10
  63 FR_Result3     = f11
  64 FR_Norm_X      = f12
  65 FR_N_float_int = f13
  66 FR_Two_N       = f14
  67 FR_Two_to_Big  = f15
  68 FR_Big         = f6
  69 FR_NBig        = f7
  70
  71 GR_N_Biased    = r15
  72 GR_Big         = r16
  73 GR_NBig        = r17
  74 GR_Scratch     = r18
  75 GR_Scratch1    = r19
  76 GR_Bias        = r20
  77 GR_N_as_int    = r21
  78
  79 GR_SAVE_B0          = r32
  80 GR_SAVE_GP          = r33
  81 GR_SAVE_PFS         = r34
  82 GR_Parameter_X      = r35
  83 GR_Parameter_Y      = r36
  84 GR_Parameter_RESULT = r37
  85 GR_Tag              = r38
  86
  87 .section .text
  88 GLOBAL_IEEE754_ENTRY(scalb)
  89
  90 //
  91 //   Is x NAN, INF, ZERO, +-?
  92 //
  93 {    .mfi
  94      alloc          r32=ar.pfs,0,3,4,0
  95      fclass.m.unc  p7,p0 = FR_Floating_X, 0xe7 //@snan | @qnan | @inf | @zero
  96      addl  GR_Scratch  = 0x019C3F,r0
  97 }
  98 //
  99 //   Is y a NAN, INF, ZERO, +-?
 100 //
 101 {    .mfi
 102      nop.m 999
 103      fclass.m.unc  p6,p0 = FR_Floating_N, 0xe7 //@snan | @qnan | @inf |  @zero
 104      addl  GR_Scratch1  = 0x063BF,r0
 105 }
 106 ;;
 107
 108 //
 109 //   Convert N to a fp integer
 110 //   Normalize x
 111 //
 112 {    .mfi
 113      nop.m 0
 114      fnorm.s1  FR_Norm_N  =   FR_Floating_N
 115      nop.i 999
 116 }
 117 {    .mfi
 118      nop.m 999
 119      fnorm.s1  FR_Norm_X  =   FR_Floating_X
 120      nop.i 999
 121 };;
 122
 123 //
 124 //   Create 2*big
 125 //   Create 2**-big
 126 //   Normalize x
 127 //   Branch on special values.
 128 //
 129 { .mib
 130      setf.exp      FR_Big = GR_Scratch
 131      nop.i 0
 132 (p6) br.cond.spnt  SCALB_NAN_INF_ZERO
 133 }
 134 { .mib
 135      setf.exp      FR_NBig = GR_Scratch1
 136      nop.i 0
 137 (p7) br.cond.spnt  SCALB_NAN_INF_ZERO
 138 };;
 139
 140 //
 141 //   Convert N to a fp integer
 142 //   Create -35000
 143 //
 144 {    .mfi
 145      addl  GR_Scratch = 1,r0
 146      fcvt.fx.trunc.s1   FR_N_float_int = FR_Norm_N
 147      addl    GR_NBig = -35000,r0
 148 }
 149 ;;
 150
 151 //
 152 //   Put N if a GP register
 153 //   Convert  N_float_int to floating point value
 154 //   Create 35000
 155 //   Build the exponent Bias
 156 //
 157 {    .mii
 158      getf.sig     GR_N_as_int = FR_N_float_int
 159      shl   GR_Scratch = GR_Scratch,63
 160      addl  GR_Big = 35000,r0
 161 }
 162 {    .mfi
 163      addl GR_Bias = 0x0FFFF,r0
 164      fcvt.xf  FR_N_float_int = FR_N_float_int
 165      nop.i 0
 166 };;
 167
 168 //
 169 //   Catch those fp values that are beyond 2**64-1
 170 //   Is N > 35000
 171 //   Is N < -35000
 172 //
 173 {     .mfi
 174      cmp.ne.unc  p9,p10 = GR_N_as_int,GR_Scratch
 175      nop.f 0
 176      nop.i 0
 177 }
 178 {     .mmi
 179      cmp.ge.unc p6, p0 = GR_N_as_int, GR_Big
 180      cmp.le.unc p8, p0 = GR_N_as_int, GR_NBig
 181      nop.i 0
 182 };;
 183
 184 //
 185 //   Is N really an int, only for those non-int indefinites?
 186 //   Create exp bias.
 187 //
 188 {    .mfi
 189      add GR_N_Biased = GR_Bias,GR_N_as_int
 190 (p9) fcmp.neq.unc.s1 p7,p0  =   FR_Norm_N, FR_N_float_int
 191      nop.i 0
 192 };;
 193
 194 //
 195 //   Branch and return if N is not an int.
 196 //   Main path, create 2**N
 197 //
 198 {    .mfi
 199      setf.exp      FR_Two_N = GR_N_Biased
 200      nop.i                      999
 201 }
 202 {    .mfb
 203      nop.m 0
 204 (p7) frcpa.s0          f8,p11     =    f0,f0
 205 (p7) br.ret.spnt    b0
 206 };;
 207
 208 //
 209 //   Set denormal on denormal input x and denormal input N
 210 //
 211 {    .mfi
 212      nop.m                      999
 213 (p10)fcmp.ge.s1    p6,p8 = FR_Norm_N,f0
 214      nop.i 0
 215 };;
 216 {    .mfi
 217      nop.m                      999
 218      fcmp.ge.s0    p0,p11 = FR_Floating_X,f0
 219      nop.i                      999
 220 }
 221 {    .mfi
 222      nop.m                      999
 223      fcmp.ge.s0    p12,p13 = FR_Floating_N,f0
 224      nop.i 0
 225 };;
 226
 227 //
 228 //   Adjust 2**N if N was very small or very large
 229 //
 230
 231 {    .mfi
 232      nop.m 0
 233 (p6) fma.s1  FR_Two_N = FR_Big,f1,f0
 234      nop.i 0
 235 }
 236 { .mlx
 237      nop.m 999
 238      movl GR_Scratch = 0x00000000000303FF
 239 };;
 240 {    .mfi
 241      nop.m 0
 242 (p8) fma.s1  FR_Two_N = FR_NBig,f1,f0
 243      nop.i 0
 244 }
 245 {    .mlx
 246      nop.m 999
 247      movl GR_Scratch1= 0x00000000000103FF
 248 };;
 249
 250 //   Set up necessary status fields
 251 //
 252 //   S0 user supplied status
 253 //   S2 user supplied status + WRE + TD  (Overflows)
 254 //   S3 user supplied status + FZ + TD   (Underflows)
 255 //
 256 {    .mfi
 257      nop.m 999
 258      fsetc.s3 0x7F,0x41
 259      nop.i 999
 260 }
 261 {    .mfi
 262      nop.m 999
 263      fsetc.s2 0x7F,0x42
 264      nop.i 999
 265 };;
 266
 267 //
 268 //   Do final operation
 269 //
 270 {    .mfi
 271      setf.exp FR_NBig = GR_Scratch
 272      fma.d.s0     FR_Result = FR_Two_N,FR_Norm_X,f0
 273      nop.i                           999
 274 }
 275 {    .mfi
 276      nop.m                           999
 277      fma.d.s3     FR_Result3 = FR_Two_N,FR_Norm_X,f0
 278      nop.i                           999
 279 };;
 280 {    .mfi
 281      setf.exp FR_Big = GR_Scratch1
 282      fma.d.s2     FR_Result2 = FR_Two_N,FR_Norm_X,f0
 283      nop.i                           999
 284 };;
 285
 286 //   Check for overflow or underflow.
 287 //
 288 //   S0 user supplied status
 289 //   S2 user supplied status + WRE + TD  (Overflow)
 290 //   S3 user supplied status + FZ + TD   (Underflow)
 291 //
 292 //
 293 //   Restore s3
 294 //   Restore s2
 295 //
 296 {    .mfi
 297      nop.m 0
 298      fsetc.s3 0x7F,0x40
 299      nop.i 999
 300 }
 301 {    .mfi
 302      nop.m 0
 303      fsetc.s2 0x7F,0x40
 304      nop.i 999
 305 };;
 306
 307 //
 308 //   Is the result zero?
 309 //
 310 {    .mfi
 311      nop.m 999
 312      fclass.m.unc   p6, p0 =  FR_Result3, 0x007
 313      nop.i 999
 314 }
 315 {    .mfi
 316      addl GR_Tag = 53, r0
 317      fcmp.ge.unc.s1 p7, p8 = FR_Result2 , FR_Big
 318      nop.i 0
 319 };;
 320
 321 //
 322 //   Detect masked underflow - Tiny + Inexact Only
 323 //
 324 {    .mfi
 325      nop.m 999
 326 (p6) fcmp.neq.unc.s1 p6, p0 = FR_Result , FR_Result2
 327      nop.i 999
 328 };;
 329
 330 //
 331 //   Is result bigger the allowed range?
 332 //   Branch out for underflow
 333 //
 334 {    .mfb
 335 (p6) addl GR_Tag = 54, r0
 336 (p8) fcmp.le.unc.s1 p9, p10 = FR_Result2 , FR_NBig
 337 (p6) br.cond.spnt SCALB_UNDERFLOW
 338 };;
 339
 340 //
 341 //   Branch out for overflow
 342 //
 343 { .mbb
 344      nop.m 0
 345 (p7) br.cond.spnt SCALB_OVERFLOW
 346 (p9) br.cond.spnt SCALB_OVERFLOW
 347 };;
 348
 349 //
 350 //   Return from main path.
 351 //
 352 {    .mfb
 353      nop.m 999
 354      nop.f 0
 355      br.ret.sptk     b0;;
 356 }
 357
 358 SCALB_NAN_INF_ZERO:
 359
 360 //
 361 //   Convert N to a fp integer
 362 //
 363 {    .mfi
 364      addl  GR_Scratch = 1,r0
 365      fcvt.fx.trunc.s1  FR_N_float_int = FR_Norm_N
 366      nop.i 999
 367 }
 368 {    .mfi
 369      nop.m 0
 370      fclass.m.unc  p6,p0 = FR_Floating_N, 0xc3 //@snan | @qnan
 371      nop.i 0
 372 };;
 373 {    .mfi
 374      nop.m 0
 375      fclass.m.unc  p7,p0 = FR_Floating_X, 0xc3 //@snan | @qnan
 376      shl   GR_Scratch = GR_Scratch,63
 377 };;
 378 {    .mfi
 379      nop.m 0
 380      fclass.m.unc  p8,p0 = FR_Floating_N, 0x21 // @inf
 381      nop.i 0
 382 }
 383   {  .mfi
 384      nop.m 0
 385      fclass.m.unc  p9,p0 = FR_Floating_N, 0x22 // @-inf
 386      nop.i 0
 387 };;
 388
 389 //
 390 //   Either X or N is a Nan, return result and possible raise invalid.
 391 //
 392 {    .mfb
 393      nop.m 0
 394 (p6) fma.d.s0     FR_Result = FR_Floating_N,FR_Floating_X,f0
 395 (p6) br.ret.spnt  b0
 396 };;
 397 {    .mfb
 398      getf.sig     GR_N_as_int = FR_N_float_int
 399 (p7) fma.d.s0     FR_Result = FR_Floating_N,FR_Floating_X,f0
 400 (p7) br.ret.spnt  b0
 401 };;
 402
 403 //
 404 //   If N + Inf do something special
 405 //   For N = -Inf, create Int
 406 //
 407 {    .mfb
 408      nop.m 0
 409 (p8) fma.d.s0    FR_Result = FR_Floating_X, FR_Floating_N,f0
 410 (p8) br.ret.spnt   b0
 411 }
 412 {    .mfi
 413      nop.m 0
 414 (p9) fnma.d.s0   FR_Floating_N = FR_Floating_N, f1, f0
 415      nop.i 0
 416 };;
 417
 418 //
 419 //   If N==-Inf,return x/(-N)
 420 //
 421 {    .mfb
 422      nop.m 0
 423 (p9) frcpa.s0        FR_Result,p6 =  FR_Floating_X,FR_Floating_N
 424 (p9) br.ret.spnt    b0
 425 };;
 426
 427 //
 428 //   Convert  N_float_int to floating point value
 429 //
 430 {     .mfi
 431      cmp.ne.unc  p9,p0     =   GR_N_as_int,GR_Scratch
 432      fcvt.xf  FR_N_float_int = FR_N_float_int
 433      nop.i  0
 434 };;
 435
 436 //
 437 //   Is N an integer.
 438 //
 439 {    .mfi
 440      nop.m 0
 441 (p9) fcmp.neq.unc.s1 p7,p0  =   FR_Norm_N, FR_N_float_int
 442      nop.i 0
 443 };;
 444
 445 //
 446 //   If N not an int, return NaN and raise invalid.
 447 //
 448 {    .mfb
 449      nop.m 0
 450 (p7) frcpa.s0        FR_Result,p6     =    f0,f0
 451 (p7) br.ret.spnt    b0
 452 };;
 453
 454 //
 455 //   Always return x in other path.
 456 //
 457 {    .mfb
 458      nop.m 0
 459      fma.d.s0      FR_Result = FR_Floating_X,f1,f0
 460      br.ret.sptk   b0
 461 };;
 462
 463 GLOBAL_IEEE754_END(scalb)
 464 __libm_error_region:
 465
 466 SCALB_OVERFLOW:
 467 SCALB_UNDERFLOW:
 468
 469 //
 470 // Get stack address of N
 471 //
 472 .prologue
 473 { .mfi
 474     add   GR_Parameter_Y=-32,sp
 475     nop.f 0
 476 .save   ar.pfs,GR_SAVE_PFS
 477     mov  GR_SAVE_PFS=ar.pfs
 478 }
 479 //
 480 // Adjust sp
 481 //
 482 { .mfi
 483 .fframe 64
 484    add sp=-64,sp
 485    nop.f 0
 486    mov GR_SAVE_GP=gp
 487 };;
 488
 489 //
 490 //  Store N on stack in correct position
 491 //  Locate the address of x on stack
 492 //
 493 { .mmi
 494    stfd [GR_Parameter_Y] = FR_Norm_N,16
 495    add GR_Parameter_X = 16,sp
 496 .save   b0, GR_SAVE_B0
 497    mov GR_SAVE_B0=b0
 498 };;
 499
 500 //
 501 // Store x on the stack.
 502 // Get address for result on stack.
 503 //
 504 .body
 505 { .mib
 506    stfd [GR_Parameter_X] = FR_Norm_X
 507    add   GR_Parameter_RESULT = 0,GR_Parameter_Y
 508    nop.b 0
 509 }
 510 { .mib
 511    stfd [GR_Parameter_Y] = FR_Result
 512    add   GR_Parameter_Y = -16,GR_Parameter_Y
 513    br.call.sptk b0=__libm_error_support#
 514 };;
 515
 516 //
 517 //  Get location of result on stack
 518 //
 519 { .mmi
 520    nop.m 0
 521    nop.m 0
 522    add   GR_Parameter_RESULT = 48,sp
 523 };;
 524
 525 //
 526 //  Get the new result
 527 //
 528 { .mmi
 529    ldfd  FR_Result = [GR_Parameter_RESULT]
 530 .restore sp
 531    add   sp = 64,sp
 532    mov   b0 = GR_SAVE_B0
 533 };;
 534
 535 //
 536 //  Restore gp, ar.pfs and return
 537 //
 538 { .mib
 539    mov   gp = GR_SAVE_GP
 540    mov   ar.pfs = GR_SAVE_PFS
 541    br.ret.sptk     b0
 542 };;
 543
 544 LOCAL_LIBM_END(__libm_error_region)
 545
 546 .type   __libm_error_support#,@function
 547 .global __libm_error_support#