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 /* fall through */ 95 96 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): 97 DPDNORMY; 98 break; 99 100 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): 101 DPDNORMX; 102 break; 103 104 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): 105 break; 106 } 107 assert(xm & DP_HIDDEN_BIT); 108 assert(ym & DP_HIDDEN_BIT); 109 110 /* provide rounding space */ 111 xm <<= 3; 112 ym <<= 3; 113 114 /* now the dirty work */ 115 116 rm = 0; 117 re = xe - ye; 118 119 for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) { 120 if (xm >= ym) { 121 xm -= ym; 122 rm |= bm; 123 if (xm == 0) 124 break; 125 } 126 xm <<= 1; 127 } 128 129 rm <<= 1; 130 if (xm) 131 rm |= 1; /* have remainder, set sticky */ 132 133 assert(rm); 134 135 /* 136 * Normalise rm to rounding precision ? 137 */ 138 while ((rm >> (DP_FBITS + 3)) == 0) { 139 rm <<= 1; 140 re--; 141 } 142 143 return ieee754dp_format(xs == ys ? 0 : 1, re, rm); 144 } 145