1 // SPDX-License-Identifier: BSD-3-Clause
2 
3 /*============================================================================
4 
5 This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
6 Package, Release 3a, by John R. Hauser.
7 
8 Copyright 2011, 2012, 2013, 2014, 2015 The Regents of the University of
9 California.  All rights reserved.
10 
11 Redistribution and use in source and binary forms, with or without
12 modification, are permitted provided that the following conditions are met:
13 
14  1. Redistributions of source code must retain the above copyright notice,
15     this list of conditions, and the following disclaimer.
16 
17  2. Redistributions in binary form must reproduce the above copyright notice,
18     this list of conditions, and the following disclaimer in the documentation
19     and/or other materials provided with the distribution.
20 
21  3. Neither the name of the University nor the names of its contributors may
22     be used to endorse or promote products derived from this software without
23     specific prior written permission.
24 
25 THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
26 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
27 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
28 DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
29 DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
30 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
31 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
32 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
33 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35 
36 =============================================================================*/
37 
38 #include <stdbool.h>
39 #include <stdint.h>
40 #include "platform.h"
41 #include "internals.h"
42 #include "specialize.h"
43 #include "softfloat.h"
44 
45 #ifdef SOFTFLOAT_FAST_INT64
46 
extF80M_sqrt(const extFloat80_t * aPtr,extFloat80_t * zPtr)47 void extF80M_sqrt( const extFloat80_t *aPtr, extFloat80_t *zPtr )
48 {
49 
50     *zPtr = extF80_sqrt( *aPtr );
51 
52 }
53 
54 #else
55 
extF80M_sqrt(const extFloat80_t * aPtr,extFloat80_t * zPtr)56 void extF80M_sqrt( const extFloat80_t *aPtr, extFloat80_t *zPtr )
57 {
58     const struct extFloat80M *aSPtr;
59     struct extFloat80M *zSPtr;
60     uint_fast16_t uiA64, signUI64;
61     int32_t expA;
62     uint64_t rem64;
63     int32_t expZ;
64     uint32_t rem[4], sig32A, recipSqrt32, sig32Z, q;
65     uint64_t sig64Z, x64;
66     uint32_t term[4], extSigZ[3];
67 
68     /*------------------------------------------------------------------------
69     *------------------------------------------------------------------------*/
70     aSPtr = (const struct extFloat80M *) aPtr;
71     zSPtr = (struct extFloat80M *) zPtr;
72     /*------------------------------------------------------------------------
73     *------------------------------------------------------------------------*/
74     uiA64 = aSPtr->signExp;
75     signUI64 = uiA64 & packToExtF80UI64( 1, 0 );
76     expA = expExtF80UI64( uiA64 );
77     rem64 = aSPtr->signif;
78     /*------------------------------------------------------------------------
79     *------------------------------------------------------------------------*/
80     if ( expA == 0x7FFF ) {
81         if ( rem64 & UINT64_C( 0x7FFFFFFFFFFFFFFF ) ) {
82             softfloat_propagateNaNExtF80M( aSPtr, 0, zSPtr );
83             return;
84         }
85         if ( signUI64 ) goto invalid;
86         rem64 = UINT64_C( 0x8000000000000000 );
87         goto copyA;
88     }
89     /*------------------------------------------------------------------------
90     *------------------------------------------------------------------------*/
91     if ( ! expA ) expA = 1;
92     if ( ! (rem64 & UINT64_C( 0x8000000000000000 )) ) {
93         if ( ! rem64 ) {
94             uiA64 = signUI64;
95             goto copyA;
96         }
97         expA += softfloat_normExtF80SigM( &rem64 );
98     }
99     if ( signUI64 ) goto invalid;
100     /*------------------------------------------------------------------------
101     *------------------------------------------------------------------------*/
102     expZ = ((expA - 0x3FFF)>>1) + 0x3FFF;
103     expA &= 1;
104     softfloat_shortShiftLeft64To96M(
105         rem64, 30 - expA, &rem[indexMultiwordHi( 4, 3 )] );
106     sig32A = rem64>>32;
107     recipSqrt32 = softfloat_approxRecipSqrt32_1( expA, sig32A );
108     sig32Z = ((uint64_t) sig32A * recipSqrt32)>>32;
109     if ( expA ) sig32Z >>= 1;
110     rem64 =
111         ((uint64_t) rem[indexWord( 4, 3 )]<<32 | rem[indexWord( 4, 2 )])
112             - (uint64_t) sig32Z * sig32Z;
113     rem[indexWord( 4, 3 )] = rem64>>32;
114     rem[indexWord( 4, 2 )] = rem64;
115     /*------------------------------------------------------------------------
116     *------------------------------------------------------------------------*/
117     q = ((uint32_t) (rem64>>2) * (uint64_t) recipSqrt32)>>32;
118     sig64Z = ((uint64_t) sig32Z<<32) + ((uint64_t) q<<3);
119     x64 = ((uint64_t) sig32Z<<32) + sig64Z;
120     term[indexWord( 3, 2 )] = 0;
121     term[indexWord( 3, 1 )] = x64>>32;
122     term[indexWord( 3, 0 )] = x64;
123     softfloat_remStep96MBy32(
124         &rem[indexMultiwordHi( 4, 3 )],
125         29,
126         term,
127         q,
128         &rem[indexMultiwordHi( 4, 3 )]
129     );
130     rem64 = (uint64_t) rem[indexWord( 4, 3 )]<<32 | rem[indexWord( 4, 2 )];
131     /*------------------------------------------------------------------------
132     *------------------------------------------------------------------------*/
133     q = (((uint32_t) (rem64>>2) * (uint64_t) recipSqrt32)>>32) + 2;
134     x64 = (uint64_t) q<<7;
135     extSigZ[indexWord( 3, 0 )] = x64;
136     x64 = (sig64Z<<1) + (x64>>32);
137     extSigZ[indexWord( 3, 2 )] = x64>>32;
138     extSigZ[indexWord( 3, 1 )] = x64;
139     /*------------------------------------------------------------------------
140     *------------------------------------------------------------------------*/
141     if ( (q & 0xFFFFFF) <= 2 ) {
142         q &= ~(uint32_t) 0xFFFF;
143         extSigZ[indexWordLo( 3 )] = q<<7;
144         x64 = sig64Z + (q>>27);
145         term[indexWord( 4, 3 )] = 0;
146         term[indexWord( 4, 2 )] = x64>>32;
147         term[indexWord( 4, 1 )] = x64;
148         term[indexWord( 4, 0 )] = q<<5;
149         rem[indexWord( 4, 0 )] = 0;
150         softfloat_remStep128MBy32( rem, 28, term, q, rem );
151         q = rem[indexWordHi( 4 )];
152         if ( q & 0x80000000 ) {
153             softfloat_sub1X96M( extSigZ );
154         } else {
155             if ( q || rem[indexWord( 4, 1 )] || rem[indexWord( 4, 2 )] ) {
156                 extSigZ[indexWordLo( 3 )] |= 1;
157             }
158         }
159     }
160     softfloat_roundPackMToExtF80M(
161         0, expZ, extSigZ, extF80_roundingPrecision, zSPtr );
162     return;
163     /*------------------------------------------------------------------------
164     *------------------------------------------------------------------------*/
165  invalid:
166     softfloat_invalidExtF80M( zSPtr );
167     return;
168     /*------------------------------------------------------------------------
169     *------------------------------------------------------------------------*/
170  copyA:
171     zSPtr->signExp = uiA64;
172     zSPtr->signif  = rem64;
173 
174 }
175 
176 #endif
177 
178