uspace/softfloat/generic/sub.c

   1 /*
   2  * Copyright (c) 2005 Josef Cejka
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  *
   9  * - Redistributions of source code must retain the above copyright
  10  *   notice, this list of conditions and the following disclaimer.
  11  * - Redistributions in binary form must reproduce the above copyright
  12  *   notice, this list of conditions and the following disclaimer in the
  13  *   documentation and/or other materials provided with the distribution.
  14  * - The name of the author may not be used to endorse or promote products
  15  *   derived from this software without specific prior written permission.
  16  *
  17  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  18  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  19  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  20  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  21  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  22  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  23  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  24  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  26  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  27  */
  28
  29  /** @addtogroup softfloat
  30  * @{
  31  */
  32 /** @file
  33  */
  34
  35 #include<sftypes.h>
  36 #include<sub.h>
  37 #include<comparison.h>
  38
  39 /** Subtract two float32 numbers with same signs
  40  */
  41 float32 subFloat32(float32 a, float32 b)
  42 {
  43         int expdiff;
  44         uint32_t exp1, exp2, frac1, frac2;
  45         float32 result;
  46
  47         result.f = 0;
  48
  49         expdiff = a.parts.exp - b.parts.exp;
  50         if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.fraction < b.parts.fraction))) {
  51                 if (isFloat32NaN(b)) {
  52                         /* TODO: fix SigNaN */
  53                         if (isFloat32SigNaN(b)) {
  54                         };
  55                         return b;
  56                 };
  57
  58                 if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
  59                         b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
  60                         return b;
  61                 }
  62
  63                 result.parts.sign = !a.parts.sign;
  64
  65                 frac1 = b.parts.fraction;
  66                 exp1 = b.parts.exp;
  67                 frac2 = a.parts.fraction;
  68                 exp2 = a.parts.exp;
  69                 expdiff *= -1;
  70         } else {
  71                 if (isFloat32NaN(a)) {
  72                         /* TODO: fix SigNaN */
  73                         if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
  74                         };
  75                         return a;
  76                 };
  77
  78                 if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
  79                         if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
  80                                 /* inf - inf => nan */
  81                                 /* TODO: fix exception */
  82                                 result.binary = FLOAT32_NAN;
  83                                 return result;
  84                         };
  85                         return a;
  86                 }
  87
  88                 result.parts.sign = a.parts.sign;
  89
  90                 frac1 = a.parts.fraction;
  91                 exp1 = a.parts.exp;
  92                 frac2 = b.parts.fraction;
  93                 exp2 = b.parts.exp;
  94         };
  95
  96         if (exp1 == 0) {
  97                 /* both are denormalized */
  98                 result.parts.fraction = frac1-frac2;
  99                 if (result.parts.fraction > frac1) {
 100                         /* TODO: underflow exception */
 101                         return result;
 102                 };
 103                 result.parts.exp = 0;
 104                 return result;
 105         };
 106
 107         /* add hidden bit */
 108         frac1 |= FLOAT32_HIDDEN_BIT_MASK;
 109
 110         if (exp2 == 0) {
 111                 /* denormalized */
 112                 --expdiff;
 113         } else {
 114                 /* normalized */
 115                 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
 116         };
 117
 118         /* create some space for rounding */
 119         frac1 <<= 6;
 120         frac2 <<= 6;
 121
 122         if (expdiff > FLOAT32_FRACTION_SIZE + 1) {
 123              goto done;
 124              };
 125
 126         frac1 = frac1 - (frac2 >> expdiff);
 127 done:
 128         /* TODO: find first nonzero digit and shift result and detect possibly underflow */
 129         while ((exp1 > 0) && (!(frac1 & (FLOAT32_HIDDEN_BIT_MASK << 6 )))) {
 130                 --exp1;
 131                 frac1 <<= 1;
 132                         /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
 133         };
 134
 135         /* rounding - if first bit after fraction is set then round up */
 136         frac1 += 0x20;
 137
 138         if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
 139                 ++exp1;
 140                 frac1 >>= 1;
 141         };
 142
 143         /*Clear hidden bit and shift */
 144         result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
 145         result.parts.exp = exp1;
 146
 147         return result;
 148 }
 149
 150 /** Subtract two float64 numbers with same signs
 151  */
 152 float64 subFloat64(float64 a, float64 b)
 153 {
 154         int expdiff;
 155         uint32_t exp1, exp2;
 156         uint64_t frac1, frac2;
 157         float64 result;
 158
 159         result.d = 0;
 160
 161         expdiff = a.parts.exp - b.parts.exp;
 162         if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.fraction < b.parts.fraction))) {
 163                 if (isFloat64NaN(b)) {
 164                         /* TODO: fix SigNaN */
 165                         if (isFloat64SigNaN(b)) {
 166                         };
 167                         return b;
 168                 };
 169
 170                 if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
 171                         b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
 172                         return b;
 173                 }
 174
 175                 result.parts.sign = !a.parts.sign;
 176
 177                 frac1 = b.parts.fraction;
 178                 exp1 = b.parts.exp;
 179                 frac2 = a.parts.fraction;
 180                 exp2 = a.parts.exp;
 181                 expdiff *= -1;
 182         } else {
 183                 if (isFloat64NaN(a)) {
 184                         /* TODO: fix SigNaN */
 185                         if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
 186                         };
 187                         return a;
 188                 };
 189
 190                 if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
 191                         if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
 192                                 /* inf - inf => nan */
 193                                 /* TODO: fix exception */
 194                                 result.binary = FLOAT64_NAN;
 195                                 return result;
 196                         };
 197                         return a;
 198                 }
 199
 200                 result.parts.sign = a.parts.sign;
 201
 202                 frac1 = a.parts.fraction;
 203                 exp1 = a.parts.exp;
 204                 frac2 = b.parts.fraction;
 205                 exp2 = b.parts.exp;
 206         };
 207
 208         if (exp1 == 0) {
 209                 /* both are denormalized */
 210                 result.parts.fraction = frac1 - frac2;
 211                 if (result.parts.fraction > frac1) {
 212                         /* TODO: underflow exception */
 213                         return result;
 214                 };
 215                 result.parts.exp = 0;
 216                 return result;
 217         };
 218
 219         /* add hidden bit */
 220         frac1 |= FLOAT64_HIDDEN_BIT_MASK;
 221
 222         if (exp2 == 0) {
 223                 /* denormalized */
 224                 --expdiff;
 225         } else {
 226                 /* normalized */
 227                 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
 228         };
 229
 230         /* create some space for rounding */
 231         frac1 <<= 6;
 232         frac2 <<= 6;
 233
 234         if (expdiff > FLOAT64_FRACTION_SIZE + 1) {
 235              goto done;
 236              };
 237
 238         frac1 = frac1 - (frac2 >> expdiff);
 239 done:
 240         /* TODO: find first nonzero digit and shift result and detect possibly underflow */
 241         while ((exp1 > 0) && (!(frac1 & (FLOAT64_HIDDEN_BIT_MASK << 6 )))) {
 242                 --exp1;
 243                 frac1 <<= 1;
 244                         /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
 245         };
 246
 247         /* rounding - if first bit after fraction is set then round up */
 248         frac1 += 0x20;
 249
 250         if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
 251                 ++exp1;
 252                 frac1 >>= 1;
 253         };
 254
 255         /*Clear hidden bit and shift */
 256         result.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK));
 257         result.parts.exp = exp1;
 258
 259         return result;
 260 }
 261
 262
 263  /** @}
 264  */
 265