sysdeps/ia64/fpu/e_hypotl.S

   1 .file "hypotl.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 //*********************************************************************
  41 //
  42 // History:
  43 // 02/02/00 hand-optimized
  44 // 04/04/00 Unwind support added
  45 // 06/20/00 new version
  46 // 08/15/00 Bundle added after call to __libm_error_support to properly
  47 //          set [the previously overwritten] GR_Parameter_RESULT.
  48 // 05/20/02 Cleaned up namespace and sf0 syntax
  49 // 02/10/03 Reordered header: .section, .global, .proc, .align
  50 //
  51 //*********************************************************************
  52 //                           ___________
  53 // Function:   hypotl(x,y) = |(x^2 + y^2) = for double extended values
  54 //             x and y
  55 //             Also provides cabsl functionality.
  56 //
  57 //*********************************************************************
  58 //
  59 // Resources Used:
  60 //
  61 //    Floating-Point Registers: f8  (Input and Return Value)
  62 //                              f9  (Input)
  63 //                              f6 -f15, f32-f34
  64 //
  65 //    General Purpose Registers:
  66 //      r2-r3   (Scratch)
  67 //      r32-r36 (Locals)
  68 //      r37-r40 (Used to pass arguments to error handling routine)
  69 //
  70 //    Predicate Registers:      p6 - p10
  71 //
  72 //*********************************************************************
  73 //
  74 // IEEE Special Conditions:
  75 //
  76 //    All faults and exceptions should be raised correctly.
  77 //    Overflow can occur.
  78 //    hypotl(Infinity and anything) = +Infinity
  79 //    hypotl(QNaN and anything) = QNaN
  80 //    hypotl(SNaN and anything ) = QNaN
  81 //
  82 //*********************************************************************
  83 //
  84 // Implementation:
  85 //  x2 = x * x   in double-extended
  86 //  y2 = y * y   in double-extended
  87 //  temp = x2 + y2   in double-extended
  88 //  sqrt(temp) rounded to double extended
  89 //
  90 //*********************************************************************
  91
  92 GR_SAVE_PFS         = r33
  93 GR_SAVE_B0          = r34
  94 GR_SAVE_GP          = r35
  95 GR_Parameter_X      = r36
  96 GR_Parameter_Y      = r37
  97 GR_Parameter_RESULT = r38
  98 GR_Parameter_TAG    = r39
  99
 100 FR_X                = f32
 101 FR_Y                = f33
 102 FR_RESULT           = f8
 103
 104 .section .text
 105
 106 LOCAL_LIBM_ENTRY(cabsl)
 107 LOCAL_LIBM_END(cabsl)
 108
 109 GLOBAL_IEEE754_ENTRY(hypotl)
 110 {.mfi
 111   alloc r32= ar.pfs,0,4,4,0
 112   // Compute x*x
 113   fma.s1 f10=f8,f8,f0
 114   // r2=bias-1
 115   mov r2=0xfffe
 116 }
 117 {.mfi
 118   nop.m 0
 119   // y*y
 120   fma.s1 f11=f9,f9,f0
 121   nop.i 0;;
 122 }
 123
 124 { .mfi
 125      nop.m 0
 126 //   Check if x is an Inf - if so return Inf even
 127 //   if y is a NaN (C9X)
 128      fclass.m.unc p7, p6 = f8, 0x023
 129      nop.i 0
 130 }
 131 {.mfi
 132         nop.m 0
 133         // if possible overflow, copy f8 to f32
 134         // set Denormal, if necessary
 135         // (p8)
 136         fma.s0 f32=f8,f1,f0
 137     nop.i 0;;
 138 }
 139 { .mfi
 140      nop.m 0
 141 //   Check if y is an Inf - if so return Inf even
 142 //   if x is a NaN (C9X)
 143      fclass.m.unc p8, p9 = f9, 0x023
 144          nop.i 0
 145 }
 146 { .mfi
 147      nop.m 999
 148 //   For x=inf, multiply y by 1 to raise invalid on y an SNaN
 149 //   (p7) fma.s0 f9=f9,f1,f0
 150      // copy f9 to f33; set Denormal, if necessary
 151          fma.s0 f33=f9,f1,f0
 152      nop.i 0;;
 153 }
 154 {.mfi
 155   nop.m 0
 156   // is y Zero ?
 157   (p6) fclass.m p6,p0=f9,0x7
 158   nop.i 0;;
 159 }
 160
 161 {.mfi
 162   // f7=0.5
 163   setf.exp f7=r2
 164   // a=x2+y2
 165   fma.s1 f12=f10,f1,f11
 166   nop.i 0
 167 }
 168 {.mfi
 169   mov r2=0x408c  //0000
 170   // dx=x*x-x2
 171   fms.s1 f13=f8,f8,f10
 172   nop.i 0;;
 173 }
 174 {.mfi
 175   nop.m 0
 176   // is x Zero ?
 177   (p9) fclass.m p9,p0=f8,0x7
 178   shl r2=r2,16
 179 }
 180 {.mfi
 181   nop.m 0
 182   // dy=y*y-y2
 183   fms.s1 f14=f9,f9,f11
 184   nop.i 0;;
 185 }
 186
 187 {.mfi
 188   nop.m 0
 189   // x not NaN ?
 190   (p6) fclass.m p7,p0=f8,0x3f
 191   nop.i 0
 192 }
 193 {.mfi
 194   nop.m 0
 195   // f6=2
 196   fma.s1 f6=f1,f1,f1
 197   nop.i 0;;
 198 }
 199
 200 {.mfi
 201   nop.m 0
 202   // f34=min(x2,y2)
 203   famin.s1 f34=f10,f11
 204   nop.i 0
 205 }
 206 {.mfb
 207   nop.m 0
 208   // f10=max(x2,y2)
 209   famax.s1 f10=f11,f10
 210   nop.b 0;; //
 211 }
 212
 213 {.mfi
 214   nop.m 0
 215   // y not NaN ?
 216   (p9) fclass.m p8,p0=f9,0x3f
 217   nop.i 0;;
 218 }
 219 {.mfb
 220   // f9=35/8
 221   setf.s f9=r2
 222   // if f8=Infinity or f9=Zero, return |f8|
 223   (p7) fmerge.s f8=f0,f32
 224   (p7) br.ret.spnt b0;;
 225 }
 226
 227
 228 {.mfi
 229   nop.m 0
 230   // z0=frsqrta(a)
 231   frsqrta.s1 f8,p6=f12
 232   nop.i 0;;
 233 }
 234 { .mfi
 235      nop.m 0
 236 //   Identify Natvals, Infs, NaNs, and Zeros
 237 //   and return result
 238      fclass.m.unc p7, p0 = f12, 0x1E7
 239      nop.i 0
 240 }
 241 {.mfi
 242   // get exponent of x^2+y^2
 243   getf.exp r3=f12
 244   // dxy=dx+dy
 245   fma.s1 f13=f13,f1,f14
 246   nop.i 0;;
 247 }
 248
 249 {.mfb
 250   // 2*emax-2
 251   mov r2=0x17ffb
 252   // if f9=Infinity or f8=Zero, return |f9|
 253   (p8) fmerge.s f8=f0,f33
 254   (p8) br.ret.spnt b0
 255 }
 256 {.mfi
 257   nop.m 0
 258   // dd=a-max(x2,y2)
 259   fnma.s1 f10=f10,f1,f12
 260   nop.i 0;;
 261 }
 262
 263 {.mfi
 264   nop.m 0
 265   // S0=a*z0
 266   (p6) fma.s1 f14=f12,f8,f0
 267   nop.i 0
 268 }
 269 {.mfi
 270   nop.m 0
 271   // H0=0.5*z0
 272   (p6) fma.s1 f15=f8,f7,f0
 273   nop.i 0;;
 274 }
 275
 276 {.mfb
 277   nop.m 0
 278   // if special case, set f8
 279   (p7) mov f8=f12
 280   (p7) br.ret.spnt b0
 281 }
 282 {.mfi
 283   nop.m 0
 284   // da=min(x2,y2)-dd
 285   fnma.s1 f10=f10,f1,f34
 286   nop.i 0;;
 287 }
 288 {.mfi
 289   nop.m 0
 290   // f6=5/2
 291   fma.s1 f6=f7,f1,f6
 292   nop.i 0
 293 }
 294 {.mfi
 295   nop.m 0
 296   // f11=3/2
 297   fma.s1 f11=f7,f1,f1
 298   nop.i 0;;
 299 }
 300
 301 {.mfi
 302   nop.m 0
 303   // d=0.5-S0*H0
 304   (p6) fnma.s1 f7=f14,f15,f7
 305   nop.i 0;;
 306 }
 307
 308 {.mfi
 309   nop.m 0
 310   // P1=3/2*d+1
 311   (p6) fma.s1 f11=f11,f7,f1
 312   nop.i 0
 313 }
 314 {.mfi
 315   nop.m 0
 316   // P2=35/8*d+5/2
 317   (p6) fma.s1 f9=f9,f7,f6
 318   nop.i 0;;
 319 }
 320 {.mfi
 321   nop.m 0
 322   // d2=d*d
 323   (p6) fma.s1 f34=f7,f7,f0
 324   nop.i 0;;
 325 }
 326
 327 {.mfi
 328   nop.m 0
 329   // T0=d*S0
 330   (p6) fma.s1 f6=f7,f14,f0
 331   nop.i 0
 332 }
 333 {.mfi
 334   nop.m 0
 335   // G0=d*H0
 336   (p6) fma.s1 f7=f7,f15,f0
 337   nop.i 0;;
 338 }
 339 {.mfi
 340   nop.m 0
 341   // P=d2*P2+P1
 342   (p6) fma.s1 f11=f34,f9,f11
 343   nop.i 0;;
 344 }
 345
 346 {.mfi
 347   nop.m 0
 348   // S1=p*T0+S0
 349   (p6) fma.s1 f14=f11,f6,f14
 350   nop.i 0
 351 }
 352 {.mfi
 353   nop.m 0
 354   // H1=p*G0+H0
 355   (p6) fma.s1 f15=f11,f7,f15
 356   nop.i 0;;
 357 }
 358
 359
 360 {.mfi
 361   nop.m 0
 362   // e1=a-S1*S1
 363   (p6) fnma.s1 f7=f14,f14,f12
 364   nop.i 0
 365 }
 366 {.mfi
 367   // Is x^2 + y^2 well less than the overflow
 368   // threshold?
 369   (p6) cmp.lt.unc p7, p8 =  r3,r2
 370   // c=dxy+da
 371   (p6) fma.s1 f13=f13,f1,f10
 372   nop.i 0;;
 373 }
 374
 375 {.mfi
 376   nop.m 0
 377   // e=e1+c
 378   (p6) fma.s1 f13=f7,f1,f13
 379   nop.i 0;;
 380 }
 381
 382 {.mfb
 383   nop.m 0
 384   // S=e*H1+S1
 385   fma.s0 f8=f13,f15,f14
 386   // No overflow in this case
 387   (p7) br.ret.sptk b0;;
 388 }
 389
 390 { .mfi
 391      nop.m 0
 392 (p8) fsetc.s2 0x7F,0x42
 393      // Possible overflow path, must detect by
 394      // Setting widest range exponent with prevailing
 395      // rounding mode.
 396      nop.i 0 ;;
 397 }
 398
 399
 400 { .mfi
 401    // bias+0x4000 (bias+EMAX+1)
 402    (p8) mov r2=0x13fff
 403    // S=e*H1+S1
 404    (p8) fma.s2 f12=f13,f15,f14
 405    nop.i 0 ;;
 406 }
 407 { .mfi
 408 (p8) setf.exp f11 = r2
 409 (p8) fsetc.s2 0x7F,0x40
 410 //   Restore Original Mode in S2
 411      nop.i 0 ;;
 412 }
 413 { .mfi
 414      nop.m 0
 415 (p8) fcmp.lt.unc.s1 p9, p10 =  f12, f11
 416      nop.i 0 ;;
 417 }
 418 { .mib
 419      nop.m 0
 420      mov GR_Parameter_TAG = 45;
 421          // No overflow
 422 (p9) br.ret.sptk b0;;
 423 }
 424 GLOBAL_IEEE754_END(hypotl)
 425
 426 LOCAL_LIBM_ENTRY(__libm_error_region)
 427 .prologue
 428 { .mfi
 429         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 430         nop.f 0
 431 .save   ar.pfs,GR_SAVE_PFS
 432         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 433 }
 434 { .mfi
 435 .fframe 64
 436         add sp=-64,sp                           // Create new stack
 437         nop.f 0
 438         mov GR_SAVE_GP=gp                       // Save gp
 439 };;
 440 { .mmi
 441         stfe [GR_Parameter_Y] = FR_Y,16         // Save Parameter 2 on stack
 442         add GR_Parameter_X = 16,sp              // Parameter 1 address
 443 .save   b0, GR_SAVE_B0
 444         mov GR_SAVE_B0=b0                       // Save b0
 445 };;
 446 .body
 447 { .mib
 448         stfe [GR_Parameter_X] = FR_X            // Store Parameter 1 on stack
 449         add   GR_Parameter_RESULT = 0,GR_Parameter_Y
 450         nop.b 0                                 // Parameter 3 address
 451 }
 452 { .mib
 453         stfe [GR_Parameter_Y] = FR_RESULT      // Store Parameter 3 on stack
 454         add   GR_Parameter_Y = -16,GR_Parameter_Y
 455         br.call.sptk b0=__libm_error_support#  // Call error handling function
 456 };;
 457 { .mmi
 458         nop.m 0
 459         nop.m 0
 460         add   GR_Parameter_RESULT = 48,sp
 461 };;
 462 { .mmi
 463         ldfe  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 464 .restore sp
 465         add   sp = 64,sp                       // Restore stack pointer
 466         mov   b0 = GR_SAVE_B0                  // Restore return address
 467 };;
 468 { .mib
 469         mov   gp = GR_SAVE_GP                  // Restore gp
 470         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 471         br.ret.sptk     b0                     // Return
 472 };;
 473 LOCAL_LIBM_END(__libm_error_region#)
 474 .type   __libm_error_support#,@function
 475 .global __libm_error_support#
 476
 477
 478