sysdeps/ia64/fpu/e_hypot.S

   1 .file "hypot.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 // 04/17/03 Added missing mutex directive
  51 //
  52 //*********************************************************************
  53 //                           ___________
  54 // Function:   hypot(x,y) = |(x^2 + y^2) = for double precision values
  55 //             x and y
  56 //             Also provides cabs functionality.
  57 //
  58 //*********************************************************************
  59 //
  60 // Resources Used:
  61 //
  62 //    Floating-Point Registers: f8  (Input and Return Value)
  63 //                              f9  (Input)
  64 //                              f6 -f15, f32-f34
  65 //
  66 //    General Purpose Registers:
  67 //      r2,r3,r29 (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 //    hypot(Infinity and anything) = +Infinity
  80 //    hypot(QNaN and anything) = QNaN
  81 //    hypot(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 double
  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                = f32
 102 FR_Y                = f33
 103 FR_RESULT           = f8
 104
 105 .section .text
 106
 107 LOCAL_LIBM_ENTRY(cabs)
 108 LOCAL_LIBM_END(cabs)
 109
 110 GLOBAL_IEEE754_ENTRY(hypot)
 111
 112 {.mfi
 113   alloc r32= ar.pfs,0,4,4,0
 114   // Compute x*x
 115   fma.s1 f10=f8,f8,f0
 116   // r2=bias-1
 117   mov r2=0xfffe
 118 }
 119 {.mfi
 120   // 63/8
 121   mov r3=0x40fc //0000
 122   // y*y
 123   fma.s1 f11=f9,f9,f0
 124   // r29=429/16
 125   mov r29=0x41d68;; //000
 126 }
 127
 128 { .mfi
 129      nop.m 0
 130 //   Check if x is an Inf - if so return Inf even
 131 //   if y is a NaN (C9X)
 132      fclass.m.unc p7, p6 = f8, 0x023
 133      shl r3=r3,16
 134 }
 135 {.mfi
 136         nop.m 0
 137   // if possible overflow, copy f8 to f32
 138   // set Denormal, if necessary
 139   // (p8)
 140   fma.d.s0 f32=f8,f1,f0
 141   nop.i 0;;
 142 }
 143 { .mfi
 144      nop.m 0
 145 //   Check if y is an Inf - if so return Inf even
 146 //   if x is a NaN (C9X)
 147      fclass.m.unc p8, p9 = f9, 0x023
 148      shl r29=r29,12
 149 }
 150 { .mfb
 151          // f7=0.5
 152          setf.exp f7=r2
 153 //   For x=inf, multiply y by 1 to raise invalid on y an SNaN
 154 //   (p7) fma.s0 f9=f9,f1,f0
 155      // copy f9 to f33; set Denormal, if necessary
 156          fma.d.s0 f33=f9,f1,f0
 157      nop.b 0;;
 158 }
 159 {.mfb
 160   // f13=63/8
 161   setf.s f13=r3
 162   // is y Zero ?
 163   (p6) fclass.m p6,p0=f9,0x7
 164   nop.b 0
 165 }
 166 {.mlx
 167   nop.m 0
 168   movl r2=0x408c0000;;
 169 }
 170
 171 {.mfi
 172   // f34=429/16
 173   setf.s f34=r29
 174   // is x Zero ?
 175   (p9) fclass.m p9,p0=f8,0x7
 176   // 231/16
 177   mov r3=0x4167;; //0000
 178 }
 179 {.mfi
 180   nop.m 0
 181   // a=x2+y2
 182   fma.s1 f12=f10,f1,f11
 183   nop.i 0;;
 184 }
 185 {.mfi
 186   nop.m 0
 187   // y not NaN ?
 188   (p9) fclass.m p8,p0=f9,0x3f
 189   shl r3=r3,16
 190 }
 191 {.mfi
 192   nop.m 0
 193   // f6=2
 194   fma.s1 f6=f1,f1,f1
 195   nop.i 0;;
 196 }
 197
 198
 199 {.mfi
 200   nop.m 0
 201   // x not NaN ?
 202   (p6) fclass.m p7,p0=f8,0x3f
 203   nop.i 0;;
 204 }
 205 {.mfi
 206   // f9=35/8
 207   setf.s f9=r2
 208   nop.f 0
 209   // 2*emax-2
 210   mov r2=0x107fb;;
 211 }
 212
 213 .pred.rel "mutex",p7,p8
 214 {.mfb
 215   nop.m 0
 216   // if f8=Infinity or f9=Zero, return |f8|
 217   (p7) fmerge.s f8=f0,f32
 218   (p7) br.ret.spnt b0
 219 }
 220 {.mfb
 221   nop.m 0
 222   // if f9=Infinity or f8=Zero, return |f9|
 223   (p8) fmerge.s f8=f0,f33
 224   (p8) br.ret.spnt b0;;
 225 }
 226
 227
 228 {.mfi
 229   // f10 =231/16
 230   setf.s f10=r3
 231   // z0=frsqrta(a)
 232   frsqrta.s1 f8,p6=f12
 233   nop.i 0;;
 234 }
 235
 236 { .mfi
 237          nop.m 0
 238 //   Identify Natvals, Infs, NaNs, and Zeros
 239 //   and return result
 240      fclass.m.unc p7, p0 = f12, 0x1E7
 241      nop.i 0;;
 242 }
 243 {.mfb
 244   // get exponent of x^2+y^2
 245   getf.exp r3=f12
 246   // if special case, set f8
 247   (p7) mov f8=f12
 248   (p7) br.ret.spnt b0;;
 249 }
 250
 251
 252 {.mfi
 253   nop.m 0
 254   // S0=a*z0
 255   (p6) fma.s1 f14=f12,f8,f0
 256   nop.i 0
 257 }
 258 {.mfi
 259   nop.m 0
 260   // H0=0.5*z0
 261   (p6) fma.s1 f15=f8,f7,f0
 262   nop.i 0;;
 263 }
 264
 265
 266 {.mfi
 267   nop.m 0
 268   // f6=5/2
 269   fma.s1 f6=f7,f1,f6
 270   nop.i 0
 271 }
 272 {.mfi
 273   nop.m 0
 274   // f11=3/2
 275   fma.s1 f11=f7,f1,f1
 276   nop.i 0;;
 277 }
 278
 279 {.mfi
 280   nop.m 0
 281   // d=0.5-S0*H0
 282   (p6) fnma.s1 f7=f14,f15,f7
 283   nop.i 0;;
 284 }
 285
 286 {.mfi
 287   nop.m 0
 288   // P67=231/16+429/16*d
 289   (p6) fma.s1 f10=f34,f7,f10
 290   nop.i 0
 291 }
 292 {.mfi
 293   nop.m 0
 294   // P45=63/8*d+35/8
 295   (p6) fma.s1 f9=f13,f7,f9
 296   nop.i 0;;
 297 }
 298 {.mfi
 299   nop.m 0
 300   // P23=5/2*d+3/2
 301   (p6) fma.s1 f11=f6,f7,f11
 302   nop.i 0
 303 }
 304 {.mfi
 305   nop.m 0
 306   // d2=d*d
 307   (p6) fma.s1 f13=f7,f7,f0
 308   nop.i 0;;
 309 }
 310
 311 {.mfi
 312   nop.m 0
 313   // P47=d2*P67+P45
 314   (p6) fma.s1 f10=f10,f13,f9
 315   nop.i 0
 316 }
 317 {.mfi
 318   nop.m 0
 319   // P13=d*P23+1
 320   (p6) fma.s1 f11=f11,f7,f1
 321   nop.i 0;;
 322 }
 323 {.mfi
 324   nop.m 0
 325   // d3=d2*d
 326   (p6) fma.s1 f13=f13,f7,f0
 327   nop.i 0;;
 328 }
 329
 330 {.mfi
 331   nop.m 0
 332   // T0=d*S0
 333   (p6) fma.s1 f15=f7,f14,f0
 334   nop.i 0
 335 }
 336 {.mfi
 337   // Is x^2 + y^2 well less than the overflow
 338   // threshold?
 339   (p6) cmp.lt.unc p7, p8 =  r3,r2
 340   // P=P13+d3*P47
 341   (p6) fma.s1 f10=f13,f10,f11
 342   nop.i 0;;
 343 }
 344
 345 {.mfb
 346   nop.m 0
 347   // S=P*T0+S0
 348   fma.d.s0 f8=f10,f15,f14
 349   // No overflow in this case
 350   (p7) br.ret.sptk b0;;
 351 }
 352
 353 { .mfi
 354      nop.m 0
 355 (p8) fsetc.s2 0x7F,0x42
 356      // Possible overflow path, must detect by
 357      // Setting widest range exponent with prevailing
 358      // rounding mode.
 359      nop.i 0 ;;
 360 }
 361
 362
 363 { .mfi
 364    // bias+0x400 (bias+EMAX+1)
 365    (p8) mov r2=0x103ff
 366    // S=P*T0+S0
 367    (p8) fma.d.s2 f12=f10,f15,f14
 368    nop.i 0 ;;
 369 }
 370 { .mfi
 371 (p8) setf.exp f11 = r2
 372 (p8) fsetc.s2 0x7F,0x40
 373 //   Restore Original Mode in S2
 374      nop.i 0 ;;
 375 }
 376 { .mfi
 377      nop.m 0
 378 (p8) fcmp.lt.unc.s1 p9, p10 =  f12, f11
 379      nop.i 0 ;;
 380 }
 381 { .mib
 382      nop.m 0
 383      mov GR_Parameter_TAG = 46
 384          // No overflow
 385 (p9) br.ret.sptk b0;;
 386 }
 387 GLOBAL_IEEE754_END(hypot)
 388
 389 LOCAL_LIBM_ENTRY(__libm_error_region)
 390 .prologue
 391 { .mfi
 392         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 393         nop.f 0
 394 .save   ar.pfs,GR_SAVE_PFS
 395         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 396 }
 397 { .mfi
 398 .fframe 64
 399         add sp=-64,sp                           // Create new stack
 400         nop.f 0
 401         mov GR_SAVE_GP=gp                       // Save gp
 402 };;
 403 { .mmi
 404         stfd [GR_Parameter_Y] = FR_Y,16         // Save Parameter 2 on stack
 405         add GR_Parameter_X = 16,sp              // Parameter 1 address
 406 .save   b0, GR_SAVE_B0
 407         mov GR_SAVE_B0=b0                       // Save b0
 408 };;
 409 .body
 410 { .mib
 411         stfd [GR_Parameter_X] = FR_X            // Store Parameter 1 on stack
 412         add   GR_Parameter_RESULT = 0,GR_Parameter_Y
 413         nop.b 0                                 // Parameter 3 address
 414 }
 415 { .mib
 416         stfd [GR_Parameter_Y] = FR_RESULT      // Store Parameter 3 on stack
 417         add   GR_Parameter_Y = -16,GR_Parameter_Y
 418         br.call.sptk b0=__libm_error_support#  // Call error handling function
 419 };;
 420 { .mmi
 421         nop.m 0
 422         nop.m 0
 423         add   GR_Parameter_RESULT = 48,sp
 424 };;
 425 { .mmi
 426         ldfd  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 427 .restore sp
 428         add   sp = 64,sp                       // Restore stack pointer
 429         mov   b0 = GR_SAVE_B0                  // Restore return address
 430 };;
 431 { .mib
 432         mov   gp = GR_SAVE_GP                  // Restore gp
 433         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 434         br.ret.sptk     b0                     // Return
 435 };;
 436 LOCAL_LIBM_END(__libm_error_region#)
 437 .type   __libm_error_support#,@function
 438 .global __libm_error_support#
 439
 440