1 // SPDX-License-Identifier: GPL-2.0-only
2 /*
3  * IEEE754 floating point arithmetic
4  * double precision: MADDF.f (Fused Multiply Add)
5  * MADDF.fmt: FPR[fd] = FPR[fd] + (FPR[fs] x FPR[ft])
6  *
7  * MIPS floating point support
8  * Copyright (C) 2015 Imagination Technologies, Ltd.
9  * Author: Markos Chandras <markos.chandras@imgtec.com>
10  */
11 
12 #include "ieee754dp.h"
13 
14 
15 /* 128 bits shift right logical with rounding. */
srl128(u64 * hptr,u64 * lptr,int count)16 static void srl128(u64 *hptr, u64 *lptr, int count)
17 {
18 	u64 low;
19 
20 	if (count >= 128) {
21 		*lptr = *hptr != 0 || *lptr != 0;
22 		*hptr = 0;
23 	} else if (count >= 64) {
24 		if (count == 64) {
25 			*lptr = *hptr | (*lptr != 0);
26 		} else {
27 			low = *lptr;
28 			*lptr = *hptr >> (count - 64);
29 			*lptr |= (*hptr << (128 - count)) != 0 || low != 0;
30 		}
31 		*hptr = 0;
32 	} else {
33 		low = *lptr;
34 		*lptr = low >> count | *hptr << (64 - count);
35 		*lptr |= (low << (64 - count)) != 0;
36 		*hptr = *hptr >> count;
37 	}
38 }
39 
_dp_maddf(union ieee754dp z,union ieee754dp x,union ieee754dp y,enum maddf_flags flags)40 static union ieee754dp _dp_maddf(union ieee754dp z, union ieee754dp x,
41 				 union ieee754dp y, enum maddf_flags flags)
42 {
43 	int re;
44 	int rs;
45 	unsigned int lxm;
46 	unsigned int hxm;
47 	unsigned int lym;
48 	unsigned int hym;
49 	u64 lrm;
50 	u64 hrm;
51 	u64 lzm;
52 	u64 hzm;
53 	u64 t;
54 	u64 at;
55 	int s;
56 
57 	COMPXDP;
58 	COMPYDP;
59 	COMPZDP;
60 
61 	EXPLODEXDP;
62 	EXPLODEYDP;
63 	EXPLODEZDP;
64 
65 	FLUSHXDP;
66 	FLUSHYDP;
67 	FLUSHZDP;
68 
69 	ieee754_clearcx();
70 
71 	rs = xs ^ ys;
72 	if (flags & MADDF_NEGATE_PRODUCT)
73 		rs ^= 1;
74 	if (flags & MADDF_NEGATE_ADDITION)
75 		zs ^= 1;
76 
77 	/*
78 	 * Handle the cases when at least one of x, y or z is a NaN.
79 	 * Order of precedence is sNaN, qNaN and z, x, y.
80 	 */
81 	if (zc == IEEE754_CLASS_SNAN)
82 		return ieee754dp_nanxcpt(z);
83 	if (xc == IEEE754_CLASS_SNAN)
84 		return ieee754dp_nanxcpt(x);
85 	if (yc == IEEE754_CLASS_SNAN)
86 		return ieee754dp_nanxcpt(y);
87 	if (zc == IEEE754_CLASS_QNAN)
88 		return z;
89 	if (xc == IEEE754_CLASS_QNAN)
90 		return x;
91 	if (yc == IEEE754_CLASS_QNAN)
92 		return y;
93 
94 	if (zc == IEEE754_CLASS_DNORM)
95 		DPDNORMZ;
96 	/* ZERO z cases are handled separately below */
97 
98 	switch (CLPAIR(xc, yc)) {
99 
100 	/*
101 	 * Infinity handling
102 	 */
103 	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
104 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
105 		ieee754_setcx(IEEE754_INVALID_OPERATION);
106 		return ieee754dp_indef();
107 
108 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
109 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
110 	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
111 	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
112 	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
113 		if ((zc == IEEE754_CLASS_INF) && (zs != rs)) {
114 			/*
115 			 * Cases of addition of infinities with opposite signs
116 			 * or subtraction of infinities with same signs.
117 			 */
118 			ieee754_setcx(IEEE754_INVALID_OPERATION);
119 			return ieee754dp_indef();
120 		}
121 		/*
122 		 * z is here either not an infinity, or an infinity having the
123 		 * same sign as product (x*y). The result must be an infinity,
124 		 * and its sign is determined only by the sign of product (x*y).
125 		 */
126 		return ieee754dp_inf(rs);
127 
128 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
129 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
130 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
131 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
132 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
133 		if (zc == IEEE754_CLASS_INF)
134 			return ieee754dp_inf(zs);
135 		if (zc == IEEE754_CLASS_ZERO) {
136 			/* Handle cases +0 + (-0) and similar ones. */
137 			if (zs == rs)
138 				/*
139 				 * Cases of addition of zeros of equal signs
140 				 * or subtraction of zeroes of opposite signs.
141 				 * The sign of the resulting zero is in any
142 				 * such case determined only by the sign of z.
143 				 */
144 				return z;
145 
146 			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
147 		}
148 		/* x*y is here 0, and z is not 0, so just return z */
149 		return z;
150 
151 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
152 		DPDNORMX;
153 		fallthrough;
154 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
155 		if (zc == IEEE754_CLASS_INF)
156 			return ieee754dp_inf(zs);
157 		DPDNORMY;
158 		break;
159 
160 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
161 		if (zc == IEEE754_CLASS_INF)
162 			return ieee754dp_inf(zs);
163 		DPDNORMX;
164 		break;
165 
166 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
167 		if (zc == IEEE754_CLASS_INF)
168 			return ieee754dp_inf(zs);
169 		/* continue to real computations */
170 	}
171 
172 	/* Finally get to do some computation */
173 
174 	/*
175 	 * Do the multiplication bit first
176 	 *
177 	 * rm = xm * ym, re = xe + ye basically
178 	 *
179 	 * At this point xm and ym should have been normalized.
180 	 */
181 	assert(xm & DP_HIDDEN_BIT);
182 	assert(ym & DP_HIDDEN_BIT);
183 
184 	re = xe + ye;
185 
186 	/* shunt to top of word */
187 	xm <<= 64 - (DP_FBITS + 1);
188 	ym <<= 64 - (DP_FBITS + 1);
189 
190 	/*
191 	 * Multiply 64 bits xm and ym to give 128 bits result in hrm:lrm.
192 	 */
193 
194 	lxm = xm;
195 	hxm = xm >> 32;
196 	lym = ym;
197 	hym = ym >> 32;
198 
199 	lrm = DPXMULT(lxm, lym);
200 	hrm = DPXMULT(hxm, hym);
201 
202 	t = DPXMULT(lxm, hym);
203 
204 	at = lrm + (t << 32);
205 	hrm += at < lrm;
206 	lrm = at;
207 
208 	hrm = hrm + (t >> 32);
209 
210 	t = DPXMULT(hxm, lym);
211 
212 	at = lrm + (t << 32);
213 	hrm += at < lrm;
214 	lrm = at;
215 
216 	hrm = hrm + (t >> 32);
217 
218 	/* Put explicit bit at bit 126 if necessary */
219 	if ((int64_t)hrm < 0) {
220 		lrm = (hrm << 63) | (lrm >> 1);
221 		hrm = hrm >> 1;
222 		re++;
223 	}
224 
225 	assert(hrm & (1 << 62));
226 
227 	if (zc == IEEE754_CLASS_ZERO) {
228 		/*
229 		 * Move explicit bit from bit 126 to bit 55 since the
230 		 * ieee754dp_format code expects the mantissa to be
231 		 * 56 bits wide (53 + 3 rounding bits).
232 		 */
233 		srl128(&hrm, &lrm, (126 - 55));
234 		return ieee754dp_format(rs, re, lrm);
235 	}
236 
237 	/* Move explicit bit from bit 52 to bit 126 */
238 	lzm = 0;
239 	hzm = zm << 10;
240 	assert(hzm & (1 << 62));
241 
242 	/* Make the exponents the same */
243 	if (ze > re) {
244 		/*
245 		 * Have to shift y fraction right to align.
246 		 */
247 		s = ze - re;
248 		srl128(&hrm, &lrm, s);
249 		re += s;
250 	} else if (re > ze) {
251 		/*
252 		 * Have to shift x fraction right to align.
253 		 */
254 		s = re - ze;
255 		srl128(&hzm, &lzm, s);
256 		ze += s;
257 	}
258 	assert(ze == re);
259 	assert(ze <= DP_EMAX);
260 
261 	/* Do the addition */
262 	if (zs == rs) {
263 		/*
264 		 * Generate 128 bit result by adding two 127 bit numbers
265 		 * leaving result in hzm:lzm, zs and ze.
266 		 */
267 		hzm = hzm + hrm + (lzm > (lzm + lrm));
268 		lzm = lzm + lrm;
269 		if ((int64_t)hzm < 0) {        /* carry out */
270 			srl128(&hzm, &lzm, 1);
271 			ze++;
272 		}
273 	} else {
274 		if (hzm > hrm || (hzm == hrm && lzm >= lrm)) {
275 			hzm = hzm - hrm - (lzm < lrm);
276 			lzm = lzm - lrm;
277 		} else {
278 			hzm = hrm - hzm - (lrm < lzm);
279 			lzm = lrm - lzm;
280 			zs = rs;
281 		}
282 		if (lzm == 0 && hzm == 0)
283 			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
284 
285 		/*
286 		 * Put explicit bit at bit 126 if necessary.
287 		 */
288 		if (hzm == 0) {
289 			/* left shift by 63 or 64 bits */
290 			if ((int64_t)lzm < 0) {
291 				/* MSB of lzm is the explicit bit */
292 				hzm = lzm >> 1;
293 				lzm = lzm << 63;
294 				ze -= 63;
295 			} else {
296 				hzm = lzm;
297 				lzm = 0;
298 				ze -= 64;
299 			}
300 		}
301 
302 		t = 0;
303 		while ((hzm >> (62 - t)) == 0)
304 			t++;
305 
306 		assert(t <= 62);
307 		if (t) {
308 			hzm = hzm << t | lzm >> (64 - t);
309 			lzm = lzm << t;
310 			ze -= t;
311 		}
312 	}
313 
314 	/*
315 	 * Move explicit bit from bit 126 to bit 55 since the
316 	 * ieee754dp_format code expects the mantissa to be
317 	 * 56 bits wide (53 + 3 rounding bits).
318 	 */
319 	srl128(&hzm, &lzm, (126 - 55));
320 
321 	return ieee754dp_format(zs, ze, lzm);
322 }
323 
ieee754dp_maddf(union ieee754dp z,union ieee754dp x,union ieee754dp y)324 union ieee754dp ieee754dp_maddf(union ieee754dp z, union ieee754dp x,
325 				union ieee754dp y)
326 {
327 	return _dp_maddf(z, x, y, 0);
328 }
329 
ieee754dp_msubf(union ieee754dp z,union ieee754dp x,union ieee754dp y)330 union ieee754dp ieee754dp_msubf(union ieee754dp z, union ieee754dp x,
331 				union ieee754dp y)
332 {
333 	return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT);
334 }
335 
ieee754dp_madd(union ieee754dp z,union ieee754dp x,union ieee754dp y)336 union ieee754dp ieee754dp_madd(union ieee754dp z, union ieee754dp x,
337 				union ieee754dp y)
338 {
339 	return _dp_maddf(z, x, y, 0);
340 }
341 
ieee754dp_msub(union ieee754dp z,union ieee754dp x,union ieee754dp y)342 union ieee754dp ieee754dp_msub(union ieee754dp z, union ieee754dp x,
343 				union ieee754dp y)
344 {
345 	return _dp_maddf(z, x, y, MADDF_NEGATE_ADDITION);
346 }
347 
ieee754dp_nmadd(union ieee754dp z,union ieee754dp x,union ieee754dp y)348 union ieee754dp ieee754dp_nmadd(union ieee754dp z, union ieee754dp x,
349 				union ieee754dp y)
350 {
351 	return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT|MADDF_NEGATE_ADDITION);
352 }
353 
ieee754dp_nmsub(union ieee754dp z,union ieee754dp x,union ieee754dp y)354 union ieee754dp ieee754dp_nmsub(union ieee754dp z, union ieee754dp x,
355 				union ieee754dp y)
356 {
357 	return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT);
358 }
359