arch/mips/math-emu/sp_maddf.c

   1 /*
   2  * IEEE754 floating point arithmetic
   3  * single precision: MADDF.f (Fused Multiply Add)
   4  * MADDF.fmt: FPR[fd] = FPR[fd] + (FPR[fs] x FPR[ft])
   5  *
   6  * MIPS floating point support
   7  * Copyright (C) 2015 Imagination Technologies, Ltd.
   8  * Author: Markos Chandras <markos.chandras@imgtec.com>
   9  *
  10  *  This program is free software; you can distribute it and/or modify it
  11  *  under the terms of the GNU General Public License as published by the
  12  *  Free Software Foundation; version 2 of the License.
  13  */
  14
  15 #include "ieee754sp.h"
  16
  17
  18 static union ieee754sp _sp_maddf(union ieee754sp z, union ieee754sp x,
  19                                  union ieee754sp y, enum maddf_flags flags)
  20 {
  21         int re;
  22         int rs;
  23         unsigned int rm;
  24         u64 rm64;
  25         u64 zm64;
  26         int s;
  27
  28         COMPXSP;
  29         COMPYSP;
  30         COMPZSP;
  31
  32         EXPLODEXSP;
  33         EXPLODEYSP;
  34         EXPLODEZSP;
  35
  36         FLUSHXSP;
  37         FLUSHYSP;
  38         FLUSHZSP;
  39
  40         ieee754_clearcx();
  41
  42         /*
  43          * Handle the cases when at least one of x, y or z is a NaN.
  44          * Order of precedence is sNaN, qNaN and z, x, y.
  45          */
  46         if (zc == IEEE754_CLASS_SNAN)
  47                 return ieee754sp_nanxcpt(z);
  48         if (xc == IEEE754_CLASS_SNAN)
  49                 return ieee754sp_nanxcpt(x);
  50         if (yc == IEEE754_CLASS_SNAN)
  51                 return ieee754sp_nanxcpt(y);
  52         if (zc == IEEE754_CLASS_QNAN)
  53                 return z;
  54         if (xc == IEEE754_CLASS_QNAN)
  55                 return x;
  56         if (yc == IEEE754_CLASS_QNAN)
  57                 return y;
  58
  59         if (zc == IEEE754_CLASS_DNORM)
  60                 SPDNORMZ;
  61         /* ZERO z cases are handled separately below */
  62
  63         switch (CLPAIR(xc, yc)) {
  64
  65
  66         /*
  67          * Infinity handling
  68          */
  69         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
  70         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
  71                 ieee754_setcx(IEEE754_INVALID_OPERATION);
  72                 return ieee754sp_indef();
  73
  74         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
  75         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
  76         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
  77         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
  78         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
  79                 if ((zc == IEEE754_CLASS_INF) &&
  80                     ((!(flags & MADDF_NEGATE_PRODUCT) && (zs != (xs ^ ys))) ||
  81                      ((flags & MADDF_NEGATE_PRODUCT) && (zs == (xs ^ ys))))) {
  82                         /*
  83                          * Cases of addition of infinities with opposite signs
  84                          * or subtraction of infinities with same signs.
  85                          */
  86                         ieee754_setcx(IEEE754_INVALID_OPERATION);
  87                         return ieee754sp_indef();
  88                 }
  89                 /*
  90                  * z is here either not an infinity, or an infinity having the
  91                  * same sign as product (x*y) (in case of MADDF.D instruction)
  92                  * or product -(x*y) (in MSUBF.D case). The result must be an
  93                  * infinity, and its sign is determined only by the value of
  94                  * (flags & MADDF_NEGATE_PRODUCT) and the signs of x and y.
  95                  */
  96                 if (flags & MADDF_NEGATE_PRODUCT)
  97                         return ieee754sp_inf(1 ^ (xs ^ ys));
  98                 else
  99                         return ieee754sp_inf(xs ^ ys);
 100
 101         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
 102         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
 103         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
 104         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
 105         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
 106                 if (zc == IEEE754_CLASS_INF)
 107                         return ieee754sp_inf(zs);
 108                 if (zc == IEEE754_CLASS_ZERO) {
 109                         /* Handle cases +0 + (-0) and similar ones. */
 110                         if ((!(flags & MADDF_NEGATE_PRODUCT)
 111                                         && (zs == (xs ^ ys))) ||
 112                             ((flags & MADDF_NEGATE_PRODUCT)
 113                                         && (zs != (xs ^ ys))))
 114                                 /*
 115                                  * Cases of addition of zeros of equal signs
 116                                  * or subtraction of zeroes of opposite signs.
 117                                  * The sign of the resulting zero is in any
 118                                  * such case determined only by the sign of z.
 119                                  */
 120                                 return z;
 121
 122                         return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
 123                 }
 124                 /* x*y is here 0, and z is not 0, so just return z */
 125                 return z;
 126
 127         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
 128                 SPDNORMX;
 129                 /* fall through */
 130
 131         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
 132                 if (zc == IEEE754_CLASS_INF)
 133                         return ieee754sp_inf(zs);
 134                 SPDNORMY;
 135                 break;
 136
 137         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
 138                 if (zc == IEEE754_CLASS_INF)
 139                         return ieee754sp_inf(zs);
 140                 SPDNORMX;
 141                 break;
 142
 143         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
 144                 if (zc == IEEE754_CLASS_INF)
 145                         return ieee754sp_inf(zs);
 146                 /* continue to real computations */
 147         }
 148
 149         /* Finally get to do some computation */
 150
 151         /*
 152          * Do the multiplication bit first
 153          *
 154          * rm = xm * ym, re = xe + ye basically
 155          *
 156          * At this point xm and ym should have been normalized.
 157          */
 158
 159         /* rm = xm * ym, re = xe+ye basically */
 160         assert(xm & SP_HIDDEN_BIT);
 161         assert(ym & SP_HIDDEN_BIT);
 162
 163         re = xe + ye;
 164         rs = xs ^ ys;
 165         if (flags & MADDF_NEGATE_PRODUCT)
 166                 rs ^= 1;
 167
 168         /* Multiple 24 bit xm and ym to give 48 bit results */
 169         rm64 = (uint64_t)xm * ym;
 170
 171         /* Shunt to top of word */
 172         rm64 = rm64 << 16;
 173
 174         /* Put explicit bit at bit 62 if necessary */
 175         if ((int64_t) rm64 < 0) {
 176                 rm64 = rm64 >> 1;
 177                 re++;
 178         }
 179
 180         assert(rm64 & (1 << 62));
 181
 182         if (zc == IEEE754_CLASS_ZERO) {
 183                 /*
 184                  * Move explicit bit from bit 62 to bit 26 since the
 185                  * ieee754sp_format code expects the mantissa to be
 186                  * 27 bits wide (24 + 3 rounding bits).
 187                  */
 188                 rm = XSPSRS64(rm64, (62 - 26));
 189                 return ieee754sp_format(rs, re, rm);
 190         }
 191
 192         /* Move explicit bit from bit 23 to bit 62 */
 193         zm64 = (uint64_t)zm << (62 - 23);
 194         assert(zm64 & (1 << 62));
 195
 196         /* Make the exponents the same */
 197         if (ze > re) {
 198                 /*
 199                  * Have to shift r fraction right to align.
 200                  */
 201                 s = ze - re;
 202                 rm64 = XSPSRS64(rm64, s);
 203                 re += s;
 204         } else if (re > ze) {
 205                 /*
 206                  * Have to shift z fraction right to align.
 207                  */
 208                 s = re - ze;
 209                 zm64 = XSPSRS64(zm64, s);
 210                 ze += s;
 211         }
 212         assert(ze == re);
 213         assert(ze <= SP_EMAX);
 214
 215         /* Do the addition */
 216         if (zs == rs) {
 217                 /*
 218                  * Generate 64 bit result by adding two 63 bit numbers
 219                  * leaving result in zm64, zs and ze.
 220                  */
 221                 zm64 = zm64 + rm64;
 222                 if ((int64_t)zm64 < 0) {        /* carry out */
 223                         zm64 = XSPSRS1(zm64);
 224                         ze++;
 225                 }
 226         } else {
 227                 if (zm64 >= rm64) {
 228                         zm64 = zm64 - rm64;
 229                 } else {
 230                         zm64 = rm64 - zm64;
 231                         zs = rs;
 232                 }
 233                 if (zm64 == 0)
 234                         return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
 235
 236                 /*
 237                  * Put explicit bit at bit 62 if necessary.
 238                  */
 239                 while ((zm64 >> 62) == 0) {
 240                         zm64 <<= 1;
 241                         ze--;
 242                 }
 243         }
 244
 245         /*
 246          * Move explicit bit from bit 62 to bit 26 since the
 247          * ieee754sp_format code expects the mantissa to be
 248          * 27 bits wide (24 + 3 rounding bits).
 249          */
 250         zm = XSPSRS64(zm64, (62 - 26));
 251
 252         return ieee754sp_format(zs, ze, zm);
 253 }
 254
 255 union ieee754sp ieee754sp_maddf(union ieee754sp z, union ieee754sp x,
 256                                 union ieee754sp y)
 257 {
 258         return _sp_maddf(z, x, y, 0);
 259 }
 260
 261 union ieee754sp ieee754sp_msubf(union ieee754sp z, union ieee754sp x,
 262                                 union ieee754sp y)
 263 {
 264         return _sp_maddf(z, x, y, MADDF_NEGATE_PRODUCT);
 265 }