update from main archive 961119
[kopensolaris-gnu/glibc.git] / sysdeps / alpha / w_sqrt.S
1 /* Copyright (C) 1996 Free Software Foundation, Inc.
2    This file is part of the GNU C Library.
3    Contributed by David Mosberger <davidm@cs.arizona.edu>, 1996.
4    Based on public-domain C source by Linus Torvalds.
5
6    The GNU C Library is free software; you can redistribute it and/or
7    modify it under the terms of the GNU Library General Public License as
8    published by the Free Software Foundation; either version 2 of the
9    License, or (at your option) any later version.
10
11    The GNU C Library is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14    Library General Public License for more details.
15
16    You should have received a copy of the GNU Library General Public
17    License along with the GNU C Library; see the file COPYING.LIB.  If not,
18    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
19    Boston, MA 02111-1307, USA.  */
20
21 /* This version is much faster than generic sqrt implementation, but
22    it doesn't handle exceptional values or the inexact flag.  Don't use
23    this if _IEEE_FP or _IEEE_FP_INEXACT is in effect. */
24
25 #ifndef _IEEE_FP
26
27 #define _ERRNO_H
28 #include <errnos.h>
29 #include <sysdep.h>
30
31         .set noreorder
32
33 #ifdef __ELF__
34         .section .rodata
35 #else
36         .rdata
37 #endif
38         .align 5        # align to cache line
39
40         /* Do all memory accesses relative to sqrtdata.  */
41 sqrtdata:
42
43 #define DN                     0x00
44 #define UP                     0x08
45 #define HALF                   0x10
46 #define ALMOST_THREE_HALF      0x18
47 #define T2                     0x20
48
49         .quad 0x3fefffffffffffff        /* DN = next(1.0) */
50         .quad 0x3ff0000000000001        /* UP = prev(1.0) */
51         .quad 0x3fe0000000000000        /* HALF = 0.5 */
52         .quad 0x3ff7ffffffc00000        /* ALMOST_THREE_HALF = 1.5-2^-30 */
53
54 /* table T2: */
55 .long   0x1500, 0x2ef8,   0x4d67,  0x6b02,  0x87be,  0xa395,  0xbe7a,  0xd866
56 .long   0xf14a, 0x1091b, 0x11fcd, 0x13552, 0x14999, 0x15c98, 0x16e34, 0x17e5f
57 .long  0x18d03, 0x19a01, 0x1a545, 0x1ae8a, 0x1b5c4, 0x1bb01, 0x1bfde, 0x1c28d
58 .long  0x1c2de, 0x1c0db, 0x1ba73, 0x1b11c, 0x1a4b5, 0x1953d, 0x18266, 0x16be0
59 .long  0x1683e, 0x179d8, 0x18a4d, 0x19992, 0x1a789, 0x1b445, 0x1bf61, 0x1c989
60 .long  0x1d16d, 0x1d77b, 0x1dddf, 0x1e2ad, 0x1e5bf, 0x1e6e8, 0x1e654, 0x1e3cd
61 .long  0x1df2a, 0x1d635, 0x1cb16, 0x1be2c, 0x1ae4e, 0x19bde, 0x1868e, 0x16e2e
62 .long  0x1527f, 0x1334a, 0x11051,  0xe951,  0xbe01,  0x8e0d,  0x5924,  0x1edd
63
64 /*
65  * Stack variables:
66  */
67 #define K      16(sp)
68 #define Y      24(sp)
69 #define FSIZE  32
70
71         .text
72
73 LEAF(__sqrt, FSIZE)
74         lda     sp, -FSIZE(sp)
75         ldgp    gp, .-__sqrt(pv)
76         stq     ra, 0(sp)
77 #ifdef PROF
78         lda     AT, _mcount
79         jsr     AT, (AT), _mcount
80 #endif
81         .prologue 1
82
83         stt     $f16, K
84         lda     t3, sqrtdata                    # load base address into t3
85
86         fblt    $f16, $negative
87
88         /* Compute initial guess.  */
89
90         .align 3
91
92         ldah    t1, 0x5fe8                      # e0    :
93         ldq     t2, K                           # .. e1 :
94         ldt     $f12, HALF(t3)                  # e0    :
95         ldt     $f18, ALMOST_THREE_HALF(t3)     # .. e1 :
96         srl     t2, 33, t0                      # e0    :
97         mult    $f16, $f12, $f11                # .. fm : $f11 = x * 0.5
98         subl    t1, t0, t1                      # e0    :
99         addt    $f12, $f12, $f17                # .. fa : $f17 = 1.0
100         srl     t1, 12, t0                      # e0    :
101         and     t0, 0xfc, t0                    # .. e1 :
102         addq    t0, t3, t0                      # e0    :
103         ldl     t0, T2(t0)                      # .. e1 :
104         addt    $f12, $f17, $f15                # fa    : $f15 = 1.5
105         subl    t1, t0, t1                      # .. e1 :
106         sll     t1, 32, t1                      # e0    :
107         ldt     $f14, DN(t3)                    # .. e1 :
108         stq     t1, Y                           # e0    :
109         ldt     $f13, Y                         # e1    :
110         addq    sp, FSIZE, sp                   # e0    :
111
112         mult    $f11, $f13, $f10        # fm    : $f10 = (x * 0.5) * y
113         mult    $f10, $f13, $f10        # fm    : $f10 = ((x * 0.5) * y) * y
114         subt    $f15, $f10, $f1         # fa    : $f1 = (1.5 - 0.5*x*y*y)
115         mult    $f13, $f1, $f13         # fm    : yp = y*(1.5 - 0.5*x*y*y)
116         mult    $f11, $f13, $f11        # fm    : $f11 = x * 0.5 * yp
117         mult    $f11, $f13, $f11        # fm    : $f11 = (x * 0.5 * yp) * yp
118         subt    $f18, $f11, $f1         # fa    : $f1= (1.5-2^-30) - 0.5*x*yp*yp
119         mult    $f13, $f1, $f13         # fm    : ypp = $f13 = yp*$f1
120         subt    $f15, $f12, $f1         # fa    : $f1 = (1.5 - 0.5)
121         ldt     $f15, UP(t3)            # .. e1 :
122         mult    $f16, $f13, $f10        # fm    : z = $f10 = x * ypp
123         mult    $f10, $f13, $f11        # fm    : $f11 = z*ypp
124         mult    $f10, $f12, $f12        # fm    : $f12 = z*0.5
125         subt    $f1, $f11, $f1          # .. fa : $f1 = 1 - z*ypp
126         mult    $f12, $f1, $f12         # fm    : $f12 = z*0.5*(1 - z*ypp)
127         addt    $f10, $f12, $f0         # fa    : zp=res=$f0= z + z*0.5*(1 - z*ypp)
128
129         mult/c  $f0, $f14, $f12         # fm    : zmi = zp * DN
130         mult/c  $f0, $f15, $f11         # fm    : zpl = zp * UP
131         mult/c  $f0, $f12, $f1          # fm    : $f1 = zp * zmi
132         mult/c  $f0, $f11, $f15         # fm    : $f15 = zp * zpl
133
134         subt    $f1, $f16, $f13         # fa    : y1 = zp*zmi - x
135         subt    $f15, $f16, $f15        # fa    : y2 = zp*zpl - x
136
137         fcmovge $f13, $f12, $f0         # res = (y1 >= 0) ? zmi : res
138         fcmovlt $f15, $f11, $f0         # res = (y2 <  0) ? zpl : res
139
140         ret
141
142 $negative:
143         lda     t1, -1
144         stq     t1, K
145         lda     t1, EDOM
146         stl     t1, errno
147 #ifdef _LIBC_REENTRANT
148         jsr     ra, __errno_location
149         lda     t1, -1
150         ldq     ra, 0(sp)
151         stl     t1, 0(v0)
152 #endif
153         ldt     $f0, K                  # res = (double) 0xffffffffffffffff
154         addq    sp, FSIZE, sp
155         ret
156
157         END(__sqrt)
158
159 weak_alias(__sqrt, sqrt)
160
161 #endif /* !_IEEE_FP */