diff options
author | Linus Torvalds <torvalds@ppc970.osdl.org> | 2005-04-16 18:20:36 -0400 |
---|---|---|
committer | Linus Torvalds <torvalds@ppc970.osdl.org> | 2005-04-16 18:20:36 -0400 |
commit | 1da177e4c3f41524e886b7f1b8a0c1fc7321cac2 (patch) | |
tree | 0bba044c4ce775e45a88a51686b5d9f90697ea9d /arch/parisc/math-emu/dfsqrt.c |
Linux-2.6.12-rc2v2.6.12-rc2
Initial git repository build. I'm not bothering with the full history,
even though we have it. We can create a separate "historical" git
archive of that later if we want to, and in the meantime it's about
3.2GB when imported into git - space that would just make the early
git days unnecessarily complicated, when we don't have a lot of good
infrastructure for it.
Let it rip!
Diffstat (limited to 'arch/parisc/math-emu/dfsqrt.c')
-rw-r--r-- | arch/parisc/math-emu/dfsqrt.c | 195 |
1 files changed, 195 insertions, 0 deletions
diff --git a/arch/parisc/math-emu/dfsqrt.c b/arch/parisc/math-emu/dfsqrt.c new file mode 100644 index 000000000000..b6ed1066f1e4 --- /dev/null +++ b/arch/parisc/math-emu/dfsqrt.c | |||
@@ -0,0 +1,195 @@ | |||
1 | /* | ||
2 | * Linux/PA-RISC Project (http://www.parisc-linux.org/) | ||
3 | * | ||
4 | * Floating-point emulation code | ||
5 | * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org> | ||
6 | * | ||
7 | * This program is free software; you can redistribute it and/or modify | ||
8 | * it under the terms of the GNU General Public License as published by | ||
9 | * the Free Software Foundation; either version 2, or (at your option) | ||
10 | * any later version. | ||
11 | * | ||
12 | * This program is distributed in the hope that it will be useful, | ||
13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
15 | * GNU General Public License for more details. | ||
16 | * | ||
17 | * You should have received a copy of the GNU General Public License | ||
18 | * along with this program; if not, write to the Free Software | ||
19 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | ||
20 | */ | ||
21 | /* | ||
22 | * BEGIN_DESC | ||
23 | * | ||
24 | * File: | ||
25 | * @(#) pa/spmath/dfsqrt.c $Revision: 1.1 $ | ||
26 | * | ||
27 | * Purpose: | ||
28 | * Double Floating-point Square Root | ||
29 | * | ||
30 | * External Interfaces: | ||
31 | * dbl_fsqrt(srcptr,nullptr,dstptr,status) | ||
32 | * | ||
33 | * Internal Interfaces: | ||
34 | * | ||
35 | * Theory: | ||
36 | * <<please update with a overview of the operation of this file>> | ||
37 | * | ||
38 | * END_DESC | ||
39 | */ | ||
40 | |||
41 | |||
42 | #include "float.h" | ||
43 | #include "dbl_float.h" | ||
44 | |||
45 | /* | ||
46 | * Double Floating-point Square Root | ||
47 | */ | ||
48 | |||
49 | /*ARGSUSED*/ | ||
50 | unsigned int | ||
51 | dbl_fsqrt( | ||
52 | dbl_floating_point *srcptr, | ||
53 | unsigned int *nullptr, | ||
54 | dbl_floating_point *dstptr, | ||
55 | unsigned int *status) | ||
56 | { | ||
57 | register unsigned int srcp1, srcp2, resultp1, resultp2; | ||
58 | register unsigned int newbitp1, newbitp2, sump1, sump2; | ||
59 | register int src_exponent; | ||
60 | register boolean guardbit = FALSE, even_exponent; | ||
61 | |||
62 | Dbl_copyfromptr(srcptr,srcp1,srcp2); | ||
63 | /* | ||
64 | * check source operand for NaN or infinity | ||
65 | */ | ||
66 | if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) { | ||
67 | /* | ||
68 | * is signaling NaN? | ||
69 | */ | ||
70 | if (Dbl_isone_signaling(srcp1)) { | ||
71 | /* trap if INVALIDTRAP enabled */ | ||
72 | if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); | ||
73 | /* make NaN quiet */ | ||
74 | Set_invalidflag(); | ||
75 | Dbl_set_quiet(srcp1); | ||
76 | } | ||
77 | /* | ||
78 | * Return quiet NaN or positive infinity. | ||
79 | * Fall thru to negative test if negative infinity. | ||
80 | */ | ||
81 | if (Dbl_iszero_sign(srcp1) || | ||
82 | Dbl_isnotzero_mantissa(srcp1,srcp2)) { | ||
83 | Dbl_copytoptr(srcp1,srcp2,dstptr); | ||
84 | return(NOEXCEPTION); | ||
85 | } | ||
86 | } | ||
87 | |||
88 | /* | ||
89 | * check for zero source operand | ||
90 | */ | ||
91 | if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) { | ||
92 | Dbl_copytoptr(srcp1,srcp2,dstptr); | ||
93 | return(NOEXCEPTION); | ||
94 | } | ||
95 | |||
96 | /* | ||
97 | * check for negative source operand | ||
98 | */ | ||
99 | if (Dbl_isone_sign(srcp1)) { | ||
100 | /* trap if INVALIDTRAP enabled */ | ||
101 | if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); | ||
102 | /* make NaN quiet */ | ||
103 | Set_invalidflag(); | ||
104 | Dbl_makequietnan(srcp1,srcp2); | ||
105 | Dbl_copytoptr(srcp1,srcp2,dstptr); | ||
106 | return(NOEXCEPTION); | ||
107 | } | ||
108 | |||
109 | /* | ||
110 | * Generate result | ||
111 | */ | ||
112 | if (src_exponent > 0) { | ||
113 | even_exponent = Dbl_hidden(srcp1); | ||
114 | Dbl_clear_signexponent_set_hidden(srcp1); | ||
115 | } | ||
116 | else { | ||
117 | /* normalize operand */ | ||
118 | Dbl_clear_signexponent(srcp1); | ||
119 | src_exponent++; | ||
120 | Dbl_normalize(srcp1,srcp2,src_exponent); | ||
121 | even_exponent = src_exponent & 1; | ||
122 | } | ||
123 | if (even_exponent) { | ||
124 | /* exponent is even */ | ||
125 | /* Add comment here. Explain why odd exponent needs correction */ | ||
126 | Dbl_leftshiftby1(srcp1,srcp2); | ||
127 | } | ||
128 | /* | ||
129 | * Add comment here. Explain following algorithm. | ||
130 | * | ||
131 | * Trust me, it works. | ||
132 | * | ||
133 | */ | ||
134 | Dbl_setzero(resultp1,resultp2); | ||
135 | Dbl_allp1(newbitp1) = 1 << (DBL_P - 32); | ||
136 | Dbl_setzero_mantissap2(newbitp2); | ||
137 | while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) { | ||
138 | Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2); | ||
139 | if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) { | ||
140 | Dbl_leftshiftby1(newbitp1,newbitp2); | ||
141 | /* update result */ | ||
142 | Dbl_addition(resultp1,resultp2,newbitp1,newbitp2, | ||
143 | resultp1,resultp2); | ||
144 | Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2); | ||
145 | Dbl_rightshiftby2(newbitp1,newbitp2); | ||
146 | } | ||
147 | else { | ||
148 | Dbl_rightshiftby1(newbitp1,newbitp2); | ||
149 | } | ||
150 | Dbl_leftshiftby1(srcp1,srcp2); | ||
151 | } | ||
152 | /* correct exponent for pre-shift */ | ||
153 | if (even_exponent) { | ||
154 | Dbl_rightshiftby1(resultp1,resultp2); | ||
155 | } | ||
156 | |||
157 | /* check for inexact */ | ||
158 | if (Dbl_isnotzero(srcp1,srcp2)) { | ||
159 | if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) { | ||
160 | Dbl_increment(resultp1,resultp2); | ||
161 | } | ||
162 | guardbit = Dbl_lowmantissap2(resultp2); | ||
163 | Dbl_rightshiftby1(resultp1,resultp2); | ||
164 | |||
165 | /* now round result */ | ||
166 | switch (Rounding_mode()) { | ||
167 | case ROUNDPLUS: | ||
168 | Dbl_increment(resultp1,resultp2); | ||
169 | break; | ||
170 | case ROUNDNEAREST: | ||
171 | /* stickybit is always true, so guardbit | ||
172 | * is enough to determine rounding */ | ||
173 | if (guardbit) { | ||
174 | Dbl_increment(resultp1,resultp2); | ||
175 | } | ||
176 | break; | ||
177 | } | ||
178 | /* increment result exponent by 1 if mantissa overflowed */ | ||
179 | if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2; | ||
180 | |||
181 | if (Is_inexacttrap_enabled()) { | ||
182 | Dbl_set_exponent(resultp1, | ||
183 | ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); | ||
184 | Dbl_copytoptr(resultp1,resultp2,dstptr); | ||
185 | return(INEXACTEXCEPTION); | ||
186 | } | ||
187 | else Set_inexactflag(); | ||
188 | } | ||
189 | else { | ||
190 | Dbl_rightshiftby1(resultp1,resultp2); | ||
191 | } | ||
192 | Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); | ||
193 | Dbl_copytoptr(resultp1,resultp2,dstptr); | ||
194 | return(NOEXCEPTION); | ||
195 | } | ||