1
2 /*============================================================================
3
4 This C source fragment is part of the SoftFloat IEC/IEEE Floating-point
5 Arithmetic Package, Release 2b.
6
7 Written by John R. Hauser. This work was made possible in part by the
8 International Computer Science Institute, located at Suite 600, 1947 Center
9 Street, Berkeley, California 94704. Funding was partially provided by the
10 National Science Foundation under grant MIP-9311980. The original version
11 of this code was written as part of a project to build a fixed-point vector
12 processor in collaboration with the University of California at Berkeley,
13 overseen by Profs. Nelson Morgan and John Wawrzynek. More information
14 is available through the Web page `http://www.cs.berkeley.edu/~jhauser/
15 arithmetic/SoftFloat.html'.
16
17 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has
18 been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
19 RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
20 AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
21 COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
22 EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
23 INSTITUTE (possibly via similar legal notice) AGAINST ALL LOSSES, COSTS, OR
24 OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
25
26 Derivative works are acceptable, even for commercial purposes, so long as
27 (1) the source code for the derivative work includes prominent notice that
28 the work is derivative, and (2) the source code includes prominent notice with
29 these four paragraphs for those parts of this code that are retained.
30
31 =============================================================================*/
32
33 /*----------------------------------------------------------------------------
34 | Shifts `a' right by the number of bits given in `count'. If any nonzero
35 | bits are shifted off, they are ``jammed'' into the least significant bit of
36 | the result by setting the least significant bit to 1. The value of `count'
37 | can be arbitrarily large; in particular, if `count' is greater than 32, the
38 | result will be either 0 or 1, depending on whether `a' is zero or nonzero.
39 | The result is stored in the location pointed to by `zPtr'.
40 *----------------------------------------------------------------------------*/
41
shift32RightJamming(bits32 a,int16 count,bits32 * zPtr)42 INLINE void shift32RightJamming( bits32 a, int16 count, bits32 *zPtr )
43 {
44 bits32 z;
45
46 if ( count == 0 ) {
47 z = a;
48 }
49 else if ( count < 32 ) {
50 z = ( a>>count ) | ( ( a<<( ( - count ) & 31 ) ) != 0 );
51 }
52 else {
53 z = ( a != 0 );
54 }
55 *zPtr = z;
56
57 }
58
59 /*----------------------------------------------------------------------------
60 | Shifts `a' right by the number of bits given in `count'. If any nonzero
61 | bits are shifted off, they are ``jammed'' into the least significant bit of
62 | the result by setting the least significant bit to 1. The value of `count'
63 | can be arbitrarily large; in particular, if `count' is greater than 64, the
64 | result will be either 0 or 1, depending on whether `a' is zero or nonzero.
65 | The result is stored in the location pointed to by `zPtr'.
66 *----------------------------------------------------------------------------*/
67
shift64RightJamming(bits64 a,int16 count,bits64 * zPtr)68 INLINE void shift64RightJamming( bits64 a, int16 count, bits64 *zPtr )
69 {
70 bits64 z;
71
72 if ( count == 0 ) {
73 z = a;
74 }
75 else if ( count < 64 ) {
76 z = ( a>>count ) | ( ( a<<( ( - count ) & 63 ) ) != 0 );
77 }
78 else {
79 z = ( a != 0 );
80 }
81 *zPtr = z;
82
83 }
84
85 /*----------------------------------------------------------------------------
86 | Shifts the 128-bit value formed by concatenating `a0' and `a1' right by 64
87 | _plus_ the number of bits given in `count'. The shifted result is at most
88 | 64 nonzero bits; this is stored at the location pointed to by `z0Ptr'. The
89 | bits shifted off form a second 64-bit result as follows: The _last_ bit
90 | shifted off is the most-significant bit of the extra result, and the other
91 | 63 bits of the extra result are all zero if and only if _all_but_the_last_
92 | bits shifted off were all zero. This extra result is stored in the location
93 | pointed to by `z1Ptr'. The value of `count' can be arbitrarily large.
94 | (This routine makes more sense if `a0' and `a1' are considered to form
95 | a fixed-point value with binary point between `a0' and `a1'. This fixed-
96 | point value is shifted right by the number of bits given in `count', and
97 | the integer part of the result is returned at the location pointed to by
98 | `z0Ptr'. The fractional part of the result may be slightly corrupted as
99 | described above, and is returned at the location pointed to by `z1Ptr'.)
100 *----------------------------------------------------------------------------*/
101
102 INLINE void
shift64ExtraRightJamming(bits64 a0,bits64 a1,int16 count,bits64 * z0Ptr,bits64 * z1Ptr)103 shift64ExtraRightJamming(
104 bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
105 {
106 bits64 z0, z1;
107 int8 negCount = ( - count ) & 63;
108
109 if ( count == 0 ) {
110 z1 = a1;
111 z0 = a0;
112 }
113 else if ( count < 64 ) {
114 z1 = ( a0<<negCount ) | ( a1 != 0 );
115 z0 = a0>>count;
116 }
117 else {
118 if ( count == 64 ) {
119 z1 = a0 | ( a1 != 0 );
120 }
121 else {
122 z1 = ( ( a0 | a1 ) != 0 );
123 }
124 z0 = 0;
125 }
126 *z1Ptr = z1;
127 *z0Ptr = z0;
128
129 }
130
131 /*----------------------------------------------------------------------------
132 | Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the
133 | number of bits given in `count'. Any bits shifted off are lost. The value
134 | of `count' can be arbitrarily large; in particular, if `count' is greater
135 | than 128, the result will be 0. The result is broken into two 64-bit pieces
136 | which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
137 *----------------------------------------------------------------------------*/
138
139 INLINE void
shift128Right(bits64 a0,bits64 a1,int16 count,bits64 * z0Ptr,bits64 * z1Ptr)140 shift128Right(
141 bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
142 {
143 bits64 z0, z1;
144 int8 negCount = ( - count ) & 63;
145
146 if ( count == 0 ) {
147 z1 = a1;
148 z0 = a0;
149 }
150 else if ( count < 64 ) {
151 z1 = ( a0<<negCount ) | ( a1>>count );
152 z0 = a0>>count;
153 }
154 else {
155 z1 = ( count < 64 ) ? ( a0>>( count & 63 ) ) : 0;
156 z0 = 0;
157 }
158 *z1Ptr = z1;
159 *z0Ptr = z0;
160
161 }
162
163 /*----------------------------------------------------------------------------
164 | Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the
165 | number of bits given in `count'. If any nonzero bits are shifted off, they
166 | are ``jammed'' into the least significant bit of the result by setting the
167 | least significant bit to 1. The value of `count' can be arbitrarily large;
168 | in particular, if `count' is greater than 128, the result will be either
169 | 0 or 1, depending on whether the concatenation of `a0' and `a1' is zero or
170 | nonzero. The result is broken into two 64-bit pieces which are stored at
171 | the locations pointed to by `z0Ptr' and `z1Ptr'.
172 *----------------------------------------------------------------------------*/
173
174 INLINE void
shift128RightJamming(bits64 a0,bits64 a1,int16 count,bits64 * z0Ptr,bits64 * z1Ptr)175 shift128RightJamming(
176 bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
177 {
178 bits64 z0, z1;
179 int8 negCount = ( - count ) & 63;
180
181 if ( count == 0 ) {
182 z1 = a1;
183 z0 = a0;
184 }
185 else if ( count < 64 ) {
186 z1 = ( a0<<negCount ) | ( a1>>count ) | ( ( a1<<negCount ) != 0 );
187 z0 = a0>>count;
188 }
189 else {
190 if ( count == 64 ) {
191 z1 = a0 | ( a1 != 0 );
192 }
193 else if ( count < 128 ) {
194 z1 = ( a0>>( count & 63 ) ) | ( ( ( a0<<negCount ) | a1 ) != 0 );
195 }
196 else {
197 z1 = ( ( a0 | a1 ) != 0 );
198 }
199 z0 = 0;
200 }
201 *z1Ptr = z1;
202 *z0Ptr = z0;
203
204 }
205
206 /*----------------------------------------------------------------------------
207 | Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' right
208 | by 64 _plus_ the number of bits given in `count'. The shifted result is
209 | at most 128 nonzero bits; these are broken into two 64-bit pieces which are
210 | stored at the locations pointed to by `z0Ptr' and `z1Ptr'. The bits shifted
211 | off form a third 64-bit result as follows: The _last_ bit shifted off is
212 | the most-significant bit of the extra result, and the other 63 bits of the
213 | extra result are all zero if and only if _all_but_the_last_ bits shifted off
214 | were all zero. This extra result is stored in the location pointed to by
215 | `z2Ptr'. The value of `count' can be arbitrarily large.
216 | (This routine makes more sense if `a0', `a1', and `a2' are considered
217 | to form a fixed-point value with binary point between `a1' and `a2'. This
218 | fixed-point value is shifted right by the number of bits given in `count',
219 | and the integer part of the result is returned at the locations pointed to
220 | by `z0Ptr' and `z1Ptr'. The fractional part of the result may be slightly
221 | corrupted as described above, and is returned at the location pointed to by
222 | `z2Ptr'.)
223 *----------------------------------------------------------------------------*/
224
225 INLINE void
shift128ExtraRightJamming(bits64 a0,bits64 a1,bits64 a2,int16 count,bits64 * z0Ptr,bits64 * z1Ptr,bits64 * z2Ptr)226 shift128ExtraRightJamming(
227 bits64 a0,
228 bits64 a1,
229 bits64 a2,
230 int16 count,
231 bits64 *z0Ptr,
232 bits64 *z1Ptr,
233 bits64 *z2Ptr
234 )
235 {
236 bits64 z0, z1, z2;
237 int8 negCount = ( - count ) & 63;
238
239 if ( count == 0 ) {
240 z2 = a2;
241 z1 = a1;
242 z0 = a0;
243 }
244 else {
245 if ( count < 64 ) {
246 z2 = a1<<negCount;
247 z1 = ( a0<<negCount ) | ( a1>>count );
248 z0 = a0>>count;
249 }
250 else {
251 if ( count == 64 ) {
252 z2 = a1;
253 z1 = a0;
254 }
255 else {
256 a2 |= a1;
257 if ( count < 128 ) {
258 z2 = a0<<negCount;
259 z1 = a0>>( count & 63 );
260 }
261 else {
262 z2 = ( count == 128 ) ? a0 : ( a0 != 0 );
263 z1 = 0;
264 }
265 }
266 z0 = 0;
267 }
268 z2 |= ( a2 != 0 );
269 }
270 *z2Ptr = z2;
271 *z1Ptr = z1;
272 *z0Ptr = z0;
273
274 }
275
276 /*----------------------------------------------------------------------------
277 | Shifts the 128-bit value formed by concatenating `a0' and `a1' left by the
278 | number of bits given in `count'. Any bits shifted off are lost. The value
279 | of `count' must be less than 64. The result is broken into two 64-bit
280 | pieces which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
281 *----------------------------------------------------------------------------*/
282
283 INLINE void
shortShift128Left(bits64 a0,bits64 a1,int16 count,bits64 * z0Ptr,bits64 * z1Ptr)284 shortShift128Left(
285 bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
286 {
287
288 *z1Ptr = a1<<count;
289 *z0Ptr =
290 ( count == 0 ) ? a0 : ( a0<<count ) | ( a1>>( ( - count ) & 63 ) );
291
292 }
293
294 /*----------------------------------------------------------------------------
295 | Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' left
296 | by the number of bits given in `count'. Any bits shifted off are lost.
297 | The value of `count' must be less than 64. The result is broken into three
298 | 64-bit pieces which are stored at the locations pointed to by `z0Ptr',
299 | `z1Ptr', and `z2Ptr'.
300 *----------------------------------------------------------------------------*/
301
302 INLINE void
shortShift192Left(bits64 a0,bits64 a1,bits64 a2,int16 count,bits64 * z0Ptr,bits64 * z1Ptr,bits64 * z2Ptr)303 shortShift192Left(
304 bits64 a0,
305 bits64 a1,
306 bits64 a2,
307 int16 count,
308 bits64 *z0Ptr,
309 bits64 *z1Ptr,
310 bits64 *z2Ptr
311 )
312 {
313 bits64 z0, z1, z2;
314 int8 negCount;
315
316 z2 = a2<<count;
317 z1 = a1<<count;
318 z0 = a0<<count;
319 if ( 0 < count ) {
320 negCount = ( ( - count ) & 63 );
321 z1 |= a2>>negCount;
322 z0 |= a1>>negCount;
323 }
324 *z2Ptr = z2;
325 *z1Ptr = z1;
326 *z0Ptr = z0;
327
328 }
329
330 /*----------------------------------------------------------------------------
331 | Adds the 128-bit value formed by concatenating `a0' and `a1' to the 128-bit
332 | value formed by concatenating `b0' and `b1'. Addition is modulo 2^128, so
333 | any carry out is lost. The result is broken into two 64-bit pieces which
334 | are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
335 *----------------------------------------------------------------------------*/
336
337 INLINE void
add128(bits64 a0,bits64 a1,bits64 b0,bits64 b1,bits64 * z0Ptr,bits64 * z1Ptr)338 add128(
339 bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 *z0Ptr, bits64 *z1Ptr )
340 {
341 bits64 z1;
342
343 z1 = a1 + b1;
344 *z1Ptr = z1;
345 *z0Ptr = a0 + b0 + ( z1 < a1 );
346
347 }
348
349 /*----------------------------------------------------------------------------
350 | Adds the 192-bit value formed by concatenating `a0', `a1', and `a2' to the
351 | 192-bit value formed by concatenating `b0', `b1', and `b2'. Addition is
352 | modulo 2^192, so any carry out is lost. The result is broken into three
353 | 64-bit pieces which are stored at the locations pointed to by `z0Ptr',
354 | `z1Ptr', and `z2Ptr'.
355 *----------------------------------------------------------------------------*/
356
357 INLINE void
add192(bits64 a0,bits64 a1,bits64 a2,bits64 b0,bits64 b1,bits64 b2,bits64 * z0Ptr,bits64 * z1Ptr,bits64 * z2Ptr)358 add192(
359 bits64 a0,
360 bits64 a1,
361 bits64 a2,
362 bits64 b0,
363 bits64 b1,
364 bits64 b2,
365 bits64 *z0Ptr,
366 bits64 *z1Ptr,
367 bits64 *z2Ptr
368 )
369 {
370 bits64 z0, z1, z2;
371 int8 carry0, carry1;
372
373 z2 = a2 + b2;
374 carry1 = ( z2 < a2 );
375 z1 = a1 + b1;
376 carry0 = ( z1 < a1 );
377 z0 = a0 + b0;
378 z1 += carry1;
379 z0 += ( z1 < carry1 );
380 z0 += carry0;
381 *z2Ptr = z2;
382 *z1Ptr = z1;
383 *z0Ptr = z0;
384
385 }
386
387 /*----------------------------------------------------------------------------
388 | Subtracts the 128-bit value formed by concatenating `b0' and `b1' from the
389 | 128-bit value formed by concatenating `a0' and `a1'. Subtraction is modulo
390 | 2^128, so any borrow out (carry out) is lost. The result is broken into two
391 | 64-bit pieces which are stored at the locations pointed to by `z0Ptr' and
392 | `z1Ptr'.
393 *----------------------------------------------------------------------------*/
394
395 INLINE void
sub128(bits64 a0,bits64 a1,bits64 b0,bits64 b1,bits64 * z0Ptr,bits64 * z1Ptr)396 sub128(
397 bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 *z0Ptr, bits64 *z1Ptr )
398 {
399
400 *z1Ptr = a1 - b1;
401 *z0Ptr = a0 - b0 - ( a1 < b1 );
402
403 }
404
405 /*----------------------------------------------------------------------------
406 | Subtracts the 192-bit value formed by concatenating `b0', `b1', and `b2'
407 | from the 192-bit value formed by concatenating `a0', `a1', and `a2'.
408 | Subtraction is modulo 2^192, so any borrow out (carry out) is lost. The
409 | result is broken into three 64-bit pieces which are stored at the locations
410 | pointed to by `z0Ptr', `z1Ptr', and `z2Ptr'.
411 *----------------------------------------------------------------------------*/
412
413 INLINE void
sub192(bits64 a0,bits64 a1,bits64 a2,bits64 b0,bits64 b1,bits64 b2,bits64 * z0Ptr,bits64 * z1Ptr,bits64 * z2Ptr)414 sub192(
415 bits64 a0,
416 bits64 a1,
417 bits64 a2,
418 bits64 b0,
419 bits64 b1,
420 bits64 b2,
421 bits64 *z0Ptr,
422 bits64 *z1Ptr,
423 bits64 *z2Ptr
424 )
425 {
426 bits64 z0, z1, z2;
427 int8 borrow0, borrow1;
428
429 z2 = a2 - b2;
430 borrow1 = ( a2 < b2 );
431 z1 = a1 - b1;
432 borrow0 = ( a1 < b1 );
433 z0 = a0 - b0;
434 z0 -= ( z1 < borrow1 );
435 z1 -= borrow1;
436 z0 -= borrow0;
437 *z2Ptr = z2;
438 *z1Ptr = z1;
439 *z0Ptr = z0;
440
441 }
442
443 /*----------------------------------------------------------------------------
444 | Multiplies `a' by `b' to obtain a 128-bit product. The product is broken
445 | into two 64-bit pieces which are stored at the locations pointed to by
446 | `z0Ptr' and `z1Ptr'.
447 *----------------------------------------------------------------------------*/
448
mul64To128(bits64 a,bits64 b,bits64 * z0Ptr,bits64 * z1Ptr)449 INLINE void mul64To128( bits64 a, bits64 b, bits64 *z0Ptr, bits64 *z1Ptr )
450 {
451 bits32 aHigh, aLow, bHigh, bLow;
452 bits64 z0, zMiddleA, zMiddleB, z1;
453
454 aLow = a;
455 aHigh = a>>32;
456 bLow = b;
457 bHigh = b>>32;
458 z1 = ( (bits64) aLow ) * bLow;
459 zMiddleA = ( (bits64) aLow ) * bHigh;
460 zMiddleB = ( (bits64) aHigh ) * bLow;
461 z0 = ( (bits64) aHigh ) * bHigh;
462 zMiddleA += zMiddleB;
463 z0 += ( ( (bits64) ( zMiddleA < zMiddleB ) )<<32 ) + ( zMiddleA>>32 );
464 zMiddleA <<= 32;
465 z1 += zMiddleA;
466 z0 += ( z1 < zMiddleA );
467 *z1Ptr = z1;
468 *z0Ptr = z0;
469
470 }
471
472 /*----------------------------------------------------------------------------
473 | Multiplies the 128-bit value formed by concatenating `a0' and `a1' by
474 | `b' to obtain a 192-bit product. The product is broken into three 64-bit
475 | pieces which are stored at the locations pointed to by `z0Ptr', `z1Ptr', and
476 | `z2Ptr'.
477 *----------------------------------------------------------------------------*/
478
479 INLINE void
mul128By64To192(bits64 a0,bits64 a1,bits64 b,bits64 * z0Ptr,bits64 * z1Ptr,bits64 * z2Ptr)480 mul128By64To192(
481 bits64 a0,
482 bits64 a1,
483 bits64 b,
484 bits64 *z0Ptr,
485 bits64 *z1Ptr,
486 bits64 *z2Ptr
487 )
488 {
489 bits64 z0, z1, z2, more1;
490
491 mul64To128( a1, b, &z1, &z2 );
492 mul64To128( a0, b, &z0, &more1 );
493 add128( z0, more1, 0, z1, &z0, &z1 );
494 *z2Ptr = z2;
495 *z1Ptr = z1;
496 *z0Ptr = z0;
497
498 }
499
500 /*----------------------------------------------------------------------------
501 | Multiplies the 128-bit value formed by concatenating `a0' and `a1' to the
502 | 128-bit value formed by concatenating `b0' and `b1' to obtain a 256-bit
503 | product. The product is broken into four 64-bit pieces which are stored at
504 | the locations pointed to by `z0Ptr', `z1Ptr', `z2Ptr', and `z3Ptr'.
505 *----------------------------------------------------------------------------*/
506
507 INLINE void
mul128To256(bits64 a0,bits64 a1,bits64 b0,bits64 b1,bits64 * z0Ptr,bits64 * z1Ptr,bits64 * z2Ptr,bits64 * z3Ptr)508 mul128To256(
509 bits64 a0,
510 bits64 a1,
511 bits64 b0,
512 bits64 b1,
513 bits64 *z0Ptr,
514 bits64 *z1Ptr,
515 bits64 *z2Ptr,
516 bits64 *z3Ptr
517 )
518 {
519 bits64 z0, z1, z2, z3;
520 bits64 more1, more2;
521
522 mul64To128( a1, b1, &z2, &z3 );
523 mul64To128( a1, b0, &z1, &more2 );
524 add128( z1, more2, 0, z2, &z1, &z2 );
525 mul64To128( a0, b0, &z0, &more1 );
526 add128( z0, more1, 0, z1, &z0, &z1 );
527 mul64To128( a0, b1, &more1, &more2 );
528 add128( more1, more2, 0, z2, &more1, &z2 );
529 add128( z0, z1, 0, more1, &z0, &z1 );
530 *z3Ptr = z3;
531 *z2Ptr = z2;
532 *z1Ptr = z1;
533 *z0Ptr = z0;
534
535 }
536
537 /*----------------------------------------------------------------------------
538 | Returns an approximation to the 64-bit integer quotient obtained by dividing
539 | `b' into the 128-bit value formed by concatenating `a0' and `a1'. The
540 | divisor `b' must be at least 2^63. If q is the exact quotient truncated
541 | toward zero, the approximation returned lies between q and q + 2 inclusive.
542 | If the exact quotient q is larger than 64 bits, the maximum positive 64-bit
543 | unsigned integer is returned.
544 *----------------------------------------------------------------------------*/
545
estimateDiv128To64(bits64 a0,bits64 a1,bits64 b)546 static bits64 estimateDiv128To64( bits64 a0, bits64 a1, bits64 b )
547 {
548 bits64 b0, b1;
549 bits64 rem0, rem1, term0, term1;
550 bits64 z;
551
552 if ( b <= a0 ) return LIT64( 0xFFFFFFFFFFFFFFFF );
553 b0 = b>>32;
554 z = ( b0<<32 <= a0 ) ? LIT64( 0xFFFFFFFF00000000 ) : ( a0 / b0 )<<32;
555 mul64To128( b, z, &term0, &term1 );
556 sub128( a0, a1, term0, term1, &rem0, &rem1 );
557 while ( ( (sbits64) rem0 ) < 0 ) {
558 z -= LIT64( 0x100000000 );
559 b1 = b<<32;
560 add128( rem0, rem1, b0, b1, &rem0, &rem1 );
561 }
562 rem0 = ( rem0<<32 ) | ( rem1>>32 );
563 z |= ( b0<<32 <= rem0 ) ? 0xFFFFFFFF : rem0 / b0;
564 return z;
565
566 }
567
568 /*----------------------------------------------------------------------------
569 | Returns an approximation to the square root of the 32-bit significand given
570 | by `a'. Considered as an integer, `a' must be at least 2^31. If bit 0 of
571 | `aExp' (the least significant bit) is 1, the integer returned approximates
572 | 2^31*sqrt(`a'/2^31), where `a' is considered an integer. If bit 0 of `aExp'
573 | is 0, the integer returned approximates 2^31*sqrt(`a'/2^30). In either
574 | case, the approximation returned lies strictly within +/-2 of the exact
575 | value.
576 *----------------------------------------------------------------------------*/
577
estimateSqrt32(int16 aExp,bits32 a)578 static bits32 estimateSqrt32( int16 aExp, bits32 a )
579 {
580 static const bits16 sqrtOddAdjustments[] = {
581 0x0004, 0x0022, 0x005D, 0x00B1, 0x011D, 0x019F, 0x0236, 0x02E0,
582 0x039C, 0x0468, 0x0545, 0x0631, 0x072B, 0x0832, 0x0946, 0x0A67
583 };
584 static const bits16 sqrtEvenAdjustments[] = {
585 0x0A2D, 0x08AF, 0x075A, 0x0629, 0x051A, 0x0429, 0x0356, 0x029E,
586 0x0200, 0x0179, 0x0109, 0x00AF, 0x0068, 0x0034, 0x0012, 0x0002
587 };
588 int8 index;
589 bits32 z;
590
591 index = ( a>>27 ) & 15;
592 if ( aExp & 1 ) {
593 z = 0x4000 + ( a>>17 ) - sqrtOddAdjustments[ index ];
594 z = ( ( a / z )<<14 ) + ( z<<15 );
595 a >>= 1;
596 }
597 else {
598 z = 0x8000 + ( a>>17 ) - sqrtEvenAdjustments[ index ];
599 z = a / z + z;
600 z = ( 0x20000 <= z ) ? 0xFFFF8000 : ( z<<15 );
601 if ( z <= a ) return (bits32) ( ( (sbits32) a )>>1 );
602 }
603 return ( (bits32) ( ( ( (bits64) a )<<31 ) / z ) ) + ( z>>1 );
604
605 }
606
607 /*----------------------------------------------------------------------------
608 | Returns the number of leading 0 bits before the most-significant 1 bit of
609 | `a'. If `a' is zero, 32 is returned.
610 *----------------------------------------------------------------------------*/
611
countLeadingZeros32(bits32 a)612 static int8 countLeadingZeros32( bits32 a )
613 {
614 static const int8 countLeadingZerosHigh[] = {
615 8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
616 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
617 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
618 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
619 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
620 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
621 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
622 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
623 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
624 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
625 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
626 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
627 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
628 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
629 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
630 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
631 };
632 int8 shiftCount;
633
634 shiftCount = 0;
635 if ( a < 0x10000 ) {
636 shiftCount += 16;
637 a <<= 16;
638 }
639 if ( a < 0x1000000 ) {
640 shiftCount += 8;
641 a <<= 8;
642 }
643 shiftCount += countLeadingZerosHigh[ a>>24 ];
644 return shiftCount;
645
646 }
647
648 /*----------------------------------------------------------------------------
649 | Returns the number of leading 0 bits before the most-significant 1 bit of
650 | `a'. If `a' is zero, 64 is returned.
651 *----------------------------------------------------------------------------*/
652
countLeadingZeros64(bits64 a)653 static int8 countLeadingZeros64( bits64 a )
654 {
655 int8 shiftCount;
656
657 shiftCount = 0;
658 if ( a < ( (bits64) 1 )<<32 ) {
659 shiftCount += 32;
660 }
661 else {
662 a >>= 32;
663 }
664 shiftCount += countLeadingZeros32( a );
665 return shiftCount;
666
667 }
668
669 /*----------------------------------------------------------------------------
670 | Returns 1 if the 128-bit value formed by concatenating `a0' and `a1'
671 | is equal to the 128-bit value formed by concatenating `b0' and `b1'.
672 | Otherwise, returns 0.
673 *----------------------------------------------------------------------------*/
674
eq128(bits64 a0,bits64 a1,bits64 b0,bits64 b1)675 INLINE flag eq128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
676 {
677
678 return ( a0 == b0 ) && ( a1 == b1 );
679
680 }
681
682 /*----------------------------------------------------------------------------
683 | Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is less
684 | than or equal to the 128-bit value formed by concatenating `b0' and `b1'.
685 | Otherwise, returns 0.
686 *----------------------------------------------------------------------------*/
687
le128(bits64 a0,bits64 a1,bits64 b0,bits64 b1)688 INLINE flag le128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
689 {
690
691 return ( a0 < b0 ) || ( ( a0 == b0 ) && ( a1 <= b1 ) );
692
693 }
694
695 /*----------------------------------------------------------------------------
696 | Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is less
697 | than the 128-bit value formed by concatenating `b0' and `b1'. Otherwise,
698 | returns 0.
699 *----------------------------------------------------------------------------*/
700
lt128(bits64 a0,bits64 a1,bits64 b0,bits64 b1)701 INLINE flag lt128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
702 {
703
704 return ( a0 < b0 ) || ( ( a0 == b0 ) && ( a1 < b1 ) );
705
706 }
707
708 /*----------------------------------------------------------------------------
709 | Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is
710 | not equal to the 128-bit value formed by concatenating `b0' and `b1'.
711 | Otherwise, returns 0.
712 *----------------------------------------------------------------------------*/
713
ne128(bits64 a0,bits64 a1,bits64 b0,bits64 b1)714 INLINE flag ne128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
715 {
716
717 return ( a0 != b0 ) || ( a1 != b1 );
718
719 }
720
721