Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
torvalds
GitHub Repository: torvalds/linux
 Path: blob/master/arch/mips/math-emu/dp_mul.c
26424 views
1
// SPDX-License-Identifier: GPL-2.0-only

2
/* IEEE754 floating point arithmetic

3
 * double precision: common utilities

4
 */

5
/*

6
 * MIPS floating point support

7
 * Copyright (C) 1994-2000 Algorithmics Ltd.

8
 */

9


10
#include "ieee754dp.h"

11


12
union ieee754dp ieee754dp_mul(union ieee754dp x, union ieee754dp y)

13
{

14
	int re;

15
	int rs;

16
	u64 rm;

17
	unsigned int lxm;

18
	unsigned int hxm;

19
	unsigned int lym;

20
	unsigned int hym;

21
	u64 lrm;

22
	u64 hrm;

23
	u64 t;

24
	u64 at;

25


26
	COMPXDP;

27
	COMPYDP;

28


29
	EXPLODEXDP;

30
	EXPLODEYDP;

31


32
	ieee754_clearcx();

33


34
	FLUSHXDP;

35
	FLUSHYDP;

36


37
	switch (CLPAIR(xc, yc)) {

38
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):

39
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):

40
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):

41
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):

42
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):

43
		return ieee754dp_nanxcpt(y);

44


45
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):

46
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):

47
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):

48
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):

49
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):

50
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):

51
		return ieee754dp_nanxcpt(x);

52


53
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):

54
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):

55
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):

56
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):

57
		return y;

58


59
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):

60
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):

61
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):

62
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):

63
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):

64
		return x;

65


66


67
	/*

68
	 * Infinity handling

69
	 */

70
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):

71
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):

72
		ieee754_setcx(IEEE754_INVALID_OPERATION);

73
		return ieee754dp_indef();

74


75
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):

76
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):

77
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):

78
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):

79
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):

80
		return ieee754dp_inf(xs ^ ys);

81


82
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):

83
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):

84
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):

85
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):

86
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):

87
		return ieee754dp_zero(xs ^ ys);

88


89


90
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):

91
		DPDNORMX;

92
		fallthrough;

93
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):

94
		DPDNORMY;

95
		break;

96


97
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):

98
		DPDNORMX;

99
		break;

100


101
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):

102
		break;

103
	}

104
	/* rm = xm * ym, re = xe+ye basically */

105
	assert(xm & DP_HIDDEN_BIT);

106
	assert(ym & DP_HIDDEN_BIT);

107


108
	re = xe + ye;

109
	rs = xs ^ ys;

110


111
	/* shunt to top of word */

112
	xm <<= 64 - (DP_FBITS + 1);

113
	ym <<= 64 - (DP_FBITS + 1);

114


115
	/*

116
	 * Multiply 64 bits xm, ym to give high 64 bits rm with stickness.

117
	 */

118


119
	lxm = xm;

120
	hxm = xm >> 32;

121
	lym = ym;

122
	hym = ym >> 32;

123


124
	lrm = DPXMULT(lxm, lym);

125
	hrm = DPXMULT(hxm, hym);

126


127
	t = DPXMULT(lxm, hym);

128


129
	at = lrm + (t << 32);

130
	hrm += at < lrm;

131
	lrm = at;

132


133
	hrm = hrm + (t >> 32);

134


135
	t = DPXMULT(hxm, lym);

136


137
	at = lrm + (t << 32);

138
	hrm += at < lrm;

139
	lrm = at;

140


141
	hrm = hrm + (t >> 32);

142


143
	rm = hrm | (lrm != 0);

144


145
	/*

146
	 * Sticky shift down to normal rounding precision.

147
	 */

148
	if ((s64) rm < 0) {

149
		rm = (rm >> (64 - (DP_FBITS + 1 + 3))) |

150
		     ((rm << (DP_FBITS + 1 + 3)) != 0);

151
		re++;

152
	} else {

153
		rm = (rm >> (64 - (DP_FBITS + 1 + 3 + 1))) |

154
		     ((rm << (DP_FBITS + 1 + 3 + 1)) != 0);

155
	}

156
	assert(rm & (DP_HIDDEN_BIT << 3));

157


158
	return ieee754dp_format(rs, re, rm);

159
}

160


161