Back to index

glibc  2.9
e_sqrt.c
Go to the documentation of this file.
00001 /* Copyright (C) 1996,1997,1998,2002,2003 Free Software Foundation, Inc.
00002    Contributed by David Mosberger (davidm@cs.arizona.edu).
00003    This file is part of the GNU C Library.
00004 
00005    The GNU C Library is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU Lesser General Public
00007    License as published by the Free Software Foundation; either
00008    version 2.1 of the License, or (at your option) any later version.
00009 
00010    The GNU C Library is distributed in the hope that it will be useful,
00011    but WITHOUT ANY WARRANTY; without even the implied warranty of
00012    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013    Lesser General Public License for more details.
00014 
00015    You should have received a copy of the GNU Lesser General Public
00016    License along with the GNU C Library; if not, write to the Free
00017    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00018    02111-1307 USA.  */
00019 
00020 #include <features.h>
00021 
00022 #if !defined(_IEEE_FP_INEXACT)
00023 
00024 /*
00025  * This version is much faster than generic sqrt implementation, but
00026  * it doesn't handle the inexact flag.  It doesn't handle exceptional
00027  * values either, but will defer to the full ieee754_sqrt routine which
00028  * can.
00029  */
00030 
00031 /* Careful with rearranging this without consulting the assembly below.  */
00032 const static struct sqrt_data_struct {
00033        unsigned long dn, up, half, almost_three_half;
00034        unsigned long one_and_a_half, two_to_minus_30, one, nan;
00035        const int T2[64];
00036 } sqrt_data __attribute__((used)) = {
00037        0x3fefffffffffffff,  /* __dn = nextafter(1,-Inf) */
00038        0x3ff0000000000001,  /* __up = nextafter(1,+Inf) */
00039        0x3fe0000000000000,  /* half */
00040        0x3ff7ffffffc00000,  /* almost_three_half = 1.5-2^-30 */
00041        0x3ff8000000000000,  /* one_and_a_half */
00042        0x3e10000000000000,  /* two_to_minus_30 */
00043        0x3ff0000000000000,  /* one */
00044        0xffffffffffffffff,  /* nan */
00045 
00046        { 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
00047        0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
00048        0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
00049        0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
00050        0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
00051        0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
00052        0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
00053        0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd }
00054 };
00055 
00056 asm ("\
00057   /* Define offsets into the structure defined in C above.  */        \n\
00058        $DN = 0*8                                               \n\
00059        $UP = 1*8                                               \n\
00060        $HALF = 2*8                                             \n\
00061        $ALMOST_THREE_HALF = 3*8                                \n\
00062        $NAN = 7*8                                              \n\
00063        $T2 = 8*8                                               \n\
00064                                                                \n\
00065   /* Stack variables.  */                                      \n\
00066        $K = 0                                                  \n\
00067        $Y = 8                                                  \n\
00068                                                                \n\
00069        .text                                                   \n\
00070        .align 5                                                \n\
00071        .globl __ieee754_sqrt                                          \n\
00072        .ent   __ieee754_sqrt                                          \n\
00073 __ieee754_sqrt:                                                       \n\
00074        ldgp   $29, 0($27)                                      \n\
00075        subq   $sp, 16, $sp                                     \n\
00076        .frame $sp, 16, $26, 0\n"
00077 #ifdef PROF
00078 "      lda    $28, _mcount                                     \n\
00079        jsr    $28, ($28), _mcount\n"
00080 #endif
00081 "      .prologue 1                                             \n\
00082                                                                \n\
00083        .align 4                                                \n\
00084        stt    $f16, $K($sp)        # e0    :                   \n\
00085        mult   $f31, $f31, $f31     # .. fm :                   \n\
00086        lda    $4, sqrt_data        # e0    :                   \n\
00087        fblt   $f16, $fixup         # .. fa :                   \n\
00088                                                                \n\
00089        ldah   $2, 0x5fe8           # e0    :                   \n\
00090        ldq    $3, $K($sp)          # .. e1 :                   \n\
00091        ldt    $f12, $HALF($4)             # e0    :                   \n\
00092        ldt    $f18, $ALMOST_THREE_HALF($4)       # .. e1 :            \n\
00093                                                                \n\
00094        sll    $3, 52, $5           # e0    :                   \n\
00095        lda    $6, 0x7fd            # .. e1 :                   \n\
00096        fnop                        # .. fa :                   \n\
00097        fnop                        # .. fm :                   \n\
00098                                                                \n\
00099        subq   $5, 1, $5            # e1    :                   \n\
00100        srl    $3, 33, $1           # .. e0 :                   \n\
00101        cmpule $5, $6, $5           # e0    :                   \n\
00102        beq    $5, $fixup           # .. e1 :                   \n\
00103                                                                \n\
00104        mult   $f16, $f12, $f11     # fm    : $f11 = x * 0.5    \n\
00105        subl   $2, $1, $2           # .. e0 :                   \n\
00106        addt   $f12, $f12, $f17     # .. fa : $f17 = 1.0        \n\
00107        srl    $2, 12, $1           # e0    :                   \n\
00108                                                                \n\
00109        and    $1, 0xfc, $1         # e0    :                   \n\
00110        addq   $1, $4, $1           # e1    :                   \n\
00111        ldl    $1, $T2($1)          # e0    :                   \n\
00112        addt   $f12, $f17, $f15     # .. fa : $f15 = 1.5        \n\
00113                                                                \n\
00114        subl   $2, $1, $2           # e0    :                   \n\
00115        ldt    $f14, $DN($4)        # .. e1 :                   \n\
00116        sll    $2, 32, $2           # e0    :                   \n\
00117        stq    $2, $Y($sp)          # e0    :                   \n\
00118                                                                \n\
00119        ldt    $f13, $Y($sp)        # e0    :                   \n\
00120        mult/su       $f11, $f13, $f10     # fm   2: $f10 = (x * 0.5) * y     \n\
00121        mult   $f10, $f13, $f10     # fm   4: $f10 = ((x*0.5)*y)*y     \n\
00122        subt   $f15, $f10, $f1             # fa   4: $f1 = (1.5-0.5*x*y*y)    \n\
00123                                                                \n\
00124        mult   $f13, $f1, $f13         # fm   4: yp = y*(1.5-0.5*x*y^2)\n\
00125        mult/su       $f11, $f13, $f1             # fm   4: $f11 = x * 0.5 * yp      \n\
00126        mult   $f1, $f13, $f11             # fm   4: $f11 = (x*0.5*yp)*yp     \n\
00127        subt   $f18, $f11, $f1             # fa   4: $f1=(1.5-2^-30)-x/2*yp^2\n\
00128                                                                \n\
00129        mult   $f13, $f1, $f13             # fm   4: ypp = $f13 = yp*$f1      \n\
00130        subt   $f15, $f12, $f1             # .. fa : $f1 = (1.5 - 0.5) \n\
00131        ldt    $f15, $UP($4)        # .. e0 :                   \n\
00132        mult/su       $f16, $f13, $f10     # fm   4: z = $f10 = x * ypp       \n\
00133                                                                \n\
00134        mult   $f10, $f13, $f11     # fm   4: $f11 = z*ypp             \n\
00135        mult   $f10, $f12, $f12     # fm    : $f12 = z*0.5             \n\
00136        subt   $f1, $f11, $f1              # fa   4: $f1 = 1 - z*ypp   \n\
00137        mult   $f12, $f1, $f12             # fm   4: $f12 = z/2*(1 - z*ypp)\n\
00138                                                                \n\
00139        addt   $f10, $f12, $f0             # fa   4: zp=res= z+z/2*(1-z*ypp)\n\
00140        mult/c $f0, $f14, $f12             # fm   4: zmi = zp * DN            \n\
00141        mult/c $f0, $f15, $f11             # fm    : zpl = zp * UP            \n\
00142        mult/c $f0, $f12, $f1              # fm    : $f1 = zp * zmi    \n\
00143                                                                \n\
00144        mult/c $f0, $f11, $f15             # fm    : $f15 = zp * zpl   \n\
00145        subt/su       $f1, $f16, $f13             # .. fa : y1 = zp*zmi - x   \n\
00146        subt/su       $f15, $f16, $f14     # fa   4: y2 = zp*zpl - x   \n\
00147        fcmovge       $f13, $f12, $f0             # fa   3: res = (y1>=0)?zmi:res    \n\
00148                                                                \n\
00149        fcmovlt       $f14, $f11, $f0             # fa   4: res = (y2<0)?zpl:res     \n\
00150        addq   $sp, 16, $sp         # .. e0 :                   \n\
00151        ret                         # .. e1 :                   \n\
00152                                                                \n\
00153        .align 4                                                \n\
00154 $fixup:                                                               \n\
00155        addq   $sp, 16, $sp                                     \n\
00156        br     __full_ieee754_sqrt  !samegp                            \n\
00157                                                                \n\
00158        .end   __ieee754_sqrt");
00159 
00160 static double __full_ieee754_sqrt(double) __attribute_used__;
00161 #define __ieee754_sqrt __full_ieee754_sqrt
00162 
00163 #endif /* _IEEE_FP_INEXACT */
00164 
00165 #include <sysdeps/ieee754/dbl-64/e_sqrt.c>