sysdeps/ia64/fpu/e_sqrt.S

   1 .file "sqrt.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 // ********************************************************************
  41 // History
  42 // ********************************************************************
  43 // 2/02/00  Initial version
  44 // 4/04/00  Unwind support added
  45 // 8/15/00  Bundle added after call to __libm_error_support to properly
  46 //          set [the previously overwritten] GR_Parameter_RESULT.
  47 //
  48 // ********************************************************************
  49 //
  50 // Function:   Combined sqrt(x), where
  51 //                        _
  52 //             sqrt(x) = |x, for double precision x values
  53 //
  54 // ********************************************************************
  55 //
  56 // Accuracy:       Correctly Rounded
  57 //
  58 // ********************************************************************
  59 //
  60 // Resources Used:
  61 //
  62 //    Floating-Point Registers: f8  (Input and Return Value)
  63 //                              f7 -f14
  64 //
  65 //    General Purpose Registers:
  66 //      r32-r36 (Locals)
  67 //      r37-r40 (Used to pass arguments to error handling routine)
  68 //
  69 //    Predicate Registers:      p6, p7, p8
  70 //
  71 // *********************************************************************
  72 //
  73 // IEEE Special Conditions:
  74 //
  75 //    All faults and exceptions should be raised correctly.
  76 //    sqrt(QNaN) = QNaN
  77 //    sqrt(SNaN) = QNaN
  78 //    sqrt(+/-0) = +/-0
  79 //    sqrt(negative) = QNaN and error handling is called
  80 //
  81 // *********************************************************************
  82 //
  83 // Implementation:
  84 //
  85 //  Modified Newton-Raphson Algorithm
  86 //
  87 // *********************************************************************
  88
  89 #include "libm_support.h"
  90
  91 GR_SAVE_PFS          = r33
  92 GR_SAVE_B0           = r34
  93 GR_SAVE_GP           = r35
  94
  95 GR_Parameter_X       = r37
  96 GR_Parameter_Y       = r38
  97 GR_Parameter_RESULT  = r39
  98
  99
 100 .section .text
 101 .proc sqrt#
 102 .global sqrt#
 103 .align 64
 104
 105 sqrt:
 106 #ifdef _LIBC
 107 .global __sqrt
 108 .type __sqrt,@function
 109 __sqrt:
 110 .global __ieee754_sqrt
 111 .type __ieee754_sqrt,@function
 112 __ieee754_sqrt:
 113 #endif
 114 { .mfi
 115   alloc r32= ar.pfs,0,5,4,0
 116   frsqrta.s0 f7,p6=f8
 117   nop.i 0
 118 } { .mlx
 119   // BEGIN DOUBLE PRECISION MINIMUM LATENCY SQUARE ROOT ALGORITHM
 120   nop.m 0
 121   // exponent of +1/2 in r2
 122   movl r2 = 0x0fffe;;
 123 } { .mmi
 124   // +1/2 in f9
 125   setf.exp f9 = r2
 126   nop.m 0
 127   nop.i 0
 128 } { .mlx
 129   nop.m 0
 130   // 3/2 in r3
 131   movl r3=0x3fc00000;;
 132 } { .mfi
 133   setf.s f10=r3
 134   // Step (1)
 135   // y0 = 1/sqrt(a) in f7
 136   fclass.m.unc p7,p8 = f8,0x3A
 137   nop.i 0;;
 138 } { .mlx
 139   nop.m 0
 140   // 5/2 in r2
 141   movl r2 = 0x40200000
 142 } { .mlx
 143   nop.m 0
 144   // 63/8 in r3
 145   movl r3 = 0x40fc0000;;
 146 } { .mfi
 147   setf.s f11=r2
 148   // Step (2)
 149   // h = +1/2 * y0 in f6
 150   (p6) fma.s1 f6=f9,f7,f0
 151   nop.i 0
 152 } { .mfi
 153   setf.s f12=r3
 154   // Step (3)
 155   // g = a * y0 in f7
 156   (p6) fma.s1 f7=f8,f7,f0
 157   nop.i 0
 158 } { .mfi
 159   nop.m 0
 160   mov   f15 = f8
 161   nop.i 0;;
 162 } { .mlx
 163   nop.m 0
 164   // 231/16 in r2
 165   movl r2 = 0x41670000;;
 166 } { .mfi
 167   setf.s f13=r2
 168   // Step (4)
 169   // e = 1/2 - g * h in f9
 170   (p6) fnma.s1 f9=f7,f6,f9
 171   nop.i 0
 172 } { .mlx
 173   nop.m 0
 174   // 35/8 in r3
 175   movl r3 = 0x408c0000;;
 176 } { .mfi
 177   setf.s f14=r3
 178   // Step (5)
 179   // S = 3/2 + 5/2 * e in f10
 180   (p6) fma.s1 f10=f11,f9,f10
 181   nop.i 0
 182 } { .mfi
 183   nop.m 0
 184   // Step (6)
 185   // e2 = e * e in f11
 186   (p6) fma.s1 f11=f9,f9,f0
 187   nop.i 0;;
 188 } { .mfi
 189   nop.m 0
 190   // Step (7)
 191   // t = 63/8 + 231/16 * e in f12
 192   (p6) fma.s1 f12=f13,f9,f12
 193   nop.i 0;;
 194 } { .mfi
 195   nop.m 0
 196   // Step (8)
 197   // S1 = e + e2 * S in f10
 198   (p6) fma.s1 f10=f11,f10,f9
 199   nop.i 0
 200 } { .mfi
 201   nop.m 0
 202   // Step (9)
 203   // e4 = e2 * e2 in f11
 204   (p6) fma.s1 f11=f11,f11,f0
 205   nop.i 0;;
 206 } { .mfi
 207   nop.m 0
 208   // Step (10)
 209   // t1 = 35/8 + e * t in f9
 210   (p6) fma.s1 f9=f9,f12,f14
 211   nop.i 0;;
 212 } { .mfi
 213   nop.m 0
 214   // Step (11)
 215   // G = g + S1 * g in f12
 216   (p6) fma.s1 f12=f10,f7,f7
 217   nop.i 0
 218 } { .mfi
 219   nop.m 0
 220   // Step (12)
 221   // E = g * e4 in f7
 222   (p6) fma.s1 f7=f7,f11,f0
 223   nop.i 0;;
 224 } { .mfi
 225   nop.m 0
 226   // Step (13)
 227   // u = S1 + e4 * t1 in f10
 228   (p6) fma.s1 f10=f11,f9,f10
 229   nop.i 0;;
 230 } { .mfi
 231   nop.m 0
 232   // Step (14)
 233   // g1 = G + t1 * E in f7
 234   (p6) fma.d.s1 f7=f9,f7,f12
 235   nop.i 0;;
 236 } { .mfi
 237   nop.m 0
 238   // Step (15)
 239   // h1 = h + u * h in f6
 240   (p6) fma.s1 f6=f10,f6,f6
 241   nop.i 0;;
 242 } { .mfi
 243   nop.m 0
 244   // Step (16)
 245   // d = a - g1 * g1 in f9
 246   (p6) fnma.s1 f9=f7,f7,f8
 247   nop.i 0;;
 248 } { .mfb
 249   nop.m 0
 250   // Step (17)
 251   // g2 = g1 + d * h1 in f7
 252   (p6) fma.d.s0 f8=f9,f6,f7
 253   (p6) br.ret.sptk b0 ;;
 254 }
 255
 256 { .mfb
 257   nop.m 0
 258   (p0) mov   f8 = f7
 259   (p8) br.ret.sptk b0 ;;
 260 }
 261 { .mfb
 262   (p7) mov   r40 = 49
 263   nop.f 0
 264   (p7) br.cond.sptk __libm_error_region ;;
 265 }
 266 // END DOUBLE PRECISION MINIMUM LATENCY SQUARE ROOT ALGORITHM
 267 .endp sqrt#
 268 ASM_SIZE_DIRECTIVE(sqrt)
 269 #ifdef _LIBC
 270 ASM_SIZE_DIRECTIVE(__sqrt)
 271 ASM_SIZE_DIRECTIVE(__ieee754_sqrt)
 272 #endif
 273
 274 // Stack operations when calling error support.
 275 //       (1)               (2)                          (3) (call)              (4)
 276 //   sp   -> +          psp -> +                     psp -> +                   sp -> +
 277 //           |                 |                            |                         |
 278 //           |                 | <- GR_Y               R3 ->| <- GR_RESULT            | -> f8
 279 //           |                 |                            |                         |
 280 //           | <-GR_Y      Y2->|                       Y2 ->| <- GR_Y                 |
 281 //           |                 |                            |                         |
 282 //           |                 | <- GR_X               X1 ->|                         |
 283 //           |                 |                            |                         |
 284 //  sp-64 -> +          sp ->  +                     sp ->  +                         +
 285 //    save ar.pfs          save b0                                               restore gp
 286 //    save gp                                                                    restore ar.pfs
 287
 288
 289 .proc __libm_error_region
 290 __libm_error_region:
 291
 292 //
 293 // This branch includes all those special values that are not negative,
 294 // with the result equal to frcpa(x)
 295 //
 296
 297 .prologue
 298 // We are distinguishing between over(under)flow and letting
 299 // __libm_error_support set ERANGE or do anything else needed.
 300
 301 // (1)
 302 { .mfi
 303         add   GR_Parameter_Y=-32,sp             // Parameter 2 value
 304         nop.f 0
 305 .save   ar.pfs,GR_SAVE_PFS
 306         mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs
 307 }
 308 { .mfi
 309 .fframe 64
 310         add sp=-64,sp                          // Create new stack
 311         nop.f 0
 312         mov GR_SAVE_GP=gp                      // Save gp
 313 };;
 314
 315
 316 // (2)
 317 { .mmi
 318         stfd [GR_Parameter_Y] = f0,16         // STORE Parameter 2 on stack
 319         add GR_Parameter_X = 16,sp            // Parameter 1 address
 320 .save   b0, GR_SAVE_B0
 321         mov GR_SAVE_B0=b0                     // Save b0
 322 };;
 323
 324 .body
 325 // (3)
 326 { .mib
 327         stfd [GR_Parameter_X] = f15                    // STORE Parameter 1 on stack
 328         add   GR_Parameter_RESULT = 0,GR_Parameter_Y   // Parameter 3 address
 329         nop.b 0
 330 }
 331 { .mib
 332         stfd [GR_Parameter_Y] = f8                     // STORE Parameter 3 on stack
 333         add   GR_Parameter_Y = -16,GR_Parameter_Y
 334         br.call.sptk b0=__libm_error_support#          // Call error handling function
 335 };;
 336 { .mmi
 337         nop.m 0
 338         nop.m 0
 339         add   GR_Parameter_RESULT = 48,sp
 340 };;
 341
 342 // (4)
 343 { .mmi
 344         ldfd  f8 = [GR_Parameter_RESULT]       // Get return result off stack
 345 .restore sp
 346         add   sp = 64,sp                       // Restore stack pointer
 347         mov   b0 = GR_SAVE_B0                  // Restore return address
 348 };;
 349 { .mib
 350         mov   gp = GR_SAVE_GP                  // Restore gp
 351         mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
 352         br.ret.sptk     b0                     // Return
 353 };;
 354
 355 .endp __libm_error_region
 356 ASM_SIZE_DIRECTIVE(__libm_error_region)
 357
 358
 359 .type   __libm_error_support#,@function
 360 .global __libm_error_support#