mips/math-emu/dp_div.c

0001 // SPDX-License-Identifier: GPL-2.0-only
0002 /* IEEE754 floating point arithmetic
0003  * double precision: common utilities
0004  */
0005 /*
0006  * MIPS floating point support
0007  * Copyright (C) 1994-2000 Algorithmics Ltd.
0008  */
0009
0010 #include "ieee754dp.h"
0011
0012 union ieee754dp ieee754dp_div(union ieee754dp x, union ieee754dp y)
0013 {
0014     u64 rm;
0015     int re;
0016     u64 bm;
0017
0018     COMPXDP;
0019     COMPYDP;
0020
0021     EXPLODEXDP;
0022     EXPLODEYDP;
0023
0024     ieee754_clearcx();
0025
0026     FLUSHXDP;
0027     FLUSHYDP;
0028
0029     switch (CLPAIR(xc, yc)) {
0030     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
0031     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
0032     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
0033     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
0034     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
0035         return ieee754dp_nanxcpt(y);
0036
0037     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
0038     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
0039     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
0040     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
0041     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
0042     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
0043         return ieee754dp_nanxcpt(x);
0044
0045     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
0046     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
0047     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
0048     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
0049         return y;
0050
0051     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
0052     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
0053     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
0054     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
0055     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
0056         return x;
0057
0058
0059     /*
0060      * Infinity handling
0061      */
0062     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
0063         ieee754_setcx(IEEE754_INVALID_OPERATION);
0064         return ieee754dp_indef();
0065
0066     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
0067     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
0068     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
0069         return ieee754dp_zero(xs ^ ys);
0070
0071     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
0072     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
0073     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
0074         return ieee754dp_inf(xs ^ ys);
0075
0076     /*
0077      * Zero handling
0078      */
0079     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
0080         ieee754_setcx(IEEE754_INVALID_OPERATION);
0081         return ieee754dp_indef();
0082
0083     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
0084     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
0085         ieee754_setcx(IEEE754_ZERO_DIVIDE);
0086         return ieee754dp_inf(xs ^ ys);
0087
0088     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
0089     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
0090         return ieee754dp_zero(xs == ys ? 0 : 1);
0091
0092     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
0093         DPDNORMX;
0094         fallthrough;
0095     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
0096         DPDNORMY;
0097         break;
0098
0099     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
0100         DPDNORMX;
0101         break;
0102
0103     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
0104         break;
0105     }
0106     assert(xm & DP_HIDDEN_BIT);
0107     assert(ym & DP_HIDDEN_BIT);
0108
0109     /* provide rounding space */
0110     xm <<= 3;
0111     ym <<= 3;
0112
0113     /* now the dirty work */
0114
0115     rm = 0;
0116     re = xe - ye;
0117
0118     for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) {
0119         if (xm >= ym) {
0120             xm -= ym;
0121             rm |= bm;
0122             if (xm == 0)
0123                 break;
0124         }
0125         xm <<= 1;
0126     }
0127
0128     rm <<= 1;
0129     if (xm)
0130         rm |= 1;    /* have remainder, set sticky */
0131
0132     assert(rm);
0133
0134     /*
0135      * Normalise rm to rounding precision ?
0136      */
0137     while ((rm >> (DP_FBITS + 3)) == 0) {
0138         rm <<= 1;
0139         re--;
0140     }
0141
0142     return ieee754dp_format(xs == ys ? 0 : 1, re, rm);
0143 }