Bug Summary

File:dconv/dragon4.cpp
Location:line 702, column 22
Description:Assigned value is garbage or undefined

Annotated Source Code

1/******************************************************************************
2 Copyright (c) 2014 Ryan Juckett
3 http://www.ryanjuckett.com/
4
5 This software is provided 'as-is', without any express or implied
6 warranty. In no event will the authors be held liable for any damages
7 arising from the use of this software.
8
9 Permission is granted to anyone to use this software for any purpose,
10 including commercial applications, and to alter it and redistribute it
11 freely, subject to the following restrictions:
12
13 1. The origin of this software must not be misrepresented; you must not
14 claim that you wrote the original software. If you use this software
15 in a product, an acknowledgment in the product documentation would be
16 appreciated but is not required.
17
18 2. Altered source versions must be plainly marked as such, and must not be
19 misrepresented as being the original software.
20
21 3. This notice may not be removed or altered from any source
22 distribution.
23*******************************************************************************
24 Copyright (c) 2015-2017 Carsten Sonne Larsen <cs@innolan.dk>
25 All rights reserved.
26
27 Redistribution and use in source and binary forms, with or without
28 modification, are permitted provided that the following conditions
29 are met:
30 1. Redistributions of source code must retain the above copyright
31 notice, this list of conditions and the following disclaimer.
32
33 2. Redistributions in binary form must reproduce the above copyright
34 notice, this list of conditions and the following disclaimer in the
35 documentation and/or other materials provided with the distribution.
36
37 THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
38 IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
39 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
40 IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
41 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
42 NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
43 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
45 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
46 THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
47
48 The origin source code can be obtained from:
49 http://www.ryanjuckett.com/
50
51******************************************************************************/
52
53#include "math.h"
54#include "dmath.h"
55#include "dragon4.h"
56
57//******************************************************************************
58// Maximum number of 32 bit blocks needed in high precision arithmetic
59// to print out 64 bit IEEE floating point values.
60//******************************************************************************
61const tU32 c_BigInt_MaxBlocks = 35;
62
63//******************************************************************************
64// This structure stores a high precision unsigned integer. It uses a buffer
65// of 32 bit integer blocks along with a length. The lowest bits of the integer
66// are stored at the start of the buffer and the length is set to the minimum
67// value that contains the integer. Thus, there are never any zero blocks at the
68// end of the buffer.
69//******************************************************************************
70struct tBigInt
71{
72 // Copy integer
73 tBigInt & operator=(const tBigInt &rhs)
74 {
75 tU32 length = rhs.m_length;
76 tU32 * pLhsCur = m_blocks;
77 for (const tU32 *pRhsCur = rhs.m_blocks, *pRhsEnd = pRhsCur + length;
78 pRhsCur != pRhsEnd;
79 ++pLhsCur, ++pRhsCur)
80 {
81 *pLhsCur = *pRhsCur;
82 }
83 m_length = length;
84 return *this;
85 }
86
87 // Data accessors
88 tU32 GetLength() const {
89 return m_length;
90 }
91 tU32 GetBlock(tU32 idx) const {
92 return m_blocks[idx];
93 }
94
95 // Zero helper functions
96 void SetZero() {
97 m_length = 0;
98 }
99 tB IsZero() const {
100 return m_length == 0;
101 }
102
103 // Basic type accessors
104 void SetU64(tU64 val)
105 {
106 if (val > 0xFFFFFFFF)
107 {
108 m_blocks[0] = val & 0xFFFFFFFF;
109 m_blocks[1] = (val >> 32) & 0xFFFFFFFF;
110 m_length = 2;
111 }
112 else if (val != 0)
113 {
114 m_blocks[0] = val & 0xFFFFFFFF;
115 m_length = 1;
116 }
117 else
118 {
119 m_length = 0;
120 }
121 }
122
123 void SetU32(tU32 val)
124 {
125 if (val != 0)
126 {
127 m_blocks[0] = val;
128 m_length = (val != 0);
129 }
130 else
131 {
132 m_length = 0;
133 }
134 }
135
136 tU32 GetU32() const {
137 return (m_length == 0) ? 0 : m_blocks[0];
138 }
139
140 // Member data
141 tU32 m_length;
142 tU32 m_blocks[c_BigInt_MaxBlocks];
143};
144
145//******************************************************************************
146// Returns 0 if (lhs = rhs), negative if (lhs < rhs), positive if (lhs > rhs)
147//******************************************************************************
148static tS32 BigInt_Compare(const tBigInt & lhs, const tBigInt & rhs)
149{
150 // A bigger length implies a bigger number.
151 tS32 lengthDiff = lhs.m_length - rhs.m_length;
152 if (lengthDiff != 0)
153 return lengthDiff;
154
155 // Compare blocks one by one from high to low.
156 for (tS32 i = lhs.m_length - 1; i >= 0; --i)
157 {
158 if (lhs.m_blocks[i] == rhs.m_blocks[i])
159 continue;
160 else if (lhs.m_blocks[i] > rhs.m_blocks[i])
161 return 1;
162 else
163 return -1;
164 }
165
166 // no blocks differed
167 return 0;
168}
169
170//******************************************************************************
171// result = lhs + rhs
172//******************************************************************************
173static void BigInt_Add(tBigInt * pResult, const tBigInt & lhs, const tBigInt & rhs)
174{
175 // determine which operand has the smaller length
176 const tBigInt * pLarge;
177 const tBigInt * pSmall;
178 if (lhs.m_length < rhs.m_length)
179 {
180 pSmall = &lhs;
181 pLarge = &rhs;
182 }
183 else
184 {
185 pSmall = &rhs;
186 pLarge = &lhs;
187 }
188
189 const tU32 largeLen = pLarge->m_length;
190 const tU32 smallLen = pSmall->m_length;
191
192 // The output will be at least as long as the largest input
193 pResult->m_length = largeLen;
194
195 // Add each block and add carry the overflow to the next block
196 tU64 carry = 0;
197 const tU32 * pLargeCur = pLarge->m_blocks;
198 const tU32 * pLargeEnd = pLargeCur + largeLen;
199 const tU32 * pSmallCur = pSmall->m_blocks;
200 const tU32 * pSmallEnd = pSmallCur + smallLen;
201 tU32 * pResultCur = pResult->m_blocks;
202 while (pSmallCur != pSmallEnd)
203 {
204 tU64 sum = carry + (tU64)(*pLargeCur) + (tU64)(*pSmallCur);
205 carry = sum >> 32;
206 (*pResultCur) = sum & 0xFFFFFFFF;
207 ++pLargeCur;
208 ++pSmallCur;
209 ++pResultCur;
210 }
211
212 // Add the carry to any blocks that only exist in the large operand
213 while (pLargeCur != pLargeEnd)
214 {
215 tU64 sum = carry + (tU64)(*pLargeCur);
216 carry = sum >> 32;
217 (*pResultCur) = sum & 0xFFFFFFFF;
218 ++pLargeCur;
219 ++pResultCur;
220 }
221
222 // If there's still a carry, append a new block
223 if (carry != 0)
224 {
225 RJ_ASSERT(carry == 1);
226 RJ_ASSERT((tU32)(pResultCur - pResult->m_blocks) == largeLen && (largeLen < c_BigInt_MaxBlocks));
227 *pResultCur = 1;
228 pResult->m_length = largeLen + 1;
229 }
230 else
231 {
232 pResult->m_length = largeLen;
233 }
234}
235
236//******************************************************************************
237// result = lhs * rhs
238//******************************************************************************
239static void BigInt_Multiply(tBigInt * pResult, const tBigInt &lhs, const tBigInt &rhs)
240{
241 RJ_ASSERT( pResult != &lhs && pResult != &rhs );
242
243 // determine which operand has the smaller length
244 const tBigInt * pLarge;
245 const tBigInt * pSmall;
246 if (lhs.m_length < rhs.m_length)
247 {
248 pSmall = &lhs;
249 pLarge = &rhs;
250 }
251 else
252 {
253 pSmall = &rhs;
254 pLarge = &lhs;
255 }
256
257 // set the maximum possible result length
258 tU32 maxResultLen = pLarge->m_length + pSmall->m_length;
259 RJ_ASSERT( maxResultLen <= c_BigInt_MaxBlocks );
260
261 // clear the result data
262 for(tU32 * pCur = pResult->m_blocks, *pEnd = pCur + maxResultLen; pCur != pEnd; ++pCur)
263 *pCur = 0;
264
265 // perform standard long multiplication
266 const tU32 *pLargeBeg = pLarge->m_blocks;
267 const tU32 *pLargeEnd = pLargeBeg + pLarge->m_length;
268
269 // for each small block
270 tU32 *pResultStart = pResult->m_blocks;
271 for(const tU32 *pSmallCur = pSmall->m_blocks, *pSmallEnd = pSmallCur + pSmall->m_length;
272 pSmallCur != pSmallEnd;
273 ++pSmallCur, ++pResultStart)
274 {
275 // if non-zero, multiply against all the large blocks and add into the result
276 const tU32 multiplier = *pSmallCur;
277 if (multiplier != 0)
278 {
279 const tU32 *pLargeCur = pLargeBeg;
280 tU32 *pResultCur = pResultStart;
281 tU64 carry = 0;
282 do
283 {
284 tU64 product = (*pResultCur) + (*pLargeCur)*(tU64)multiplier + carry;
285 carry = product >> 32;
286 *pResultCur = product & 0xFFFFFFFF;
287 ++pLargeCur;
288 ++pResultCur;
289 } while(pLargeCur != pLargeEnd);
290
291 RJ_ASSERT(pResultCur < pResult->m_blocks + maxResultLen);
292 *pResultCur = (tU32)(carry & 0xFFFFFFFF);
293 }
294 }
295
296 // check if the terminating block has no set bits
297 if (maxResultLen > 0 && pResult->m_blocks[maxResultLen - 1] == 0)
298 pResult->m_length = maxResultLen-1;
299 else
300 pResult->m_length = maxResultLen;
301}
302
303//******************************************************************************
304// result = lhs * rhs
305//******************************************************************************
306static void BigInt_Multiply(tBigInt * pResult, const tBigInt & lhs, tU32 rhs)
307{
308 // perform long multiplication
309 tU32 carry = 0;
310 tU32 *pResultCur = pResult->m_blocks;
311 const tU32 *pLhsCur = lhs.m_blocks;
312 const tU32 *pLhsEnd = lhs.m_blocks + lhs.m_length;
313 for ( ; pLhsCur != pLhsEnd; ++pLhsCur, ++pResultCur )
314 {
315 tU64 product = (tU64)(*pLhsCur) * rhs + carry;
316 *pResultCur = (tU32)(product & 0xFFFFFFFF);
317 carry = product >> 32;
318 }
319
320 // if there is a remaining carry, grow the array
321 if (carry != 0)
322 {
323 // grow the array
324 RJ_ASSERT(lhs.m_length + 1 <= c_BigInt_MaxBlocks);
325 *pResultCur = (tU32)carry;
326 pResult->m_length = lhs.m_length + 1;
327 }
328 else
329 {
330 pResult->m_length = lhs.m_length;
331 }
332}
333
334//******************************************************************************
335// result = in * 2
336//******************************************************************************
337static void BigInt_Multiply2(tBigInt * pResult, const tBigInt &in)
338{
339 // shift all the blocks by one
340 tU32 carry = 0;
341
342 tU32 *pResultCur = pResult->m_blocks;
343 const tU32 *pLhsCur = in.m_blocks;
344 const tU32 *pLhsEnd = in.m_blocks + in.m_length;
345 for ( ; pLhsCur != pLhsEnd; ++pLhsCur, ++pResultCur )
346 {
347 tU32 cur = *pLhsCur;
348 *pResultCur = (cur << 1) | carry;
349 carry = cur >> 31;
350 }
351
352 if (carry != 0)
353 {
354 // grow the array
355 RJ_ASSERT(in.m_length + 1 <= c_BigInt_MaxBlocks);
356 *pResultCur = carry;
357 pResult->m_length = in.m_length + 1;
358 }
359 else
360 {
361 pResult->m_length = in.m_length;
362 }
363}
364
365//******************************************************************************
366// result = result * 2
367//******************************************************************************
368static void BigInt_Multiply2(tBigInt * pResult)
369{
370 // shift all the blocks by one
371 tU32 carry = 0;
372
373 tU32 *pCur = pResult->m_blocks;
374 tU32 *pEnd = pResult->m_blocks + pResult->m_length;
375 for ( ; pCur != pEnd; ++pCur )
376 {
377 tU32 cur = *pCur;
378 *pCur = (cur << 1) | carry;
379 carry = cur >> 31;
380 }
381
382 if (carry != 0)
383 {
384 // grow the array
385 RJ_ASSERT(pResult->m_length + 1 <= c_BigInt_MaxBlocks);
386 *pCur = carry;
387 ++pResult->m_length;
388 }
389}
390
391//******************************************************************************
392// result = result * 10
393//******************************************************************************
394static void BigInt_Multiply10(tBigInt * pResult)
395{
396 // multiply all the blocks
397 tU64 carry = 0;
398
399 tU32 *pCur = pResult->m_blocks;
400 tU32 *pEnd = pResult->m_blocks + pResult->m_length;
401 for ( ; pCur != pEnd; ++pCur )
402 {
403 tU64 product = (tU64)(*pCur) * 10ull + carry;
404 (*pCur) = (tU32)(product & 0xFFFFFFFF);
405 carry = product >> 32;
406 }
407
408 if (carry != 0)
409 {
410 // grow the array
411 RJ_ASSERT(pResult->m_length + 1 <= c_BigInt_MaxBlocks);
412 *pCur = (tU32)carry;
413 ++pResult->m_length;
414 }
415}
416
417//******************************************************************************
418//******************************************************************************
419static tU32 g_PowerOf10_U32[] =
420{
421 1, // 10 ^ 0
422 10, // 10 ^ 1
423 100, // 10 ^ 2
424 1000, // 10 ^ 3
425 10000, // 10 ^ 4
426 100000, // 10 ^ 5
427 1000000, // 10 ^ 6
428 10000000, // 10 ^ 7
429};
430
431//******************************************************************************
432// Note: This has a lot of wasted space in the big integer structures of the
433// early table entries. It wouldn't be terribly hard to make the multiply
434// function work on integer pointers with an array length instead of
435// the tBigInt struct which would allow us to store a minimal amount of
436// data here.
437//******************************************************************************
438static tBigInt g_PowerOf10_Big[] =
439{
440 // 10 ^ 8
441 { 1, { 100000000 } },
442 // 10 ^ 16
443 { 2, { 0x6fc10000, 0x002386f2 } },
444 // 10 ^ 32
445 { 4, { 0x00000000, 0x85acef81, 0x2d6d415b, 0x000004ee, } },
446 // 10 ^ 64
447 { 7, { 0x00000000, 0x00000000, 0xbf6a1f01, 0x6e38ed64, 0xdaa797ed, 0xe93ff9f4, 0x00184f03, } },
448 // 10 ^ 128
449 { 14, { 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x2e953e01, 0x03df9909, 0x0f1538fd,
450 0x2374e42f, 0xd3cff5ec, 0xc404dc08, 0xbccdb0da, 0xa6337f19, 0xe91f2603, 0x0000024e,
451 }
452 },
453 // 10 ^ 256
454 { 27, { 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
455 0x00000000, 0x982e7c01, 0xbed3875b, 0xd8d99f72, 0x12152f87, 0x6bde50c6, 0xcf4a6e70,
456 0xd595d80f, 0x26b2716e, 0xadc666b0, 0x1d153624, 0x3c42d35a, 0x63ff540e, 0xcc5573c0,
457 0x65f9ef17, 0x55bc28f2, 0x80dcc7f7, 0xf46eeddc, 0x5fdcefce, 0x000553f7,
458 }
459 }
460};
461
462//******************************************************************************
463// result = 10^exponent
464//******************************************************************************
465static void BigInt_Pow10(tBigInt * pResult, tU32 exponent)
466{
467 // make sure the exponent is within the bounds of the lookup table data
468 RJ_ASSERT(exponent < 512);
469
470 // create two temporary values to reduce large integer copy operations
471 tBigInt temp1;
472 tBigInt temp2;
473 tBigInt *pCurTemp = &temp1;
474 tBigInt *pNextTemp = &temp2;
475
476 // initialize the result by looking up a 32-bit power of 10 corresponding to the first 3 bits
477 tU32 smallExponent = exponent & 0x7;
478 pCurTemp->SetU32(g_PowerOf10_U32[smallExponent]);
479
480 // remove the low bits that we used for the 32-bit lookup table
481 exponent >>= 3;
482 tU32 tableIdx = 0;
483
484 // while there are remaining bits in the exponent to be processed
485 while (exponent != 0)
486 {
487 // if the current bit is set, multiply it with the corresponding power of 10
488 if(exponent & 1)
489 {
490 // multiply into the next temporary
491 BigInt_Multiply( pNextTemp, *pCurTemp, g_PowerOf10_Big[tableIdx] );
492
493 // swap to the next temporary
494 tBigInt * pSwap = pCurTemp;
495 pCurTemp = pNextTemp;
496 pNextTemp = pSwap;
497 }
498
499 // advance to the next bit
500 ++tableIdx;
501 exponent >>= 1;
502 }
503
504 // output the result
505 *pResult = *pCurTemp;
506}
507
508//******************************************************************************
509// result = in * 10^exponent
510//******************************************************************************
511static void BigInt_MultiplyPow10(tBigInt * pResult, const tBigInt & in, tU32 exponent)
512{
513 // make sure the exponent is within the bounds of the lookup table data
514 RJ_ASSERT(exponent < 512);
515
516 // create two temporary values to reduce large integer copy operations
517 tBigInt temp1;
518 tBigInt temp2;
519 tBigInt *pCurTemp = &temp1;
520 tBigInt *pNextTemp = &temp2;
521
522 // initialize the result by looking up a 32-bit power of 10 corresponding to the first 3 bits
523 tU32 smallExponent = exponent & 0x7;
524 if (smallExponent != 0)
525 {
526 BigInt_Multiply( pCurTemp, in, g_PowerOf10_U32[smallExponent] );
527 }
528 else
529 {
530 *pCurTemp = in;
531 }
532
533 // remove the low bits that we used for the 32-bit lookup table
534 exponent >>= 3;
535 tU32 tableIdx = 0;
536
537 // while there are remaining bits in the exponent to be processed
538 while (exponent != 0)
539 {
540 // if the current bit is set, multiply it with the corresponding power of 10
541 if(exponent & 1)
542 {
543 // multiply into the next temporary
544 BigInt_Multiply( pNextTemp, *pCurTemp, g_PowerOf10_Big[tableIdx] );
545
546 // swap to the next temporary
547 tBigInt * pSwap = pCurTemp;
548 pCurTemp = pNextTemp;
549 pNextTemp = pSwap;
550 }
551
552 // advance to the next bit
553 ++tableIdx;
554 exponent >>= 1;
555 }
556
557 // output the result
558 *pResult = *pCurTemp;
559}
560
561//******************************************************************************
562// result = 2^exponent
563//******************************************************************************
564static inline void BigInt_Pow2(tBigInt * pResult, tU32 exponent)
565{
566 tU32 blockIdx = exponent / 32;
567 RJ_ASSERT( blockIdx < c_BigInt_MaxBlocks );
568
569 for ( tU32 i = 0; i <= blockIdx; ++i)
570 pResult->m_blocks[i] = 0;
571
572 pResult->m_length = blockIdx + 1;
573
574 tU32 bitIdx = (exponent % 32);
575 pResult->m_blocks[blockIdx] |= (1 << bitIdx);
576}
577
578//******************************************************************************
579// This function will divide two large numbers under the assumption that the
580// result is within the range [0,10) and the input numbers have been shifted
581// to satisfy:
582// - The highest block of the divisor is greater than or equal to 8 such that
583// there is enough precision to make an accurate first guess at the quotient.
584// - The highest block of the divisor is less than the maximum value on an
585// unsigned 32-bit integer such that we can safely increment without overflow.
586// - The dividend does not contain more blocks than the divisor such that we
587// can estimate the quotient by dividing the equivalently placed high blocks.
588//
589// quotient = floor(dividend / divisor)
590// remainder = dividend - quotient*divisor
591//
592// pDividend is updated to be the remainder and the quotient is returned.
593//******************************************************************************
594static tU32 BigInt_DivideWithRemainder_MaxQuotient9(tBigInt * pDividend, const tBigInt & divisor)
595{
596 // Check that the divisor has been correctly shifted into range and that it is not
597 // smaller than the dividend in length.
598 RJ_ASSERT( !divisor.IsZero() &&
599 divisor.m_blocks[divisor.m_length-1] >= 8 &&
600 divisor.m_blocks[divisor.m_length-1] < 0xFFFFFFFF &&
601 pDividend->m_length <= divisor.m_length );
602
603 // If the dividend is smaller than the divisor, the quotient is zero and the divisor is already
604 // the remainder.
605 tU32 length = divisor.m_length;
606 if (pDividend->m_length < divisor.m_length)
607 return 0;
608
609 const tU32 * pFinalDivisorBlock = divisor.m_blocks + length - 1;
610 tU32 * pFinalDividendBlock = pDividend->m_blocks + length - 1;
611
612 // Compute an estimated quotient based on the high block value. This will either match the actual quotient or
613 // undershoot by one.
614 tU32 quotient = *pFinalDividendBlock / (*pFinalDivisorBlock + 1);
615 RJ_ASSERT(quotient <= 9);
616
617 // Divide out the estimated quotient
618 if (quotient != 0)
619 {
620 // dividend = dividend - divisor*quotient
621 const tU32 *pDivisorCur = divisor.m_blocks;
622 tU32 *pDividendCur = pDividend->m_blocks;
623
624 tU64 borrow = 0;
625 tU64 carry = 0;
626 do
627 {
628 tU64 product = (tU64)*pDivisorCur * (tU64)quotient + carry;
629 carry = product >> 32;
630
631 tU64 difference = (tU64)*pDividendCur - (product & 0xFFFFFFFF) - borrow;
632 borrow = (difference >> 32) & 1;
633
634 *pDividendCur = difference & 0xFFFFFFFF;
635
636 ++pDivisorCur;
637 ++pDividendCur;
638 } while(pDivisorCur <= pFinalDivisorBlock);
639
640 // remove all leading zero blocks from dividend
641 while (length > 0 && pDividend->m_blocks[length - 1] == 0)
642 --length;
643
644 pDividend->m_length = length;
645 }
646
647 // If the dividend is still larger than the divisor, we overshot our estimate quotient. To correct,
648 // we increment the quotient and subtract one more divisor from the dividend.
649 if ( BigInt_Compare(*pDividend, divisor) >= 0 )
650 {
651 ++quotient;
652
653 // dividend = dividend - divisor
654 const tU32 *pDivisorCur = divisor.m_blocks;
655 tU32 *pDividendCur = pDividend->m_blocks;
656
657 tU64 borrow = 0;
658 do
659 {
660 tU64 difference = (tU64)*pDividendCur - (tU64)*pDivisorCur - borrow;
661 borrow = (difference >> 32) & 1;
662
663 *pDividendCur = difference & 0xFFFFFFFF;
664
665 ++pDivisorCur;
666 ++pDividendCur;
667 } while(pDivisorCur <= pFinalDivisorBlock);
668
669 // remove all leading zero blocks from dividend
670 while (length > 0 && pDividend->m_blocks[length - 1] == 0)
671 --length;
672
673 pDividend->m_length = length;
674 }
675
676 return quotient;
677}
678
679//******************************************************************************
680// result = result << shift
681//******************************************************************************
682static void BigInt_ShiftLeft(tBigInt * pResult, tU32 shift)
683{
684 RJ_ASSERT( shift != 0 );
685
686 tU32 shiftBlocks = shift / 32;
687 tU32 shiftBits = shift % 32;
688
689 // process blocks high to low so that we can safely process in place
690 const tU32 * pInBlocks = pResult->m_blocks;
691 tS32 inLength = pResult->m_length;
692 RJ_ASSERT( inLength + shiftBlocks < c_BigInt_MaxBlocks );
693
694 // check if the shift is block aligned
695 if (shiftBits == 0)
21
←
Assuming 'shiftBits' is equal to 0
→
22
←
Taking true branch
→
696 {
697 // copy blcoks from high to low
698 for (tU32 * pInCur = pResult->m_blocks + inLength, *pOutCur = pInCur + shiftBlocks;
23
←
Loop condition is true. Entering loop body
→
699 pInCur >= pInBlocks;
700 --pInCur, --pOutCur)
701 {
702 *pOutCur = *pInCur;
24
←
Assigned value is garbage or undefined
703 }
704
705 // zero the remaining low blocks
706 for ( tU32 i = 0; i < shiftBlocks; ++i)
707 pResult->m_blocks[i] = 0;
708
709 pResult->m_length += shiftBlocks;
710 }
711 // else we need to shift partial blocks
712 else
713 {
714 tS32 inBlockIdx = inLength - 1;
715 tU32 outBlockIdx = inLength + shiftBlocks;
716
717 // set the length to hold the shifted blocks
718 RJ_ASSERT( outBlockIdx < c_BigInt_MaxBlocks );
719 pResult->m_length = outBlockIdx + 1;
720
721 // output the initial blocks
722 const tU32 lowBitsShift = (32 - shiftBits);
723 tU32 highBits = 0;
724 tU32 block = pResult->m_blocks[inBlockIdx];
725 tU32 lowBits = block >> lowBitsShift;
726 while ( inBlockIdx > 0 )
727 {
728 pResult->m_blocks[outBlockIdx] = highBits | lowBits;
729 highBits = block << shiftBits;
730
731 --inBlockIdx;
732 --outBlockIdx;
733
734 block = pResult->m_blocks[inBlockIdx];
735 lowBits = block >> lowBitsShift;
736 }
737
738 // output the final blocks
739 RJ_ASSERT( outBlockIdx == shiftBlocks + 1 );
740 pResult->m_blocks[outBlockIdx] = highBits | lowBits;
741 pResult->m_blocks[outBlockIdx-1] = block << shiftBits;
742
743 // zero the remaining low blocks
744 for ( tU32 i = 0; i < shiftBlocks; ++i)
745 pResult->m_blocks[i] = 0;
746
747 // check if the terminating block has no set bits
748 if (pResult->m_blocks[pResult->m_length - 1] == 0)
749 --pResult->m_length;
750 }
751}
752
753//******************************************************************************
754// This is an implementation the Dragon4 algorithm to convert a binary number
755// in floating point format to a decimal number in string format. The function
756// returns the number of digits written to the output buffer and the output is
757// not NUL terminated.
758//
759// The floating point input value is (mantissa * 2^exponent).
760//
761// See the following papers for more information on the algorithm:
762// "How to Print Floating-Point Numbers Accurately"
763// Steele and White
764// http://kurtstephens.com/files/p372-steele.pdf
765// "Printing Floating-Point Numbers Quickly and Accurately"
766// Burger and Dybvig
767// http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.72.4656&rep=rep1&type=pdf
768//******************************************************************************
769tU32 Dragon4
770(
771 const tU64 mantissa, // value significand
772 const tS32 exponent, // value exponent in base 2
773 const tU32 mantissaHighBitIdx, // index of the highest set mantissa bit
774 const tB hasUnequalMargins, // is the high margin twice as large as the low margin
775 const tCutoffMode cutoffMode, // how to determine output length
776 tU32 cutoffNumber, // parameter to the selected cutoffMode
777 tC8 * pOutBuffer, // buffer to output into
778 tU32 bufferSize, // maximum characters that can be printed to pOutBuffer
779 tS32 * pOutExponent // the base 10 exponent of the first digit
780)
781{
782 tC8 * pCurDigit = pOutBuffer;
783
784 RJ_ASSERT( bufferSize > 0 );
785
786 // if the mantissa is zero, the value is zero regardless of the exponent
787 if (mantissa == 0)
1
Assuming 'mantissa' is not equal to 0
→
2
←
Taking false branch
→
788 {
789 *pCurDigit = '0';
790 *pOutExponent = 0;
791 return 1;
792 }
793
794 // compute the initial state in integral form such that
795 // value = scaledValue / scale
796 // marginLow = scaledMarginLow / scale
797 tBigInt scale; // positive scale applied to value and margin such that they can be
798 // represented as whole numbers
799 tBigInt scaledValue; // scale * mantissa
800 tBigInt scaledMarginLow; // scale * 0.5 * (distance between this floating-point number and its
801 // immediate lower value)
802
803 // For normalized IEEE floating point values, each time the exponent is incremented the margin also
804 // doubles. That creates a subset of transition numbers where the high margin is twice the size of
805 // the low margin.
806 tBigInt * pScaledMarginHigh;
807 tBigInt optionalMarginHigh;
808
809 if ( hasUnequalMargins )
3
←
Assuming 'hasUnequalMargins' is 0
→
4
←
Taking false branch
→
810 {
811 // if we have no fractional component
812 if (exponent > 0)
813 {
814 // 1) Expand the input value by multiplying out the mantissa and exponent. This represents
815 // the input value in its whole number representation.
816 // 2) Apply an additional scale of 2 such that later comparisons against the margin values
817 // are simplified.
818 // 3) Set the margin value to the lowest mantissa bit's scale.
819
820 // scaledValue = 2 * 2 * mantissa*2^exponent
821 scaledValue.SetU64( 4 * mantissa );
822 BigInt_ShiftLeft( &scaledValue, exponent );
823
824 // scale = 2 * 2 * 1
825 scale.SetU32( 4 );
826
827 // scaledMarginLow = 2 * 2^(exponent-1)
828 BigInt_Pow2( &scaledMarginLow, exponent );
829
830 // scaledMarginHigh = 2 * 2 * 2^(exponent-1)
831 BigInt_Pow2( &optionalMarginHigh, exponent + 1 );
832 }
833 // else we have a fractional exponent
834 else
835 {
836 // In order to track the mantissa data as an integer, we store it as is with a large scale
837
838 // scaledValue = 2 * 2 * mantissa
839 scaledValue.SetU64( 4 * mantissa );
840
841 // scale = 2 * 2 * 2^(-exponent)
842 BigInt_Pow2(&scale, -exponent + 2 );
843
844 // scaledMarginLow = 2 * 2^(-1)
845 scaledMarginLow.SetU32( 1 );
846
847 // scaledMarginHigh = 2 * 2 * 2^(-1)
848 optionalMarginHigh.SetU32( 2 );
849 }
850
851 // the high and low margins are different
852 pScaledMarginHigh = &optionalMarginHigh;
853 }
854 else
855 {
856 // if we have no fractional component
857 if (exponent > 0)
5
←
Assuming 'exponent' is <= 0
→
6
←
Taking false branch
→
858 {
859 // 1) Expand the input value by multiplying out the mantissa and exponent. This represents
860 // the input value in its whole number representation.
861 // 2) Apply an additional scale of 2 such that later comparisons against the margin values
862 // are simplified.
863 // 3) Set the margin value to the lowest mantissa bit's scale.
864
865 // scaledValue = 2 * mantissa*2^exponent
866 scaledValue.SetU64( 2 * mantissa );
867 BigInt_ShiftLeft( &scaledValue, exponent );
868
869 // scale = 2 * 1
870 scale.SetU32( 2 );
871
872 // scaledMarginLow = 2 * 2^(exponent-1)
873 BigInt_Pow2( &scaledMarginLow, exponent );
874 }
875 // else we have a fractional exponent
876 else
877 {
878 // In order to track the mantissa data as an integer, we store it as is with a large scale
879
880 // scaledValue = 2 * mantissa
881 scaledValue.SetU64( 2 * mantissa );
882
883 // scale = 2 * 2^(-exponent)
884 BigInt_Pow2(&scale, -exponent + 1 );
885
886 // scaledMarginLow = 2 * 2^(-1)
887 scaledMarginLow.SetU32( 1 );
888 }
889
890 // the high and low margins are equal
891 pScaledMarginHigh = &scaledMarginLow;
892 }
893
894 // Compute an estimate for digitExponent that will be correct or undershoot by one.
895 // This optimization is based on the paper "Printing Floating-Point Numbers Quickly and Accurately"
896 // by Burger and Dybvig http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.72.4656&rep=rep1&type=pdf
897 // We perform an additional subtraction of 0.69 to increase the frequency of a failed estimate
898 // because that lets us take a faster branch in the code. 0.69 is chosen because 0.69 + log10(2) is
899 // less than one by a reasonable epsilon that will account for any floating point error.
900 //
901 // We want to set digitExponent to floor(log10(v)) + 1
902 // v = mantissa*2^exponent
903 // log2(v) = log2(mantissa) + exponent;
904 // log10(v) = log2(v) * log10(2)
905 // floor(log2(v)) = mantissaHighBitIdx + exponent;
906 // log10(v) - log10(2) < (mantissaHighBitIdx + exponent) * log10(2) <= log10(v)
907 // log10(v) < (mantissaHighBitIdx + exponent) * log10(2) + log10(2) <= log10(v) + log10(2)
908 // floor( log10(v) ) < ceil( (mantissaHighBitIdx + exponent) * log10(2) ) <= floor( log10(v) ) + 1
909 const tF64 log10_2 = 0.30102999566398119521373889472449;
910 tS32 digitExponent = (tS32)(ceil(tF64((tS32)mantissaHighBitIdx + exponent) * log10_2 - 0.69));
911
912 // if the digit exponent is smaller than the smallest desired digit for fractional cutoff,
913 // pull the digit back into legal range at which point we will round to the appropriate value.
914 // Note that while our value for digitExponent is still an estimate, this is safe because it
915 // only increases the number. This will either correct digitExponent to an accurate value or it
916 // will clamp it above the accurate value.
917 if (cutoffMode == CutoffMode_FractionLength && digitExponent <= -(tS32)cutoffNumber)
7
←
Assuming 'cutoffMode' is not equal to CutoffMode_FractionLength
→
918 {
919 digitExponent = -(tS32)cutoffNumber + 1;
920 }
921
922 // Divide value by 10^digitExponent.
923 if (digitExponent > 0)
8
←
Assuming 'digitExponent' is <= 0
→
9
←
Taking false branch
→
924 {
925 // The exponent is positive creating a division so we multiply up the scale.
926 tBigInt temp;
927 BigInt_MultiplyPow10( &temp, scale, digitExponent );
928 scale = temp;
929 }
930 else if (digitExponent < 0)
10
←
Assuming 'digitExponent' is >= 0
→
11
←
Taking false branch
→
931 {
932 // The exponent is negative creating a multiplication so we multiply up the scaledValue,
933 // scaledMarginLow and scaledMarginHigh.
934 tBigInt pow10;
935 BigInt_Pow10( &pow10, -digitExponent);
936
937 tBigInt temp;
938 BigInt_Multiply( &temp, scaledValue, pow10);
939 scaledValue = temp;
940
941 BigInt_Multiply( &temp, scaledMarginLow, pow10);
942 scaledMarginLow = temp;
943
944 if (pScaledMarginHigh != &scaledMarginLow)
945 BigInt_Multiply2( pScaledMarginHigh, scaledMarginLow );
946 }
947
948 // If (value + marginHigh) >= 1, our estimate for digitExponent was too low
949 tBigInt scaledValueHigh;
950 BigInt_Add( &scaledValueHigh, scaledValue, *pScaledMarginHigh );
951 if( BigInt_Compare(scaledValueHigh,scale) >= 0 )
12
←
Taking false branch
→
952 {
953 // The exponent estimate was incorrect.
954 // Increment the exponent and don't perform the premultiply needed
955 // for the first loop iteration.
956 digitExponent = digitExponent + 1;
957 }
958 else
959 {
960 // The exponent estimate was correct.
961 // Multiply larger by the output base to prepare for the first loop iteration.
962 BigInt_Multiply10( &scaledValue );
963 BigInt_Multiply10( &scaledMarginLow );
964 if (pScaledMarginHigh != &scaledMarginLow)
13
←
Taking false branch
→
965 BigInt_Multiply2( pScaledMarginHigh, scaledMarginLow );
966 }
967
968 // Compute the cutoff exponent (the exponent of the final digit to print).
969 // Default to the maximum size of the output buffer.
970 tS32 cutoffExponent = digitExponent - bufferSize;
971 switch(cutoffMode)
14
←
Control jumps to 'case CutoffMode_TotalLength:' at line 978
→
972 {
973 // print digits until we pass the accuracy margin limits or buffer size
974 case CutoffMode_Unique:
975 break;
976
977 // print cutoffNumber of digits or until we reach the buffer size
978 case CutoffMode_TotalLength:
979 {
980 tS32 desiredCutoffExponent = digitExponent - (tS32)cutoffNumber;
981 if (desiredCutoffExponent > cutoffExponent)
15
←
Taking false branch
→
982 cutoffExponent = desiredCutoffExponent;
983 }
984 break;
16
←
 Execution continues on line 997
→
985
986 // print cutoffNumber digits past the decimal point or until we reach the buffer size
987 case CutoffMode_FractionLength:
988 {
989 tS32 desiredCutoffExponent = -(tS32)cutoffNumber;
990 if (desiredCutoffExponent > cutoffExponent)
991 cutoffExponent = desiredCutoffExponent;
992 }
993 break;
994 }
995
996 // Output the exponent of the first digit we will print
997 *pOutExponent = digitExponent-1;
998
999 // In preparation for calling BigInt_DivideWithRemainder_MaxQuotient9(),
1000 // we need to scale up our values such that the highest block of the denominator
1001 // is greater than or equal to 8. We also need to guarantee that the numerator
1002 // can never have a length greater than the denominator after each loop iteration.
1003 // This requires the highest block of the denominator to be less than or equal to
1004 // 429496729 which is the highest number that can be multiplied by 10 without
1005 // overflowing to a new block.
1006 RJ_ASSERT( scale.GetLength() > 0 );
1007 tU32 hiBlock = scale.GetBlock( scale.GetLength() - 1 );
1008 if (hiBlock < 8 || hiBlock > 429496729)
17
←
Assuming 'hiBlock' is >= 8
→
18
←
Assuming 'hiBlock' is > 429496729
→
19
←
Taking true branch
→
1009 {
1010 // Perform a bit shift on all values to get the highest block of the denominator into
1011 // the range [8,429496729]. We are more likely to make accurate quotient estimations
1012 // in BigInt_DivideWithRemainder_MaxQuotient9() with higher denominator values so
1013 // we shift the denominator to place the highest bit at index 27 of the highest block.
1014 // This is safe because (2^28 - 1) = 268435455 which is less than 429496729. This means
1015 // that all values with a highest bit at index 27 are within range.
1016 tU32 hiBlockLog2 = LogBase2(hiBlock);
1017 RJ_ASSERT(hiBlockLog2 < 3 || hiBlockLog2 > 27);
1018 tU32 shift = (32 + 27 - hiBlockLog2) % 32;
1019
1020 BigInt_ShiftLeft( &scale, shift );
20
←
Calling 'BigInt_ShiftLeft'
→
1021 BigInt_ShiftLeft( &scaledValue, shift);
1022 BigInt_ShiftLeft( &scaledMarginLow, shift);
1023 if (pScaledMarginHigh != &scaledMarginLow)
1024 BigInt_Multiply2( pScaledMarginHigh, scaledMarginLow );
1025 }
1026
1027 // These values are used to inspect why the print loop terminated so we can properly
1028 // round the final digit.
1029 tB low; // did the value get within marginLow distance from zero
1030 tB high; // did the value get within marginHigh distance from one
1031 tU32 outputDigit; // current digit being output
1032
1033 if (cutoffMode == CutoffMode_Unique)
1034 {
1035 // For the unique cutoff mode, we will try to print until we have reached a level of
1036 // precision that uniquely distinguishes this value from its neighbors. If we run
1037 // out of space in the output buffer, we terminate early.
1038 for (;;)
1039 {
1040 digitExponent = digitExponent-1;
1041
1042 // divide out the scale to extract the digit
1043 outputDigit = BigInt_DivideWithRemainder_MaxQuotient9(&scaledValue, scale);
1044 RJ_ASSERT( outputDigit < 10 );
1045
1046 // update the high end of the value
1047 BigInt_Add( &scaledValueHigh, scaledValue, *pScaledMarginHigh );
1048
1049 // stop looping if we are far enough away from our neighboring values
1050 // or if we have reached the cutoff digit
1051 low = BigInt_Compare(scaledValue, scaledMarginLow) < 0;
1052 high = BigInt_Compare(scaledValueHigh, scale) > 0;
1053 if (low | high | (digitExponent == cutoffExponent))
1054 break;
1055
1056 // store the output digit
1057 *pCurDigit = (tC8)('0' + outputDigit);
1058 ++pCurDigit;
1059
1060 // multiply larger by the output base
1061 BigInt_Multiply10( &scaledValue );
1062 BigInt_Multiply10( &scaledMarginLow );
1063 if (pScaledMarginHigh != &scaledMarginLow)
1064 BigInt_Multiply2( pScaledMarginHigh, scaledMarginLow );
1065 }
1066 }
1067 else
1068 {
1069 // For length based cutoff modes, we will try to print until we
1070 // have exhausted all precision (i.e. all remaining digits are zeros) or
1071 // until we reach the desired cutoff digit.
1072 low = false;
1073 high = false;
1074
1075 for (;;)
1076 {
1077 digitExponent = digitExponent-1;
1078
1079 // divide out the scale to extract the digit
1080 outputDigit = BigInt_DivideWithRemainder_MaxQuotient9(&scaledValue, scale);
1081 RJ_ASSERT( outputDigit < 10 );
1082
1083 if ( scaledValue.IsZero() | (digitExponent == cutoffExponent) )
1084 break;
1085
1086 // store the output digit
1087 *pCurDigit = (tC8)('0' + outputDigit);
1088 ++pCurDigit;
1089
1090 // multiply larger by the output base
1091 BigInt_Multiply10(&scaledValue);
1092 }
1093 }
1094
1095 // round off the final digit
1096 // default to rounding down if value got too close to 0
1097 tB roundDown = low;
1098
1099 // if it is legal to round up and down
1100 if (low == high)
1101 {
1102 // round to the closest digit by comparing value with 0.5. To do this we need to convert
1103 // the inequality to large integer values.
1104 // compare( value, 0.5 )
1105 // compare( scale * value, scale * 0.5 )
1106 // compare( 2 * scale * value, scale )
1107 BigInt_Multiply2(&scaledValue);
1108 tS32 compare = BigInt_Compare(scaledValue, scale);
1109 roundDown = compare < 0;
1110
1111 // if we are directly in the middle, round towards the even digit (i.e. IEEE rouding rules)
1112 if (compare == 0)
1113 roundDown = (outputDigit & 1) == 0;
1114 }
1115
1116 // print the rounded digit
1117 if (roundDown)
1118 {
1119 *pCurDigit = (tC8)('0' + outputDigit);
1120 ++pCurDigit;
1121 }
1122 else
1123 {
1124 // handle rounding up
1125 if (outputDigit == 9)
1126 {
1127 // find the first non-nine prior digit
1128 for (;;)
1129 {
1130 // if we are at the first digit
1131 if (pCurDigit == pOutBuffer)
1132 {
1133 // output 1 at the next highest exponent
1134 *pCurDigit = '1';
1135 ++pCurDigit;
1136 *pOutExponent += 1;
1137 break;
1138 }
1139
1140 --pCurDigit;
1141 if (*pCurDigit != '9')
1142 {
1143 // increment the digit
1144 *pCurDigit += 1;
1145 ++pCurDigit;
1146 break;
1147 }
1148 }
1149 }
1150 else
1151 {
1152 // values in the range [0,8] can perform a simple round up
1153 *pCurDigit = (tC8)('0' + outputDigit + 1);
1154 ++pCurDigit;
1155 }
1156 }
1157
1158 // return the number of digits output
1159 RJ_ASSERT(pCurDigit - pOutBuffer <= (tPtrDiff)bufferSize);
1160 return pCurDigit - pOutBuffer;
1161}