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/mips/math-emu/dp_mul.c |
Linux-2.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/mips/math-emu/dp_mul.c')
-rw-r--r-- | arch/mips/math-emu/dp_mul.c | 177 |
1 files changed, 177 insertions, 0 deletions
diff --git a/arch/mips/math-emu/dp_mul.c b/arch/mips/math-emu/dp_mul.c new file mode 100644 index 00000000000..f2373902f52 --- /dev/null +++ b/arch/mips/math-emu/dp_mul.c | |||
@@ -0,0 +1,177 @@ | |||
1 | /* IEEE754 floating point arithmetic | ||
2 | * double precision: common utilities | ||
3 | */ | ||
4 | /* | ||
5 | * MIPS floating point support | ||
6 | * Copyright (C) 1994-2000 Algorithmics Ltd. | ||
7 | * http://www.algor.co.uk | ||
8 | * | ||
9 | * ######################################################################## | ||
10 | * | ||
11 | * This program is free software; you can distribute it and/or modify it | ||
12 | * under the terms of the GNU General Public License (Version 2) as | ||
13 | * published by the Free Software Foundation. | ||
14 | * | ||
15 | * This program is distributed in the hope it will be useful, but WITHOUT | ||
16 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | ||
17 | * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | ||
18 | * for more details. | ||
19 | * | ||
20 | * You should have received a copy of the GNU General Public License along | ||
21 | * with this program; if not, write to the Free Software Foundation, Inc., | ||
22 | * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. | ||
23 | * | ||
24 | * ######################################################################## | ||
25 | */ | ||
26 | |||
27 | |||
28 | #include "ieee754dp.h" | ||
29 | |||
30 | ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y) | ||
31 | { | ||
32 | COMPXDP; | ||
33 | COMPYDP; | ||
34 | |||
35 | EXPLODEXDP; | ||
36 | EXPLODEYDP; | ||
37 | |||
38 | CLEARCX; | ||
39 | |||
40 | FLUSHXDP; | ||
41 | FLUSHYDP; | ||
42 | |||
43 | switch (CLPAIR(xc, yc)) { | ||
44 | case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): | ||
45 | case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): | ||
46 | case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN): | ||
47 | case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN): | ||
48 | case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN): | ||
49 | case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN): | ||
50 | case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN): | ||
51 | case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO): | ||
52 | case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM): | ||
53 | case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM): | ||
54 | case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF): | ||
55 | SETCX(IEEE754_INVALID_OPERATION); | ||
56 | return ieee754dp_nanxcpt(ieee754dp_indef(), "mul", x, y); | ||
57 | |||
58 | case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN): | ||
59 | case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN): | ||
60 | case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN): | ||
61 | case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN): | ||
62 | return y; | ||
63 | |||
64 | case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN): | ||
65 | case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO): | ||
66 | case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM): | ||
67 | case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM): | ||
68 | case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF): | ||
69 | return x; | ||
70 | |||
71 | |||
72 | /* Infinity handling */ | ||
73 | |||
74 | case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO): | ||
75 | case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF): | ||
76 | SETCX(IEEE754_INVALID_OPERATION); | ||
77 | return ieee754dp_xcpt(ieee754dp_indef(), "mul", x, y); | ||
78 | |||
79 | case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF): | ||
80 | case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF): | ||
81 | case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM): | ||
82 | case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM): | ||
83 | case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF): | ||
84 | return ieee754dp_inf(xs ^ ys); | ||
85 | |||
86 | case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): | ||
87 | case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM): | ||
88 | case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM): | ||
89 | case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO): | ||
90 | case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO): | ||
91 | return ieee754dp_zero(xs ^ ys); | ||
92 | |||
93 | |||
94 | case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): | ||
95 | DPDNORMX; | ||
96 | |||
97 | case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): | ||
98 | DPDNORMY; | ||
99 | break; | ||
100 | |||
101 | case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): | ||
102 | DPDNORMX; | ||
103 | break; | ||
104 | |||
105 | case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): | ||
106 | break; | ||
107 | } | ||
108 | /* rm = xm * ym, re = xe+ye basicly */ | ||
109 | assert(xm & DP_HIDDEN_BIT); | ||
110 | assert(ym & DP_HIDDEN_BIT); | ||
111 | { | ||
112 | int re = xe + ye; | ||
113 | int rs = xs ^ ys; | ||
114 | u64 rm; | ||
115 | |||
116 | /* shunt to top of word */ | ||
117 | xm <<= 64 - (DP_MBITS + 1); | ||
118 | ym <<= 64 - (DP_MBITS + 1); | ||
119 | |||
120 | /* multiply 32bits xm,ym to give high 32bits rm with stickness | ||
121 | */ | ||
122 | |||
123 | /* 32 * 32 => 64 */ | ||
124 | #define DPXMULT(x,y) ((u64)(x) * (u64)y) | ||
125 | |||
126 | { | ||
127 | unsigned lxm = xm; | ||
128 | unsigned hxm = xm >> 32; | ||
129 | unsigned lym = ym; | ||
130 | unsigned hym = ym >> 32; | ||
131 | u64 lrm; | ||
132 | u64 hrm; | ||
133 | |||
134 | lrm = DPXMULT(lxm, lym); | ||
135 | hrm = DPXMULT(hxm, hym); | ||
136 | |||
137 | { | ||
138 | u64 t = DPXMULT(lxm, hym); | ||
139 | { | ||
140 | u64 at = | ||
141 | lrm + (t << 32); | ||
142 | hrm += at < lrm; | ||
143 | lrm = at; | ||
144 | } | ||
145 | hrm = hrm + (t >> 32); | ||
146 | } | ||
147 | |||
148 | { | ||
149 | u64 t = DPXMULT(hxm, lym); | ||
150 | { | ||
151 | u64 at = | ||
152 | lrm + (t << 32); | ||
153 | hrm += at < lrm; | ||
154 | lrm = at; | ||
155 | } | ||
156 | hrm = hrm + (t >> 32); | ||
157 | } | ||
158 | rm = hrm | (lrm != 0); | ||
159 | } | ||
160 | |||
161 | /* | ||
162 | * sticky shift down to normal rounding precision | ||
163 | */ | ||
164 | if ((s64) rm < 0) { | ||
165 | rm = | ||
166 | (rm >> (64 - (DP_MBITS + 1 + 3))) | | ||
167 | ((rm << (DP_MBITS + 1 + 3)) != 0); | ||
168 | re++; | ||
169 | } else { | ||
170 | rm = | ||
171 | (rm >> (64 - (DP_MBITS + 1 + 3 + 1))) | | ||
172 | ((rm << (DP_MBITS + 1 + 3 + 1)) != 0); | ||
173 | } | ||
174 | assert(rm & (DP_HIDDEN_BIT << 3)); | ||
175 | DPNORMRET2(rs, re, rm, "mul", x, y); | ||
176 | } | ||
177 | } | ||