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