softfloat/s_mulAddF64.c

   1
   2 /*============================================================================
   3
   4 This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
   5 Package, Release 3a, by John R. Hauser.
   6
   7 Copyright 2011, 2012, 2013, 2014 The Regents of the University of California.
   8 All rights reserved.
   9
  10 Redistribution and use in source and binary forms, with or without
  11 modification, are permitted provided that the following conditions are met:
  12
  13  1. Redistributions of source code must retain the above copyright notice,
  14     this list of conditions, and the following disclaimer.
  15
  16  2. Redistributions in binary form must reproduce the above copyright notice,
  17     this list of conditions, and the following disclaimer in the documentation
  18     and/or other materials provided with the distribution.
  19
  20  3. Neither the name of the University nor the names of its contributors may
  21     be used to endorse or promote products derived from this software without
  22     specific prior written permission.
  23
  24 THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
  25 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  26 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
  27 DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
  28 DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  29 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  30 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  31 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  32 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  33 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  34
  35 =============================================================================*/
  36
  37 #include <stdbool.h>
  38 #include <stdint.h>
  39 #include "platform.h"
  40 #include "internals.h"
  41 #include "specialize.h"
  42 #include "softfloat.h"
  43
  44 #ifdef SOFTFLOAT_FAST_INT64
  45
  46 float64_t
  47  softfloat_mulAddF64(
  48      uint_fast64_t uiA, uint_fast64_t uiB, uint_fast64_t uiC, uint_fast8_t op )
  49 {
  50     bool signA;
  51     int_fast16_t expA;
  52     uint_fast64_t sigA;
  53     bool signB;
  54     int_fast16_t expB;
  55     uint_fast64_t sigB;
  56     bool signC;
  57     int_fast16_t expC;
  58     uint_fast64_t sigC;
  59     bool signZ;
  60     uint_fast64_t magBits, uiZ;
  61     struct exp16_sig64 normExpSig;
  62     int_fast16_t expZ;
  63     struct uint128 sig128Z;
  64     uint_fast64_t sigZ;
  65     int_fast16_t expDiff;
  66     struct uint128 sig128C;
  67     int_fast8_t shiftCount;
  68     union ui64_f64 uZ;
  69
  70     /*------------------------------------------------------------------------
  71     *------------------------------------------------------------------------*/
  72     signA = signF64UI( uiA );
  73     expA  = expF64UI( uiA );
  74     sigA  = fracF64UI( uiA );
  75     signB = signF64UI( uiB );
  76     expB  = expF64UI( uiB );
  77     sigB  = fracF64UI( uiB );
  78     signC = signF64UI( uiC ) ^ (op == softfloat_mulAdd_subC);
  79     expC  = expF64UI( uiC );
  80     sigC  = fracF64UI( uiC );
  81     signZ = signA ^ signB ^ (op == softfloat_mulAdd_subProd);
  82     /*------------------------------------------------------------------------
  83     *------------------------------------------------------------------------*/
  84     if ( expA == 0x7FF ) {
  85         if ( sigA || ((expB == 0x7FF) && sigB) ) goto propagateNaN_ABC;
  86         magBits = expB | sigB;
  87         goto infProdArg;
  88     }
  89     if ( expB == 0x7FF ) {
  90         if ( sigB ) goto propagateNaN_ABC;
  91         magBits = expA | sigA;
  92         goto infProdArg;
  93     }
  94     if ( expC == 0x7FF ) {
  95         if ( sigC ) {
  96             uiZ = 0;
  97             goto propagateNaN_ZC;
  98         }
  99         uiZ = uiC;
 100         goto uiZ;
 101     }
 102     /*------------------------------------------------------------------------
 103     *------------------------------------------------------------------------*/
 104     if ( ! expA ) {
 105         if ( ! sigA ) goto zeroProd;
 106         normExpSig = softfloat_normSubnormalF64Sig( sigA );
 107         expA = normExpSig.exp;
 108         sigA = normExpSig.sig;
 109     }
 110     if ( ! expB ) {
 111         if ( ! sigB ) goto zeroProd;
 112         normExpSig = softfloat_normSubnormalF64Sig( sigB );
 113         expB = normExpSig.exp;
 114         sigB = normExpSig.sig;
 115     }
 116     /*------------------------------------------------------------------------
 117     *------------------------------------------------------------------------*/
 118     expZ = expA + expB - 0x3FE;
 119     sigA = (sigA | UINT64_C( 0x0010000000000000 ))<<10;
 120     sigB = (sigB | UINT64_C( 0x0010000000000000 ))<<10;
 121     sig128Z = softfloat_mul64To128( sigA, sigB );
 122     if ( sig128Z.v64 < UINT64_C( 0x2000000000000000 ) ) {
 123         --expZ;
 124         sig128Z =
 125             softfloat_add128(
 126                 sig128Z.v64, sig128Z.v0, sig128Z.v64, sig128Z.v0 );
 127     }
 128     if ( ! expC ) {
 129         if ( ! sigC ) {
 130             --expZ;
 131             sigZ = sig128Z.v64<<1 | (sig128Z.v0 != 0);
 132             goto roundPack;
 133         }
 134         normExpSig = softfloat_normSubnormalF64Sig( sigC );
 135         expC = normExpSig.exp;
 136         sigC = normExpSig.sig;
 137     }
 138     sigC = (sigC | UINT64_C( 0x0010000000000000 ))<<9;
 139     /*------------------------------------------------------------------------
 140     *------------------------------------------------------------------------*/
 141     expDiff = expZ - expC;
 142     if ( expDiff < 0 ) {
 143         expZ = expC;
 144         if ( (signZ == signC) || (expDiff < -1) ) {
 145             sig128Z.v64 = softfloat_shiftRightJam64( sig128Z.v64, -expDiff );
 146         } else {
 147             sig128Z =
 148                 softfloat_shortShiftRightJam128( sig128Z.v64, sig128Z.v0, 1 );
 149         }
 150     } else if ( expDiff ) {
 151         sig128C = softfloat_shiftRightJam128( sigC, 0, expDiff );
 152     }
 153     /*------------------------------------------------------------------------
 154     *------------------------------------------------------------------------*/
 155     if ( signZ == signC ) {
 156         /*--------------------------------------------------------------------
 157         *--------------------------------------------------------------------*/
 158         if ( expDiff <= 0 ) {
 159             sigZ = (sigC + sig128Z.v64) | (sig128Z.v0 != 0);
 160         } else {
 161             sig128Z =
 162                 softfloat_add128(
 163                     sig128Z.v64, sig128Z.v0, sig128C.v64, sig128C.v0 );
 164             sigZ = sig128Z.v64 | (sig128Z.v0 != 0);
 165         }
 166         if ( sigZ < UINT64_C( 0x4000000000000000 ) ) {
 167             --expZ;
 168             sigZ <<= 1;
 169         }
 170     } else {
 171         /*--------------------------------------------------------------------
 172         *--------------------------------------------------------------------*/
 173         if ( expDiff < 0 ) {
 174             signZ = signC;
 175             sig128Z = softfloat_sub128( sigC, 0, sig128Z.v64, sig128Z.v0 );
 176         } else if ( ! expDiff ) {
 177             sig128Z.v64 = sig128Z.v64 - sigC;
 178             if ( ! (sig128Z.v64 | sig128Z.v0) ) goto completeCancellation;
 179             if ( sig128Z.v64 & UINT64_C( 0x8000000000000000 ) ) {
 180                 signZ ^= 1;
 181                 sig128Z = softfloat_sub128( 0, 0, sig128Z.v64, sig128Z.v0 );
 182             }
 183         } else {
 184             sig128Z =
 185                 softfloat_sub128(
 186                     sig128Z.v64, sig128Z.v0, sig128C.v64, sig128C.v0 );
 187         }
 188         /*--------------------------------------------------------------------
 189         *--------------------------------------------------------------------*/
 190         if ( ! sig128Z.v64 ) {
 191             expZ -= 64;
 192             sig128Z.v64 = sig128Z.v0;
 193             sig128Z.v0 = 0;
 194         }
 195         shiftCount = softfloat_countLeadingZeros64( sig128Z.v64 ) - 1;
 196         expZ -= shiftCount;
 197         if ( shiftCount < 0 ) {
 198             sigZ = softfloat_shortShiftRightJam64( sig128Z.v64, -shiftCount );
 199         } else {
 200             sig128Z =
 201                 softfloat_shortShiftLeft128(
 202                     sig128Z.v64, sig128Z.v0, shiftCount );
 203             sigZ = sig128Z.v64;
 204         }
 205         sigZ |= (sig128Z.v0 != 0);
 206     }
 207  roundPack:
 208     return softfloat_roundPackToF64( signZ, expZ, sigZ );
 209     /*------------------------------------------------------------------------
 210     *------------------------------------------------------------------------*/
 211  propagateNaN_ABC:
 212     uiZ = softfloat_propagateNaNF64UI( uiA, uiB );
 213     goto propagateNaN_ZC;
 214     /*------------------------------------------------------------------------
 215     *------------------------------------------------------------------------*/
 216  infProdArg:
 217     if ( magBits ) {
 218         uiZ = packToF64UI( signZ, 0x7FF, 0 );
 219         if ( expC != 0x7FF ) goto uiZ;
 220         if ( sigC ) goto propagateNaN_ZC;
 221         if ( signZ == signC ) goto uiZ;
 222     }
 223  invalid:
 224     softfloat_raiseFlags( softfloat_flag_invalid );
 225     uiZ = defaultNaNF64UI;
 226  propagateNaN_ZC:
 227     uiZ = softfloat_propagateNaNF64UI( uiZ, uiC );
 228     goto uiZ;
 229     /*------------------------------------------------------------------------
 230     *------------------------------------------------------------------------*/
 231  zeroProd:
 232     uiZ = uiC;
 233     if ( ! (expC | sigC) && (signZ != signC) ) {
 234  completeCancellation:
 235         uiZ =
 236             packToF64UI( softfloat_roundingMode == softfloat_round_min, 0, 0 );
 237     }
 238  uiZ:
 239     uZ.ui = uiZ;
 240     return uZ.f;
 241
 242 }
 243
 244 #else
 245
 246 float64_t
 247  softfloat_mulAddF64(
 248      uint_fast64_t uiA, uint_fast64_t uiB, uint_fast64_t uiC, uint_fast8_t op )
 249 {
 250     bool signA;
 251     int_fast16_t expA;
 252     uint64_t sigA;
 253     bool signB;
 254     int_fast16_t expB;
 255     uint64_t sigB;
 256     bool signC;
 257     int_fast16_t expC;
 258     uint64_t sigC;
 259     bool signZ;
 260     uint64_t magBits, uiZ;
 261     struct exp16_sig64 normExpSig;
 262     int_fast16_t expZ;
 263     uint32_t sig128Z[4];
 264     uint64_t sigZ;
 265     int_fast16_t shiftCount, expDiff;
 266     uint32_t sig128C[4];
 267     union ui64_f64 uZ;
 268
 269     /*------------------------------------------------------------------------
 270     *------------------------------------------------------------------------*/
 271     signA = signF64UI( uiA );
 272     expA  = expF64UI( uiA );
 273     sigA  = fracF64UI( uiA );
 274     signB = signF64UI( uiB );
 275     expB  = expF64UI( uiB );
 276     sigB  = fracF64UI( uiB );
 277     signC = signF64UI( uiC ) ^ (op == softfloat_mulAdd_subC);
 278     expC  = expF64UI( uiC );
 279     sigC  = fracF64UI( uiC );
 280     signZ = signA ^ signB ^ (op == softfloat_mulAdd_subProd);
 281     /*------------------------------------------------------------------------
 282     *------------------------------------------------------------------------*/
 283     if ( expA == 0x7FF ) {
 284         if ( sigA || ((expB == 0x7FF) && sigB) ) goto propagateNaN_ABC;
 285         magBits = expB | sigB;
 286         goto infProdArg;
 287     }
 288     if ( expB == 0x7FF ) {
 289         if ( sigB ) goto propagateNaN_ABC;
 290         magBits = expA | sigA;
 291         goto infProdArg;
 292     }
 293     if ( expC == 0x7FF ) {
 294         if ( sigC ) {
 295             uiZ = 0;
 296             goto propagateNaN_ZC;
 297         }
 298         uiZ = uiC;
 299         goto uiZ;
 300     }
 301     /*------------------------------------------------------------------------
 302     *------------------------------------------------------------------------*/
 303     if ( ! expA ) {
 304         if ( ! sigA ) goto zeroProd;
 305         normExpSig = softfloat_normSubnormalF64Sig( sigA );
 306         expA = normExpSig.exp;
 307         sigA = normExpSig.sig;
 308     }
 309     if ( ! expB ) {
 310         if ( ! sigB ) goto zeroProd;
 311         normExpSig = softfloat_normSubnormalF64Sig( sigB );
 312         expB = normExpSig.exp;
 313         sigB = normExpSig.sig;
 314     }
 315     /*------------------------------------------------------------------------
 316     *------------------------------------------------------------------------*/
 317     expZ = expA + expB - 0x3FE;
 318     sigA = (sigA | UINT64_C( 0x0010000000000000 ))<<10;
 319     sigB = (sigB | UINT64_C( 0x0010000000000000 ))<<11;
 320     softfloat_mul64To128M( sigA, sigB, sig128Z );
 321     sigZ =
 322         (uint64_t) sig128Z[indexWord( 4, 3 )]<<32 | sig128Z[indexWord( 4, 2 )];
 323     shiftCount = 0;
 324     if ( ! (sigZ & UINT64_C( 0x4000000000000000 )) ) {
 325         --expZ;
 326         shiftCount = -1;
 327     }
 328     if ( ! expC ) {
 329         if ( ! sigC ) {
 330             if ( shiftCount ) sigZ <<= 1;
 331             goto sigZ;
 332         }
 333         normExpSig = softfloat_normSubnormalF64Sig( sigC );
 334         expC = normExpSig.exp;
 335         sigC = normExpSig.sig;
 336     }
 337     sigC = (sigC | UINT64_C( 0x0010000000000000 ))<<10;
 338     /*------------------------------------------------------------------------
 339     *------------------------------------------------------------------------*/
 340     expDiff = expZ - expC;
 341     if ( expDiff < 0 ) {
 342         expZ = expC;
 343         if ( (signZ == signC) || (expDiff < -1) ) {
 344             shiftCount -= expDiff;
 345             if ( shiftCount) {
 346                 sigZ = softfloat_shiftRightJam64( sigZ, shiftCount );
 347             }
 348         } else {
 349             if ( ! shiftCount ) {
 350                 softfloat_shortShiftRight128M( sig128Z, 1, sig128Z );
 351             }
 352         }
 353     } else {
 354         if ( shiftCount ) softfloat_add128M( sig128Z, sig128Z, sig128Z );
 355         if ( ! expDiff ) {
 356             sigZ =
 357                 (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
 358                     | sig128Z[indexWord( 4, 2 )];
 359         } else {
 360             sig128C[indexWord( 4, 3 )] = sigC>>32;
 361             sig128C[indexWord( 4, 2 )] = sigC;
 362             sig128C[indexWord( 4, 1 )] = 0;
 363             sig128C[indexWord( 4, 0 )] = 0;
 364             softfloat_shiftRightJam128M( sig128C, expDiff, sig128C );
 365         }
 366     }
 367     /*------------------------------------------------------------------------
 368     *------------------------------------------------------------------------*/
 369     if ( signZ == signC ) {
 370         /*--------------------------------------------------------------------
 371         *--------------------------------------------------------------------*/
 372         if ( expDiff <= 0 ) {
 373             sigZ += sigC;
 374         } else {
 375             softfloat_add128M( sig128Z, sig128C, sig128Z );
 376             sigZ =
 377                 (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
 378                     | sig128Z[indexWord( 4, 2 )];
 379         }
 380         if ( sigZ & UINT64_C( 0x8000000000000000 ) ) {
 381             ++expZ;
 382             sigZ = softfloat_shortShiftRightJam64( sigZ, 1 );
 383         }
 384     } else {
 385         /*--------------------------------------------------------------------
 386         *--------------------------------------------------------------------*/
 387         if ( expDiff < 0 ) {
 388             signZ = signC;
 389             if ( expDiff < -1 ) {
 390                 sigZ = sigC - sigZ;
 391                 if (
 392                     sig128Z[indexWord( 4, 1 )] || sig128Z[indexWord( 4, 0 )]
 393                 ) {
 394                     sigZ = (sigZ - 1) | 1;
 395                 }
 396                 if ( ! (sigZ & UINT64_C( 0x4000000000000000 )) ) {
 397                     --expZ;
 398                     sigZ <<= 1;
 399                 }
 400                 goto roundPack;
 401             } else {
 402                 sig128C[indexWord( 4, 3 )] = sigC>>32;
 403                 sig128C[indexWord( 4, 2 )] = sigC;
 404                 sig128C[indexWord( 4, 1 )] = 0;
 405                 sig128C[indexWord( 4, 0 )] = 0;
 406                 softfloat_sub128M( sig128C, sig128Z, sig128Z );
 407             }
 408         } else if ( ! expDiff ) {
 409             sigZ -= sigC;
 410             if (
 411                 ! sigZ && ! sig128Z[indexWord( 4, 1 )]
 412                     && ! sig128Z[indexWord( 4, 0 )]
 413             ) {
 414                 goto completeCancellation;
 415             }
 416             sig128Z[indexWord( 4, 3 )] = sigZ>>32;
 417             sig128Z[indexWord( 4, 2 )] = sigZ;
 418             if ( sigZ & UINT64_C( 0x8000000000000000 ) ) {
 419                 signZ ^= 1;
 420                 softfloat_negX128M( sig128Z );
 421             }
 422         } else {
 423             softfloat_sub128M( sig128Z, sig128C, sig128Z );
 424             if ( 1 < expDiff ) {
 425                 sigZ =
 426                     (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
 427                         | sig128Z[indexWord( 4, 2 )];
 428                 if ( ! (sigZ & UINT64_C( 0x4000000000000000 )) ) {
 429                     --expZ;
 430                     sigZ <<= 1;
 431                 }
 432                 goto sigZ;
 433             }
 434         }
 435         /*--------------------------------------------------------------------
 436         *--------------------------------------------------------------------*/
 437         shiftCount = 0;
 438         sigZ =
 439             (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
 440                 | sig128Z[indexWord( 4, 2 )];
 441         if ( ! sigZ ) {
 442             shiftCount = 64;
 443             sigZ =
 444                 (uint64_t) sig128Z[indexWord( 4, 1 )]<<32
 445                     | sig128Z[indexWord( 4, 0 )];
 446         }
 447         shiftCount += softfloat_countLeadingZeros64( sigZ ) - 1;
 448         if ( shiftCount ) {
 449             expZ -= shiftCount;
 450             softfloat_shiftLeft128M( sig128Z, shiftCount, sig128Z );
 451             sigZ =
 452                 (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
 453                     | sig128Z[indexWord( 4, 2 )];
 454         }
 455     }
 456  sigZ:
 457     if ( sig128Z[indexWord( 4, 1 )] || sig128Z[indexWord( 4, 0 )] ) sigZ |= 1;
 458  roundPack:
 459     return softfloat_roundPackToF64( signZ, expZ - 1, sigZ );
 460     /*------------------------------------------------------------------------
 461     *------------------------------------------------------------------------*/
 462  propagateNaN_ABC:
 463     uiZ = softfloat_propagateNaNF64UI( uiA, uiB );
 464     goto propagateNaN_ZC;
 465     /*------------------------------------------------------------------------
 466     *------------------------------------------------------------------------*/
 467  infProdArg:
 468     if ( magBits ) {
 469         uiZ = packToF64UI( signZ, 0x7FF, 0 );
 470         if ( expC != 0x7FF ) goto uiZ;
 471         if ( sigC ) goto propagateNaN_ZC;
 472         if ( signZ == signC ) goto uiZ;
 473     }
 474  invalid:
 475     softfloat_raiseFlags( softfloat_flag_invalid );
 476     uiZ = defaultNaNF64UI;
 477  propagateNaN_ZC:
 478     uiZ = softfloat_propagateNaNF64UI( uiZ, uiC );
 479     goto uiZ;
 480     /*------------------------------------------------------------------------
 481     *------------------------------------------------------------------------*/
 482  zeroProd:
 483     uiZ = uiC;
 484     if ( ! (expC | sigC) && (signZ != signC) ) {
 485  completeCancellation:
 486         uiZ =
 487             packToF64UI( softfloat_roundingMode == softfloat_round_min, 0, 0 );
 488     }
 489  uiZ:
 490     uZ.ui = uiZ;
 491     return uZ.f;
 492
 493 }
 494
 495 #endif
 496