Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
torvalds
GitHub Repository: torvalds/linux
 Path: blob/master/arch/mips/math-emu/dp_div.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_div(union ieee754dp x, union ieee754dp y)

13
{

14
	u64 rm;

15
	int re;

16
	u64 bm;

17


18
	COMPXDP;

19
	COMPYDP;

20


21
	EXPLODEXDP;

22
	EXPLODEYDP;

23


24
	ieee754_clearcx();

25


26
	FLUSHXDP;

27
	FLUSHYDP;

28


29
	switch (CLPAIR(xc, yc)) {

30
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):

31
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):

32
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):

33
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):

34
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):

35
		return ieee754dp_nanxcpt(y);

36


37
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):

38
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):

39
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):

40
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):

41
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):

42
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):

43
		return ieee754dp_nanxcpt(x);

44


45
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):

46
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):

47
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):

48
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):

49
		return y;

50


51
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):

52
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):

53
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):

54
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):

55
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):

56
		return x;

57


58


59
	/*

60
	 * Infinity handling

61
	 */

62
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):

63
		ieee754_setcx(IEEE754_INVALID_OPERATION);

64
		return ieee754dp_indef();

65


66
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):

67
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):

68
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):

69
		return ieee754dp_zero(xs ^ ys);

70


71
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):

72
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):

73
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):

74
		return ieee754dp_inf(xs ^ ys);

75


76
	/*

77
	 * Zero handling

78
	 */

79
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):

80
		ieee754_setcx(IEEE754_INVALID_OPERATION);

81
		return ieee754dp_indef();

82


83
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):

84
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):

85
		ieee754_setcx(IEEE754_ZERO_DIVIDE);

86
		return ieee754dp_inf(xs ^ ys);

87


88
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):

89
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):

90
		return ieee754dp_zero(xs == ys ? 0 : 1);

91


92
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):

93
		DPDNORMX;

94
		fallthrough;

95
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):

96
		DPDNORMY;

97
		break;

98


99
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):

100
		DPDNORMX;

101
		break;

102


103
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):

104
		break;

105
	}

106
	assert(xm & DP_HIDDEN_BIT);

107
	assert(ym & DP_HIDDEN_BIT);

108


109
	/* provide rounding space */

110
	xm <<= 3;

111
	ym <<= 3;

112


113
	/* now the dirty work */

114


115
	rm = 0;

116
	re = xe - ye;

117


118
	for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) {

119
		if (xm >= ym) {

120
			xm -= ym;

121
			rm |= bm;

122
			if (xm == 0)

123
				break;

124
		}

125
		xm <<= 1;

126
	}

127


128
	rm <<= 1;

129
	if (xm)

130
		rm |= 1;	/* have remainder, set sticky */

131


132
	assert(rm);

133


134
	/*

135
	 * Normalise rm to rounding precision ?

136
	 */

137
	while ((rm >> (DP_FBITS + 3)) == 0) {

138
		rm <<= 1;

139
		re--;

140
	}

141


142
	return ieee754dp_format(xs == ys ? 0 : 1, re, rm);

143
}

144


145