uspace/softfloat/generic/div.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<add.h>
  37 #include<div.h>
  38 #include<comparison.h>
  39 #include<mul.h>
  40 #include<common.h>
  41
  42
  43 float32 divFloat32(float32 a, float32 b)
  44 {
  45         float32 result;
  46         int32_t aexp, bexp, cexp;
  47         uint64_t afrac, bfrac, cfrac;
  48
  49         result.parts.sign = a.parts.sign ^ b.parts.sign;
  50
  51         if (isFloat32NaN(a)) {
  52                 if (isFloat32SigNaN(a)) {
  53                         /*FIXME: SigNaN*/
  54                 }
  55                 /*NaN*/
  56                 return a;
  57         }
  58
  59         if (isFloat32NaN(b)) {
  60                 if (isFloat32SigNaN(b)) {
  61                         /*FIXME: SigNaN*/
  62                 }
  63                 /*NaN*/
  64                 return b;
  65         }
  66
  67         if (isFloat32Infinity(a)) {
  68                 if (isFloat32Infinity(b)) {
  69                         /*FIXME: inf / inf */
  70                         result.binary = FLOAT32_NAN;
  71                         return result;
  72                 }
  73                 /* inf / num */
  74                 result.parts.exp = a.parts.exp;
  75                 result.parts.fraction = a.parts.fraction;
  76                 return result;
  77         }
  78
  79         if (isFloat32Infinity(b)) {
  80                 if (isFloat32Zero(a)) {
  81                         /* FIXME 0 / inf */
  82                         result.parts.exp = 0;
  83                         result.parts.fraction = 0;
  84                         return result;
  85                 }
  86                 /* FIXME: num / inf*/
  87                 result.parts.exp = 0;
  88                 result.parts.fraction = 0;
  89                 return result;
  90         }
  91
  92         if (isFloat32Zero(b)) {
  93                 if (isFloat32Zero(a)) {
  94                         /*FIXME: 0 / 0*/
  95                         result.binary = FLOAT32_NAN;
  96                         return result;
  97                 }
  98                 /* FIXME: division by zero */
  99                 result.parts.exp = 0;
 100                 result.parts.fraction = 0;
 101                 return result;
 102         }
 103
 104
 105         afrac = a.parts.fraction;
 106         aexp = a.parts.exp;
 107         bfrac = b.parts.fraction;
 108         bexp = b.parts.exp;
 109
 110         /* denormalized numbers */
 111         if (aexp == 0) {
 112                 if (afrac == 0) {
 113                 result.parts.exp = 0;
 114                 result.parts.fraction = 0;
 115                 return result;
 116                 }
 117                 /* normalize it*/
 118
 119                 afrac <<= 1;
 120                         /* afrac is nonzero => it must stop */
 121                 while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
 122                         afrac <<= 1;
 123                         aexp--;
 124                 }
 125         }
 126
 127         if (bexp == 0) {
 128                 bfrac <<= 1;
 129                         /* bfrac is nonzero => it must stop */
 130                 while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
 131                         bfrac <<= 1;
 132                         bexp--;
 133                 }
 134         }
 135
 136         afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
 137         bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
 138
 139         if ( bfrac <= (afrac << 1) ) {
 140                 afrac >>= 1;
 141                 aexp++;
 142         }
 143
 144         cexp = aexp - bexp + FLOAT32_BIAS - 2;
 145
 146         cfrac = (afrac << 32) / bfrac;
 147         if ((  cfrac & 0x3F ) == 0) {
 148                 cfrac |= ( bfrac * cfrac != afrac << 32 );
 149         }
 150
 151         /* pack and round */
 152
 153         /* find first nonzero digit and shift result and detect possibly underflow */
 154         while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
 155                 cexp--;
 156                 cfrac <<= 1;
 157                         /* TODO: fix underflow */
 158         };
 159
 160         cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
 161
 162         if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
 163                 ++cexp;
 164                 cfrac >>= 1;
 165                 }
 166
 167         /* check overflow */
 168         if (cexp >= FLOAT32_MAX_EXPONENT ) {
 169                 /* FIXME: overflow, return infinity */
 170                 result.parts.exp = FLOAT32_MAX_EXPONENT;
 171                 result.parts.fraction = 0;
 172                 return result;
 173         }
 174
 175         if (cexp < 0) {
 176                 /* FIXME: underflow */
 177                 result.parts.exp = 0;
 178                 if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
 179                         result.parts.fraction = 0;
 180                         return result;
 181                 }
 182                 cfrac >>= 1;
 183                 while (cexp < 0) {
 184                         cexp ++;
 185                         cfrac >>= 1;
 186                 }
 187
 188         } else {
 189                 result.parts.exp = (uint32_t)cexp;
 190         }
 191
 192         result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
 193
 194         return result;
 195 }
 196
 197 float64 divFloat64(float64 a, float64 b)
 198 {
 199         float64 result;
 200         int64_t aexp, bexp, cexp;
 201         uint64_t afrac, bfrac, cfrac;
 202         uint64_t remlo, remhi;
 203
 204         result.parts.sign = a.parts.sign ^ b.parts.sign;
 205
 206         if (isFloat64NaN(a)) {
 207
 208                 if (isFloat64SigNaN(b)) {
 209                         /*FIXME: SigNaN*/
 210                         return b;
 211                 }
 212
 213                 if (isFloat64SigNaN(a)) {
 214                         /*FIXME: SigNaN*/
 215                 }
 216                 /*NaN*/
 217                 return a;
 218         }
 219
 220         if (isFloat64NaN(b)) {
 221                 if (isFloat64SigNaN(b)) {
 222                         /*FIXME: SigNaN*/
 223                 }
 224                 /*NaN*/
 225                 return b;
 226         }
 227
 228         if (isFloat64Infinity(a)) {
 229                 if (isFloat64Infinity(b) || isFloat64Zero(b)) {
 230                         /*FIXME: inf / inf */
 231                         result.binary = FLOAT64_NAN;
 232                         return result;
 233                 }
 234                 /* inf / num */
 235                 result.parts.exp = a.parts.exp;
 236                 result.parts.fraction = a.parts.fraction;
 237                 return result;
 238         }
 239
 240         if (isFloat64Infinity(b)) {
 241                 if (isFloat64Zero(a)) {
 242                         /* FIXME 0 / inf */
 243                         result.parts.exp = 0;
 244                         result.parts.fraction = 0;
 245                         return result;
 246                 }
 247                 /* FIXME: num / inf*/
 248                 result.parts.exp = 0;
 249                 result.parts.fraction = 0;
 250                 return result;
 251         }
 252
 253         if (isFloat64Zero(b)) {
 254                 if (isFloat64Zero(a)) {
 255                         /*FIXME: 0 / 0*/
 256                         result.binary = FLOAT64_NAN;
 257                         return result;
 258                 }
 259                 /* FIXME: division by zero */
 260                 result.parts.exp = 0;
 261                 result.parts.fraction = 0;
 262                 return result;
 263         }
 264
 265
 266         afrac = a.parts.fraction;
 267         aexp = a.parts.exp;
 268         bfrac = b.parts.fraction;
 269         bexp = b.parts.exp;
 270
 271         /* denormalized numbers */
 272         if (aexp == 0) {
 273                 if (afrac == 0) {
 274                         result.parts.exp = 0;
 275                         result.parts.fraction = 0;
 276                         return result;
 277                 }
 278                 /* normalize it*/
 279
 280                 aexp++;
 281                         /* afrac is nonzero => it must stop */
 282                 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
 283                         afrac <<= 1;
 284                         aexp--;
 285                 }
 286         }
 287
 288         if (bexp == 0) {
 289                 bexp++;
 290                         /* bfrac is nonzero => it must stop */
 291                 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
 292                         bfrac <<= 1;
 293                         bexp--;
 294                 }
 295         }
 296
 297         afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
 298         bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
 299
 300         if ( bfrac <= (afrac << 1) ) {
 301                 afrac >>= 1;
 302                 aexp++;
 303         }
 304
 305         cexp = aexp - bexp + FLOAT64_BIAS - 2;
 306
 307         cfrac = divFloat64estim(afrac, bfrac);
 308
 309         if ((  cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
 310                 mul64integers( bfrac, cfrac, &remlo, &remhi);
 311                 /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/
 312                 remhi = afrac - remhi - ( remlo > 0);
 313                 remlo = - remlo;
 314
 315                 while ((int64_t) remhi < 0) {
 316                         cfrac--;
 317                         remlo += bfrac;
 318                         remhi += ( remlo < bfrac );
 319                 }
 320                 cfrac |= ( remlo != 0 );
 321         }
 322
 323         /* round and shift */
 324         result = finishFloat64(cexp, cfrac, result.parts.sign);
 325         return result;
 326
 327 }
 328
 329 uint64_t divFloat64estim(uint64_t a, uint64_t b)
 330 {
 331         uint64_t bhi;
 332         uint64_t remhi, remlo;
 333         uint64_t result;
 334
 335         if ( b <= a ) {
 336                 return 0xFFFFFFFFFFFFFFFFull;
 337         }
 338
 339         bhi = b >> 32;
 340         result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
 341         mul64integers(b, result, &remlo, &remhi);
 342
 343         remhi = a - remhi - (remlo > 0);
 344         remlo = - remlo;
 345
 346         b <<= 32;
 347         while ( (int64_t) remhi < 0 ) {
 348                         result -= 0x1ll << 32;
 349                         remlo += b;
 350                         remhi += bhi + ( remlo < b );
 351                 }
 352         remhi = (remhi << 32) | (remlo >> 32);
 353         if (( bhi << 32) <= remhi) {
 354                 result |= 0xFFFFFFFF;
 355         } else {
 356                 result |= remhi / bhi;
 357         }
 358
 359
 360         return result;
 361 }
 362
 363 /** @}
 364  */