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 enum maddf_flags {
  18         maddf_negate_product    = 1 << 0,
  19 };
  20
  21 static union ieee754sp _sp_maddf(union ieee754sp z, union ieee754sp x,
  22                                  union ieee754sp y, enum maddf_flags flags)
  23 {
  24         int re;
  25         int rs;
  26         unsigned rm;
  27         unsigned short lxm;
  28         unsigned short hxm;
  29         unsigned short lym;
  30         unsigned short hym;
  31         unsigned lrm;
  32         unsigned hrm;
  33         unsigned t;
  34         unsigned at;
  35         int s;
  36
  37         COMPXSP;
  38         COMPYSP;
  39         COMPZSP;
  40
  41         EXPLODEXSP;
  42         EXPLODEYSP;
  43         EXPLODEZSP;
  44
  45         FLUSHXSP;
  46         FLUSHYSP;
  47         FLUSHZSP;
  48
  49         ieee754_clearcx();
  50
  51         switch (zc) {
  52         case IEEE754_CLASS_SNAN:
  53                 ieee754_setcx(IEEE754_INVALID_OPERATION);
  54                 return ieee754sp_nanxcpt(z);
  55         case IEEE754_CLASS_DNORM:
  56                 SPDNORMZ;
  57         /* QNAN is handled separately below */
  58         }
  59
  60         switch (CLPAIR(xc, yc)) {
  61         case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
  62         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
  63         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
  64         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
  65         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
  66                 return ieee754sp_nanxcpt(y);
  67
  68         case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
  69         case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
  70         case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
  71         case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
  72         case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
  73         case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
  74                 return ieee754sp_nanxcpt(x);
  75
  76         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
  77         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
  78         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
  79         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
  80                 return y;
  81
  82         case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
  83         case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
  84         case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
  85         case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
  86         case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
  87                 return x;
  88
  89         /*
  90          * Infinity handling
  91          */
  92         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
  93         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
  94                 if (zc == IEEE754_CLASS_QNAN)
  95                         return z;
  96                 ieee754_setcx(IEEE754_INVALID_OPERATION);
  97                 return ieee754sp_indef();
  98
  99         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
 100         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
 101         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
 102         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
 103         case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
 104                 if (zc == IEEE754_CLASS_QNAN)
 105                         return z;
 106                 return ieee754sp_inf(xs ^ ys);
 107
 108         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
 109         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
 110         case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
 111         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
 112         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
 113                 if (zc == IEEE754_CLASS_INF)
 114                         return ieee754sp_inf(zs);
 115                 /* Multiplication is 0 so just return z */
 116                 return z;
 117
 118         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
 119                 SPDNORMX;
 120
 121         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
 122                 if (zc == IEEE754_CLASS_QNAN)
 123                         return z;
 124                 else if (zc == IEEE754_CLASS_INF)
 125                         return ieee754sp_inf(zs);
 126                 SPDNORMY;
 127                 break;
 128
 129         case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
 130                 if (zc == IEEE754_CLASS_QNAN)
 131                         return z;
 132                 else if (zc == IEEE754_CLASS_INF)
 133                         return ieee754sp_inf(zs);
 134                 SPDNORMX;
 135                 break;
 136
 137         case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
 138                 if (zc == IEEE754_CLASS_QNAN)
 139                         return z;
 140                 else if (zc == IEEE754_CLASS_INF)
 141                         return ieee754sp_inf(zs);
 142                 /* fall through to real computations */
 143         }
 144
 145         /* Finally get to do some computation */
 146
 147         /*
 148          * Do the multiplication bit first
 149          *
 150          * rm = xm * ym, re = xe + ye basically
 151          *
 152          * At this point xm and ym should have been normalized.
 153          */
 154
 155         /* rm = xm * ym, re = xe+ye basically */
 156         assert(xm & SP_HIDDEN_BIT);
 157         assert(ym & SP_HIDDEN_BIT);
 158
 159         re = xe + ye;
 160         rs = xs ^ ys;
 161         if (flags & maddf_negate_product)
 162                 rs ^= 1;
 163
 164         /* shunt to top of word */
 165         xm <<= 32 - (SP_FBITS + 1);
 166         ym <<= 32 - (SP_FBITS + 1);
 167
 168         /*
 169          * Multiply 32 bits xm, ym to give high 32 bits rm with stickness.
 170          */
 171         lxm = xm & 0xffff;
 172         hxm = xm >> 16;
 173         lym = ym & 0xffff;
 174         hym = ym >> 16;
 175
 176         lrm = lxm * lym;        /* 16 * 16 => 32 */
 177         hrm = hxm * hym;        /* 16 * 16 => 32 */
 178
 179         t = lxm * hym; /* 16 * 16 => 32 */
 180         at = lrm + (t << 16);
 181         hrm += at < lrm;
 182         lrm = at;
 183         hrm = hrm + (t >> 16);
 184
 185         t = hxm * lym; /* 16 * 16 => 32 */
 186         at = lrm + (t << 16);
 187         hrm += at < lrm;
 188         lrm = at;
 189         hrm = hrm + (t >> 16);
 190
 191         rm = hrm | (lrm != 0);
 192
 193         /*
 194          * Sticky shift down to normal rounding precision.
 195          */
 196         if ((int) rm < 0) {
 197                 rm = (rm >> (32 - (SP_FBITS + 1 + 3))) |
 198                     ((rm << (SP_FBITS + 1 + 3)) != 0);
 199                 re++;
 200         } else {
 201                 rm = (rm >> (32 - (SP_FBITS + 1 + 3 + 1))) |
 202                      ((rm << (SP_FBITS + 1 + 3 + 1)) != 0);
 203         }
 204         assert(rm & (SP_HIDDEN_BIT << 3));
 205
 206         /* And now the addition */
 207
 208         assert(zm & SP_HIDDEN_BIT);
 209
 210         /*
 211          * Provide guard,round and stick bit space.
 212          */
 213         zm <<= 3;
 214
 215         if (ze > re) {
 216                 /*
 217                  * Have to shift r fraction right to align.
 218                  */
 219                 s = ze - re;
 220                 rm = XSPSRS(rm, s);
 221                 re += s;
 222         } else if (re > ze) {
 223                 /*
 224                  * Have to shift z fraction right to align.
 225                  */
 226                 s = re - ze;
 227                 zm = XSPSRS(zm, s);
 228                 ze += s;
 229         }
 230         assert(ze == re);
 231         assert(ze <= SP_EMAX);
 232
 233         if (zs == rs) {
 234                 /*
 235                  * Generate 28 bit result of adding two 27 bit numbers
 236                  * leaving result in zm, zs and ze.
 237                  */
 238                 zm = zm + rm;
 239
 240                 if (zm >> (SP_FBITS + 1 + 3)) { /* carry out */
 241                         zm = XSPSRS1(zm);
 242                         ze++;
 243                 }
 244         } else {
 245                 if (zm >= rm) {
 246                         zm = zm - rm;
 247                 } else {
 248                         zm = rm - zm;
 249                         zs = rs;
 250                 }
 251                 if (zm == 0)
 252                         return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
 253
 254                 /*
 255                  * Normalize in extended single precision
 256                  */
 257                 while ((zm >> (SP_MBITS + 3)) == 0) {
 258                         zm <<= 1;
 259                         ze--;
 260                 }
 261
 262         }
 263         return ieee754sp_format(zs, ze, zm);
 264 }
 265
 266 union ieee754sp ieee754sp_maddf(union ieee754sp z, union ieee754sp x,
 267                                 union ieee754sp y)
 268 {
 269         return _sp_maddf(z, x, y, 0);
 270 }
 271
 272 union ieee754sp ieee754sp_msubf(union ieee754sp z, union ieee754sp x,
 273                                 union ieee754sp y)
 274 {
 275         return _sp_maddf(z, x, y, maddf_negate_product);
 276 }