sysdeps/ia64/fpu/e_scalb.S

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