sysdeps/ia64/fpu/e_hypotf.S

   1 .file "hypotf.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/26/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 // 04/17/03 Added missing mutex directive
  51 //
  52 //*********************************************************************
  53 //                           ___________
  54 // Function:   hypotf(x,y) = |(x^2 + y^2) = for single precision values
  55 //             x and y
  56 //             Also provides cabsf functionality.
  57 //
  58 //*********************************************************************
  59 //
  60 // Resources Used:
  61 //
  62 //    Floating-Point Registers: f8  (Input and Return Value)
  63 //                              f9  (Input)
  64 //                              f6 -f15
  65 //
  66 //    General Purpose Registers:
  67 //      r2-r3   (Scratch)
  68 //      r32-r36 (Locals)
  69 //      r37-r40 (Used to pass arguments to error handling routine)
  70 //
  71 //    Predicate Registers:      p6 - p10
  72 //
  73 //*********************************************************************
  74 //
  75 // IEEE Special Conditions:
  76 //
  77 //    All faults and exceptions should be raised correctly.
  78 //    Overflow can occur.
  79 //    hypotf(Infinity and anything) = +Infinity
  80 //    hypotf(QNaN and anything) = QNaN
  81 //    hypotf(SNaN and anything ) = QNaN
  82 //
  83 //*********************************************************************
  84 //
  85 // Implementation:
  86 //  x2 = x * x   in double-extended
  87 //  y2 = y * y   in double-extended
  88 //  temp = x2 + y2   in double-extended
  89 //  sqrt(temp) rounded to single precision
  90 //
  91 //*********************************************************************
  92
  93 GR_SAVE_PFS         = r33
  94 GR_SAVE_B0          = r34
  95 GR_SAVE_GP          = r35
  96 GR_Parameter_X      = r36
  97 GR_Parameter_Y      = r37
  98 GR_Parameter_RESULT = r38
  99 GR_Parameter_TAG    = r39
 100
 101 FR_X                = f14
 102 FR_Y                = f15
 103 FR_RESULT           = f8
 104
 105 .section .text
 106
 107 LOCAL_LIBM_ENTRY(cabsf)
 108 LOCAL_LIBM_END(cabsf)
 109
 110 GLOBAL_IEEE754_ENTRY(hypotf)
 111 {.mfi
 112   alloc r32= ar.pfs,0,4,4,0
 113   // Compute x*x
 114   fma.s1 f10=f8,f8,f0
 115   // r2=bias-1
 116   mov r2=0xfffe
 117 }
 118 {.mfi
 119   nop.m 0
 120   // y*y
 121   fma.s1 f11=f9,f9,f0
 122   nop.i 0;;
 123 }
 124
 125 { .mfi
 126      nop.m 0
 127 //   Check if x is an Inf - if so return Inf even
 128 //   if y is a NaN (C9X)
 129      fclass.m.unc p7, p6 = f8, 0x023
 130      nop.i 0
 131 }
 132 {.mfi
 133   nop.m 0
 134   // if possible overflow, copy f8 to f14
 135   // set Denormal, if necessary
 136   // (p8)
 137   fma.s.s0 f14=f8,f1,f0
 138   nop.i 0;;
 139 }
 140
 141 { .mfi
 142      nop.m 0
 143 //   Check if y is an Inf - if so return Inf even
 144 //   if x is a NaN (C9X)
 145      fclass.m.unc p8, p9 = f9, 0x023
 146          nop.i 0
 147 }
 148 { .mfi
 149      nop.m 0
 150 //   For x=inf, multiply y by 1 to raise invalid on y an SNaN
 151 //   (p7) fma.s0 f9=f9,f1,f0
 152      // copy f9 to f15; set Denormal, if necessary
 153          fma.s.s0 f15=f9,f1,f0
 154      nop.i 0;;
 155 }
 156 {.mfi
 157   nop.m 0
 158   // is y Zero ?
 159   (p6) fclass.m p6,p0=f9,0x7
 160   nop.i 0;;
 161 }
 162 {.mfi
 163   nop.m 0
 164   // is x Zero ?
 165   (p9) fclass.m p9,p0=f8,0x7
 166   nop.i 0;;
 167 }
 168
 169 {.mfi
 170   // f7=0.5
 171   setf.exp f7=r2
 172   // a=x2+y2
 173   fma.s1 f12=f10,f1,f11
 174   nop.i 0;;
 175 }
 176
 177 {.mfi
 178   nop.m 0
 179   // x not NaN ?
 180   (p6) fclass.m p7,p0=f8,0x3f
 181   nop.i 0
 182 }
 183 {.mfi
 184   // 2*emax-2
 185   mov r2=0x100fb
 186   // f6=2
 187   fma.s1 f6=f1,f1,f1
 188   nop.i 0;;
 189 }
 190
 191 {.mfi
 192   nop.m 0
 193   // y not NaN ?
 194   (p9) fclass.m p8,p0=f9,0x3f
 195   nop.i 0;;
 196 }
 197
 198 .pred.rel "mutex",p7,p8
 199 {.mfb
 200   nop.m 0
 201   // if f8=Infinity or f9=Zero, return |f8|
 202   (p7) fmerge.s f8=f0,f14
 203   (p7) br.ret.spnt b0
 204 }
 205 {.mfb
 206   nop.m 0
 207   // if f9=Infinity or f8=Zero, return |f9|
 208   (p8) fmerge.s f8=f0,f15
 209   (p8) br.ret.spnt b0;;
 210 }
 211
 212 { .mfi
 213          nop.m 0
 214 //   Identify Natvals, Infs, NaNs, and Zeros
 215 //   and return result
 216      fclass.m.unc p7, p0 = f12, 0x1E7
 217      nop.i 0
 218 }
 219 {.mfi
 220   nop.m 0
 221   // z0=frsqrta(a)
 222   frsqrta.s1 f8,p6=f12
 223   nop.i 0;;
 224 }
 225
 226 {.mfb
 227   // get exponent of x^2+y^2
 228   getf.exp r3=f12
 229   // if special case, set f8
 230   (p7) mov f8=f12
 231   (p7) br.ret.spnt b0;;
 232 }
 233
 234
 235 {.mfi
 236   nop.m 0
 237   // S0=a*z0
 238   (p6) fma.s1 f12=f12,f8,f0
 239   nop.i 0
 240 }
 241 {.mfi
 242   nop.m 0
 243   // H0=0.5*z0
 244   (p6) fma.s1 f10=f8,f7,f0
 245   nop.i 0;;
 246 }
 247
 248
 249 {.mfi
 250   nop.m 0
 251   // f6=5/2
 252   fma.s1 f6=f7,f1,f6
 253   nop.i 0
 254 }
 255 {.mfi
 256   nop.m 0
 257   // f11=3/2
 258   fma.s1 f11=f7,f1,f1
 259   nop.i 0;;
 260 }
 261
 262 {.mfi
 263   nop.m 0
 264   // d=0.5-S0*H0
 265   (p6) fnma.s1 f7=f12,f10,f7
 266   nop.i 0;;
 267 }
 268
 269 {.mfi
 270   nop.m 0
 271   // P01=d+1
 272   (p6) fma.s1 f10=f1,f7,f1
 273   nop.i 0
 274 }
 275 {.mfi
 276   nop.m 0
 277   // P23=5/2*d+3/2
 278   (p6) fma.s1 f11=f6,f7,f11
 279   nop.i 0;;
 280 }
 281 {.mfi
 282   nop.m 0
 283   // d2=d*d
 284   (p6) fma.s1 f7=f7,f7,f0
 285   nop.i 0;;
 286 }
 287
 288
 289 {.mfi
 290   // Is x^2 + y^2 well less than the overflow
 291   // threshold?
 292   (p6) cmp.lt.unc p7, p8 =  r3,r2
 293   // P=P01+d2*P23
 294   (p6) fma.s1 f10=f7,f11,f10
 295   nop.i 0;;
 296 }
 297
 298 {.mfb
 299   nop.m 0
 300   // S=P*S0
 301   fma.s.s0 f8=f10,f12,f0
 302   // No overflow in this case
 303   (p7) br.ret.sptk b0;;
 304 }
 305
 306 { .mfi
 307      nop.m 0
 308 (p8) fsetc.s2 0x7F,0x42
 309      // Possible overflow path, must detect by
 310      // Setting widest range exponent with prevailing
 311      // rounding mode.
 312      nop.i 0 ;;
 313 }
 314
 315
 316 { .mfi
 317    // bias+0x400 (bias+EMAX+1)
 318    (p8) mov r2=0x1007f
 319    // S=P*S0
 320    (p8) fma.s.s2 f12=f10,f12,f0
 321    nop.i 0 ;;
 322 }
 323 { .mfi
 324 (p8) setf.exp f11 = r2
 325 (p8) fsetc.s2 0x7F,0x40
 326 //   Restore Original Mode in S2
 327      nop.i 0 ;;
 328 }
 329 { .mfi
 330      nop.m 0
 331 (p8) fcmp.lt.unc.s1 p9, p10 =  f12, f11
 332      nop.i 0 ;;
 333 }
 334 { .mib
 335      nop.m 0
 336      mov GR_Parameter_TAG = 47
 337          // No overflow
 338 (p9) br.ret.sptk b0;;
 339 }
 340 GLOBAL_IEEE754_END(hypotf)
 341
 342 LOCAL_LIBM_ENTRY(__libm_error_region)
 343 .prologue
 344 { .mii
 345         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 346         mov   GR_Parameter_TAG = 47
 347 .save   ar.pfs,GR_SAVE_PFS
 348         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 349 }
 350 { .mfi
 351 .fframe 64
 352         add sp=-64,sp                           // Create new stack
 353         nop.f 0
 354         mov GR_SAVE_GP=gp                       // Save gp
 355 };;
 356 { .mmi
 357         stfs [GR_Parameter_Y] = FR_Y,16         // Store Parameter 2 on stack
 358         add GR_Parameter_X = 16,sp              // Parameter 1 address
 359 .save   b0, GR_SAVE_B0
 360         mov GR_SAVE_B0=b0                       // Save b0
 361 };;
 362 .body
 363 { .mib
 364         stfs [GR_Parameter_X] = FR_X            // Store Parameter 1 on stack
 365         add   GR_Parameter_RESULT = 0,GR_Parameter_Y
 366         nop.b 0                                 // Parameter 3 address
 367 }
 368 { .mib
 369         stfs [GR_Parameter_Y] = FR_RESULT       // Store Parameter 3 on stack
 370         add   GR_Parameter_Y = -16,GR_Parameter_Y
 371         br.call.sptk b0=__libm_error_support#   // Call error handling function
 372 };;
 373 { .mmi
 374         nop.m 0
 375         nop.m 0
 376         add   GR_Parameter_RESULT = 48,sp
 377 };;
 378 { .mmi
 379         ldfs  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 380 .restore sp
 381         add   sp = 64,sp                       // Restore stack pointer
 382         mov   b0 = GR_SAVE_B0                  // Restore return address
 383 };;
 384 { .mib
 385         mov   gp = GR_SAVE_GP                  // Restore gp
 386         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 387         br.ret.sptk     b0                     // Return
 388 };;
 389
 390 LOCAL_LIBM_END(__libm_error_region)
 391
 392
 393 .type   __libm_error_support#,@function
 394 .global __libm_error_support#
 395