https://github.com/torvalds/linux
Tip revision: 7c53f6b671f4aba70ff15e1b05148b10d58c2837 authored by Linus Torvalds on 10 January 2021, 22:34:50 UTC
Linux 5.11-rc3
Linux 5.11-rc3
Tip revision: 7c53f6b
dp_div.c
// SPDX-License-Identifier: GPL-2.0-only
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
*/
#include "ieee754dp.h"
union ieee754dp ieee754dp_div(union ieee754dp x, union ieee754dp y)
{
u64 rm;
int re;
u64 bm;
COMPXDP;
COMPYDP;
EXPLODEXDP;
EXPLODEYDP;
ieee754_clearcx();
FLUSHXDP;
FLUSHYDP;
switch (CLPAIR(xc, yc)) {
case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
return ieee754dp_nanxcpt(y);
case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
return ieee754dp_nanxcpt(x);
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
return y;
case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
return x;
/*
* Infinity handling
*/
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
ieee754_setcx(IEEE754_INVALID_OPERATION);
return ieee754dp_indef();
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
return ieee754dp_zero(xs ^ ys);
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
return ieee754dp_inf(xs ^ ys);
/*
* Zero handling
*/
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
ieee754_setcx(IEEE754_INVALID_OPERATION);
return ieee754dp_indef();
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
ieee754_setcx(IEEE754_ZERO_DIVIDE);
return ieee754dp_inf(xs ^ ys);
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
return ieee754dp_zero(xs == ys ? 0 : 1);
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
DPDNORMX;
fallthrough;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
DPDNORMY;
break;
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
DPDNORMX;
break;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
break;
}
assert(xm & DP_HIDDEN_BIT);
assert(ym & DP_HIDDEN_BIT);
/* provide rounding space */
xm <<= 3;
ym <<= 3;
/* now the dirty work */
rm = 0;
re = xe - ye;
for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) {
if (xm >= ym) {
xm -= ym;
rm |= bm;
if (xm == 0)
break;
}
xm <<= 1;
}
rm <<= 1;
if (xm)
rm |= 1; /* have remainder, set sticky */
assert(rm);
/*
* Normalise rm to rounding precision ?
*/
while ((rm >> (DP_FBITS + 3)) == 0) {
rm <<= 1;
re--;
}
return ieee754dp_format(xs == ys ? 0 : 1, re, rm);
}