softfloat/s_mulAddF128.c

   1
   2 /*============================================================================
   3
   4 This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
   5 Package, Release 3d, by John R. Hauser.
   6
   7 Copyright 2011, 2012, 2013, 2014, 2015, 2016 The Regents of the University of
   8 California.  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 float128_t
  45  softfloat_mulAddF128(
  46      uint_fast64_t uiA64,
  47      uint_fast64_t uiA0,
  48      uint_fast64_t uiB64,
  49      uint_fast64_t uiB0,
  50      uint_fast64_t uiC64,
  51      uint_fast64_t uiC0,
  52      uint_fast8_t op
  53  )
  54 {
  55     bool signA;
  56     int_fast32_t expA;
  57     struct uint128 sigA;
  58     bool signB;
  59     int_fast32_t expB;
  60     struct uint128 sigB;
  61     bool signC;
  62     int_fast32_t expC;
  63     struct uint128 sigC;
  64     bool signZ;
  65     uint_fast64_t magBits;
  66     struct uint128 uiZ;
  67     struct exp32_sig128 normExpSig;
  68     int_fast32_t expZ;
  69     uint64_t sig256Z[4];
  70     struct uint128 sigZ;
  71     int_fast32_t shiftDist, expDiff;
  72     struct uint128 x128;
  73     uint64_t sig256C[4];
  74     static uint64_t zero256[4] = INIT_UINTM4( 0, 0, 0, 0 );
  75     uint_fast64_t sigZExtra, sig256Z0;
  76     union ui128_f128 uZ;
  77
  78     /*------------------------------------------------------------------------
  79     *------------------------------------------------------------------------*/
  80     signA = signF128UI64( uiA64 );
  81     expA  = expF128UI64( uiA64 );
  82     sigA.v64 = fracF128UI64( uiA64 );
  83     sigA.v0  = uiA0;
  84     signB = signF128UI64( uiB64 );
  85     expB  = expF128UI64( uiB64 );
  86     sigB.v64 = fracF128UI64( uiB64 );
  87     sigB.v0  = uiB0;
  88     signC = signF128UI64( uiC64 ) ^ (op == softfloat_mulAdd_subC);
  89     expC  = expF128UI64( uiC64 );
  90     sigC.v64 = fracF128UI64( uiC64 );
  91     sigC.v0  = uiC0;
  92     signZ = signA ^ signB ^ (op == softfloat_mulAdd_subProd);
  93     /*------------------------------------------------------------------------
  94     *------------------------------------------------------------------------*/
  95     if ( expA == 0x7FFF ) {
  96         if (
  97             (sigA.v64 | sigA.v0) || ((expB == 0x7FFF) && (sigB.v64 | sigB.v0))
  98         ) {
  99             goto propagateNaN_ABC;
 100         }
 101         magBits = expB | sigB.v64 | sigB.v0;
 102         goto infProdArg;
 103     }
 104     if ( expB == 0x7FFF ) {
 105         if ( sigB.v64 | sigB.v0 ) goto propagateNaN_ABC;
 106         magBits = expA | sigA.v64 | sigA.v0;
 107         goto infProdArg;
 108     }
 109     if ( expC == 0x7FFF ) {
 110         if ( sigC.v64 | sigC.v0 ) {
 111             uiZ.v64 = 0;
 112             uiZ.v0  = 0;
 113             goto propagateNaN_ZC;
 114         }
 115         uiZ.v64 = uiC64;
 116         uiZ.v0  = uiC0;
 117         goto uiZ;
 118     }
 119     /*------------------------------------------------------------------------
 120     *------------------------------------------------------------------------*/
 121     if ( ! expA ) {
 122         if ( ! (sigA.v64 | sigA.v0) ) goto zeroProd;
 123         normExpSig = softfloat_normSubnormalF128Sig( sigA.v64, sigA.v0 );
 124         expA = normExpSig.exp;
 125         sigA = normExpSig.sig;
 126     }
 127     if ( ! expB ) {
 128         if ( ! (sigB.v64 | sigB.v0) ) goto zeroProd;
 129         normExpSig = softfloat_normSubnormalF128Sig( sigB.v64, sigB.v0 );
 130         expB = normExpSig.exp;
 131         sigB = normExpSig.sig;
 132     }
 133     /*------------------------------------------------------------------------
 134     *------------------------------------------------------------------------*/
 135     expZ = expA + expB - 0x3FFE;
 136     sigA.v64 |= UINT64_C( 0x0001000000000000 );
 137     sigB.v64 |= UINT64_C( 0x0001000000000000 );
 138     sigA = softfloat_shortShiftLeft128( sigA.v64, sigA.v0, 8 );
 139     sigB = softfloat_shortShiftLeft128( sigB.v64, sigB.v0, 15 );
 140     softfloat_mul128To256M( sigA.v64, sigA.v0, sigB.v64, sigB.v0, sig256Z );
 141     sigZ.v64 = sig256Z[indexWord( 4, 3 )];
 142     sigZ.v0  = sig256Z[indexWord( 4, 2 )];
 143     shiftDist = 0;
 144     if ( ! (sigZ.v64 & UINT64_C( 0x0100000000000000 )) ) {
 145         --expZ;
 146         shiftDist = -1;
 147     }
 148     if ( ! expC ) {
 149         if ( ! (sigC.v64 | sigC.v0) ) {
 150             shiftDist += 8;
 151             goto sigZ;
 152         }
 153         normExpSig = softfloat_normSubnormalF128Sig( sigC.v64, sigC.v0 );
 154         expC = normExpSig.exp;
 155         sigC = normExpSig.sig;
 156     }
 157     sigC.v64 |= UINT64_C( 0x0001000000000000 );
 158     sigC = softfloat_shortShiftLeft128( sigC.v64, sigC.v0, 8 );
 159     /*------------------------------------------------------------------------
 160     *------------------------------------------------------------------------*/
 161     expDiff = expZ - expC;
 162     if ( expDiff < 0 ) {
 163         expZ = expC;
 164         if ( (signZ == signC) || (expDiff < -1) ) {
 165             shiftDist -= expDiff;
 166             if ( shiftDist ) {
 167                 sigZ =
 168                     softfloat_shiftRightJam128( sigZ.v64, sigZ.v0, shiftDist );
 169             }
 170         } else {
 171             if ( ! shiftDist ) {
 172                 x128 =
 173                     softfloat_shortShiftRight128(
 174                         sig256Z[indexWord( 4, 1 )], sig256Z[indexWord( 4, 0 )],
 175                         1
 176                     );
 177                 sig256Z[indexWord( 4, 1 )] = (sigZ.v0<<63) | x128.v64;
 178                 sig256Z[indexWord( 4, 0 )] = x128.v0;
 179                 sigZ = softfloat_shortShiftRight128( sigZ.v64, sigZ.v0, 1 );
 180                 sig256Z[indexWord( 4, 3 )] = sigZ.v64;
 181                 sig256Z[indexWord( 4, 2 )] = sigZ.v0;
 182             }
 183         }
 184     } else {
 185         if ( shiftDist ) softfloat_add256M( sig256Z, sig256Z, sig256Z );
 186         if ( ! expDiff ) {
 187             sigZ.v64 = sig256Z[indexWord( 4, 3 )];
 188             sigZ.v0  = sig256Z[indexWord( 4, 2 )];
 189         } else {
 190             sig256C[indexWord( 4, 3 )] = sigC.v64;
 191             sig256C[indexWord( 4, 2 )] = sigC.v0;
 192             sig256C[indexWord( 4, 1 )] = 0;
 193             sig256C[indexWord( 4, 0 )] = 0;
 194             softfloat_shiftRightJam256M( sig256C, expDiff, sig256C );
 195         }
 196     }
 197     /*------------------------------------------------------------------------
 198     *------------------------------------------------------------------------*/
 199     shiftDist = 8;
 200     if ( signZ == signC ) {
 201         /*--------------------------------------------------------------------
 202         *--------------------------------------------------------------------*/
 203         if ( expDiff <= 0 ) {
 204             sigZ = softfloat_add128( sigC.v64, sigC.v0, sigZ.v64, sigZ.v0 );
 205         } else {
 206             softfloat_add256M( sig256Z, sig256C, sig256Z );
 207             sigZ.v64 = sig256Z[indexWord( 4, 3 )];
 208             sigZ.v0  = sig256Z[indexWord( 4, 2 )];
 209         }
 210         if ( sigZ.v64 & UINT64_C( 0x0200000000000000 ) ) {
 211             ++expZ;
 212             shiftDist = 9;
 213         }
 214     } else {
 215         /*--------------------------------------------------------------------
 216         *--------------------------------------------------------------------*/
 217         if ( expDiff < 0 ) {
 218             signZ = signC;
 219             if ( expDiff < -1 ) {
 220                 sigZ =
 221                     softfloat_sub128( sigC.v64, sigC.v0, sigZ.v64, sigZ.v0 );
 222                 sigZExtra =
 223                     sig256Z[indexWord( 4, 1 )] | sig256Z[indexWord( 4, 0 )];
 224                 if ( sigZExtra ) {
 225                     sigZ = softfloat_sub128( sigZ.v64, sigZ.v0, 0, 1 );
 226                 }
 227                 if ( ! (sigZ.v64 & UINT64_C( 0x0100000000000000 )) ) {
 228                     --expZ;
 229                     shiftDist = 7;
 230                 }
 231                 goto shiftRightRoundPack;
 232             } else {
 233                 sig256C[indexWord( 4, 3 )] = sigC.v64;
 234                 sig256C[indexWord( 4, 2 )] = sigC.v0;
 235                 sig256C[indexWord( 4, 1 )] = 0;
 236                 sig256C[indexWord( 4, 0 )] = 0;
 237                 softfloat_sub256M( sig256C, sig256Z, sig256Z );
 238             }
 239         } else if ( ! expDiff ) {
 240             sigZ = softfloat_sub128( sigZ.v64, sigZ.v0, sigC.v64, sigC.v0 );
 241             if (
 242                 ! (sigZ.v64 | sigZ.v0) && ! sig256Z[indexWord( 4, 1 )]
 243                     && ! sig256Z[indexWord( 4, 0 )]
 244             ) {
 245                 goto completeCancellation;
 246             }
 247             sig256Z[indexWord( 4, 3 )] = sigZ.v64;
 248             sig256Z[indexWord( 4, 2 )] = sigZ.v0;
 249             if ( sigZ.v64 & UINT64_C( 0x8000000000000000 ) ) {
 250                 signZ = ! signZ;
 251                 softfloat_sub256M( zero256, sig256Z, sig256Z );
 252             }
 253         } else {
 254             softfloat_sub256M( sig256Z, sig256C, sig256Z );
 255             if ( 1 < expDiff ) {
 256                 sigZ.v64 = sig256Z[indexWord( 4, 3 )];
 257                 sigZ.v0  = sig256Z[indexWord( 4, 2 )];
 258                 if ( ! (sigZ.v64 & UINT64_C( 0x0100000000000000 )) ) {
 259                     --expZ;
 260                     shiftDist = 7;
 261                 }
 262                 goto sigZ;
 263             }
 264         }
 265         /*--------------------------------------------------------------------
 266         *--------------------------------------------------------------------*/
 267         sigZ.v64  = sig256Z[indexWord( 4, 3 )];
 268         sigZ.v0   = sig256Z[indexWord( 4, 2 )];
 269         sigZExtra = sig256Z[indexWord( 4, 1 )];
 270         sig256Z0  = sig256Z[indexWord( 4, 0 )];
 271         if ( sigZ.v64 ) {
 272             if ( sig256Z0 ) sigZExtra |= 1;
 273         } else {
 274             expZ -= 64;
 275             sigZ.v64  = sigZ.v0;
 276             sigZ.v0   = sigZExtra;
 277             sigZExtra = sig256Z0;
 278             if ( ! sigZ.v64 ) {
 279                 expZ -= 64;
 280                 sigZ.v64  = sigZ.v0;
 281                 sigZ.v0   = sigZExtra;
 282                 sigZExtra = 0;
 283                 if ( ! sigZ.v64 ) {
 284                     expZ -= 64;
 285                     sigZ.v64 = sigZ.v0;
 286                     sigZ.v0  = 0;
 287                 }
 288             }
 289         }
 290         shiftDist = softfloat_countLeadingZeros64( sigZ.v64 );
 291         expZ += 7 - shiftDist;
 292         shiftDist = 15 - shiftDist;
 293         if ( 0 < shiftDist ) goto shiftRightRoundPack;
 294         if ( shiftDist ) {
 295             shiftDist = -shiftDist;
 296             sigZ = softfloat_shortShiftLeft128( sigZ.v64, sigZ.v0, shiftDist );
 297             x128 = softfloat_shortShiftLeft128( 0, sigZExtra, shiftDist );
 298             sigZ.v0 |= x128.v64;
 299             sigZExtra = x128.v0;
 300         }
 301         goto roundPack;
 302     }
 303  sigZ:
 304     sigZExtra = sig256Z[indexWord( 4, 1 )] | sig256Z[indexWord( 4, 0 )];
 305  shiftRightRoundPack:
 306     sigZExtra = (uint64_t) (sigZ.v0<<(64 - shiftDist)) | (sigZExtra != 0);
 307     sigZ = softfloat_shortShiftRight128( sigZ.v64, sigZ.v0, shiftDist );
 308  roundPack:
 309     return
 310         softfloat_roundPackToF128(
 311             signZ, expZ - 1, sigZ.v64, sigZ.v0, sigZExtra );
 312     /*------------------------------------------------------------------------
 313     *------------------------------------------------------------------------*/
 314  propagateNaN_ABC:
 315     uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
 316     goto propagateNaN_ZC;
 317     /*------------------------------------------------------------------------
 318     *------------------------------------------------------------------------*/
 319  infProdArg:
 320     if ( magBits ) {
 321         uiZ.v64 = packToF128UI64( signZ, 0x7FFF, 0 );
 322         uiZ.v0 = 0;
 323         if ( expC != 0x7FFF ) goto uiZ;
 324         if ( sigC.v64 | sigC.v0 ) goto propagateNaN_ZC;
 325         if ( signZ == signC ) goto uiZ;
 326     }
 327     softfloat_raiseFlags( softfloat_flag_invalid );
 328     uiZ.v64 = defaultNaNF128UI64;
 329     uiZ.v0  = defaultNaNF128UI0;
 330  propagateNaN_ZC:
 331     uiZ = softfloat_propagateNaNF128UI( uiZ.v64, uiZ.v0, uiC64, uiC0 );
 332     goto uiZ;
 333     /*------------------------------------------------------------------------
 334     *------------------------------------------------------------------------*/
 335  zeroProd:
 336     uiZ.v64 = uiC64;
 337     uiZ.v0  = uiC0;
 338     if ( ! (expC | sigC.v64 | sigC.v0) && (signZ != signC) ) {
 339  completeCancellation:
 340         uiZ.v64 =
 341             packToF128UI64(
 342                 (softfloat_roundingMode == softfloat_round_min), 0, 0 );
 343         uiZ.v0 = 0;
 344     }
 345  uiZ:
 346     uZ.ui = uiZ;
 347     return uZ.f;
 348
 349 }
 350