1
2/*============================================================================
3
4This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
5Package, Release 3e, by John R. Hauser.
6
7Copyright 2011, 2012, 2013, 2014 The Regents of the University of California.
8All rights reserved.
9
10Redistribution and use in source and binary forms, with or without
11modification, 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
24THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
25EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
26WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
27DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
28DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
29(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
31ON 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
33SOFTWARE, 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
46void
47 extF80M_rem(
48     const extFloat80_t *aPtr, const extFloat80_t *bPtr, extFloat80_t *zPtr )
49{
50
51    *zPtr = extF80_rem( *aPtr, *bPtr );
52
53}
54
55#else
56
57void
58 extF80M_rem(
59     const extFloat80_t *aPtr, const extFloat80_t *bPtr, extFloat80_t *zPtr )
60{
61    const struct extFloat80M *aSPtr, *bSPtr;
62    struct extFloat80M *zSPtr;
63    uint_fast16_t uiA64;
64    int32_t expA, expB;
65    uint64_t x64;
66    bool signRem;
67    uint64_t sigA;
68    int32_t expDiff;
69    uint32_t rem[3], x[3], sig32B, q, recip32, rem2[3], *remPtr, *altRemPtr;
70    uint32_t *newRemPtr, wordMeanRem;
71
72    /*------------------------------------------------------------------------
73    *------------------------------------------------------------------------*/
74    aSPtr = (const struct extFloat80M *) aPtr;
75    bSPtr = (const struct extFloat80M *) bPtr;
76    zSPtr = (struct extFloat80M *) zPtr;
77    /*------------------------------------------------------------------------
78    *------------------------------------------------------------------------*/
79    uiA64 = aSPtr->signExp;
80    expA = expExtF80UI64( uiA64 );
81    expB = expExtF80UI64( bSPtr->signExp );
82    /*------------------------------------------------------------------------
83    *------------------------------------------------------------------------*/
84    if ( (expA == 0x7FFF) || (expB == 0x7FFF) ) {
85        if ( softfloat_tryPropagateNaNExtF80M( aSPtr, bSPtr, zSPtr ) ) return;
86        if ( expA == 0x7FFF ) goto invalid;
87        /*--------------------------------------------------------------------
88        | If we get here, then argument b is an infinity and `expB' is 0x7FFF;
89        | Doubling `expB' is an easy way to ensure that `expDiff' later is
90        | less than -1, which will result in returning a canonicalized version
91        | of argument a.
92        *--------------------------------------------------------------------*/
93        expB += expB;
94    }
95    /*------------------------------------------------------------------------
96    *------------------------------------------------------------------------*/
97    if ( ! expB ) expB = 1;
98    x64 = bSPtr->signif;
99    if ( ! (x64 & UINT64_C( 0x8000000000000000 )) ) {
100        if ( ! x64 ) goto invalid;
101        expB += softfloat_normExtF80SigM( &x64 );
102    }
103    signRem = signExtF80UI64( uiA64 );
104    if ( ! expA ) expA = 1;
105    sigA = aSPtr->signif;
106    if ( ! (sigA & UINT64_C( 0x8000000000000000 )) ) {
107        if ( ! sigA ) {
108            expA = 0;
109            goto copyA;
110        }
111        expA += softfloat_normExtF80SigM( &sigA );
112    }
113    /*------------------------------------------------------------------------
114    *------------------------------------------------------------------------*/
115    expDiff = expA - expB;
116    if ( expDiff < -1 ) goto copyA;
117    rem[indexWord( 3, 2 )] = sigA>>34;
118    rem[indexWord( 3, 1 )] = sigA>>2;
119    rem[indexWord( 3, 0 )] = (uint32_t) sigA<<30;
120    x[indexWord( 3, 0 )] = (uint32_t) x64<<30;
121    sig32B = x64>>32;
122    x64 >>= 2;
123    x[indexWord( 3, 2 )] = x64>>32;
124    x[indexWord( 3, 1 )] = x64;
125    if ( expDiff < 1 ) {
126        if ( expDiff ) {
127            --expB;
128            softfloat_add96M( x, x, x );
129            q = 0;
130        } else {
131            q = (softfloat_compare96M( x, rem ) <= 0);
132            if ( q ) softfloat_sub96M( rem, x, rem );
133        }
134    } else {
135        recip32 = softfloat_approxRecip32_1( sig32B );
136        expDiff -= 30;
137        for (;;) {
138            x64 = (uint64_t) rem[indexWordHi( 3 )] * recip32;
139            if ( expDiff < 0 ) break;
140            q = (x64 + 0x80000000)>>32;
141            softfloat_remStep96MBy32( rem, 29, x, q, rem );
142            if ( rem[indexWordHi( 3 )] & 0x80000000 ) {
143                softfloat_add96M( rem, x, rem );
144            }
145            expDiff -= 29;
146        }
147        /*--------------------------------------------------------------------
148        | (`expDiff' cannot be less than -29 here.)
149        *--------------------------------------------------------------------*/
150        q = (uint32_t) (x64>>32)>>(~expDiff & 31);
151        softfloat_remStep96MBy32( rem, expDiff + 30, x, q, rem );
152        if ( rem[indexWordHi( 3 )] & 0x80000000 ) {
153            remPtr = rem;
154            altRemPtr = rem2;
155            softfloat_add96M( remPtr, x, altRemPtr );
156            goto selectRem;
157        }
158    }
159    /*------------------------------------------------------------------------
160    *------------------------------------------------------------------------*/
161    remPtr = rem;
162    altRemPtr = rem2;
163    do {
164        ++q;
165        newRemPtr = altRemPtr;
166        softfloat_sub96M( remPtr, x, newRemPtr );
167        altRemPtr = remPtr;
168        remPtr = newRemPtr;
169    } while ( ! (remPtr[indexWordHi( 3 )] & 0x80000000) );
170 selectRem:
171    softfloat_add96M( remPtr, altRemPtr, x );
172    wordMeanRem = x[indexWordHi( 3 )];
173    if (
174        (wordMeanRem & 0x80000000)
175            || (! wordMeanRem && (q & 1) && ! x[indexWord( 3, 0 )]
176                    && ! x[indexWord( 3, 1 )])
177    ) {
178        remPtr = altRemPtr;
179    }
180    if ( remPtr[indexWordHi( 3 )] & 0x80000000 ) {
181        signRem = ! signRem;
182        softfloat_negX96M( remPtr );
183    }
184    softfloat_normRoundPackMToExtF80M( signRem, expB + 2, remPtr, 80, zSPtr );
185    return;
186    /*------------------------------------------------------------------------
187    *------------------------------------------------------------------------*/
188 invalid:
189    softfloat_invalidExtF80M( zSPtr );
190    return;
191    /*------------------------------------------------------------------------
192    *------------------------------------------------------------------------*/
193 copyA:
194    if ( expA < 1 ) {
195        sigA >>= 1 - expA;
196        expA = 0;
197    }
198    zSPtr->signExp = packToExtF80UI64( signRem, expA );
199    zSPtr->signif = sigA;
200
201}
202
203#endif
204
205