### Bug Summary

File: | dconv/dragon4.cpp
Location: | line 631, column 51
Description: | The left operand of '-' is a garbage value

### 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 | //******************************************************************************
61 | const 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 | //******************************************************************************
70 | struct 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 | //******************************************************************************
148 | static 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 | //******************************************************************************
173 | static 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 | //******************************************************************************
239 | static 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 | //******************************************************************************
306 | static 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 | //******************************************************************************
337 | static 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 | //******************************************************************************
368 | static 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 | //******************************************************************************
394 | static 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 )

14 |

Loop condition is false. Execution continues on line 408

402 {
403 tU64 product = (tU64)(*pCur) * 10ull + carry;
404 (*pCur) = (tU32)(product & 0xFFFFFFFF);
405 carry = product >> 32;
406 }
407
408 if (carry != 0)

15 |

Taking false branch
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)

34 |

Taking false branch
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)

35 |

Assuming 'quotient' is not equal to 0

36 |

Taking true branch
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

37 |

Loop condition is true. Execution continues on line 628
627 {
628 tU64 product = (tU64)*pDivisorCur * (tU64)quotient + carry;
629 carry = product >> 32;
630
631 tU64 difference = (tU64)*pDividendCur - (product & 0xFFFFFFFF) - borrow;

38 |

The left operand of '-' is a garbage value
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)

25 |

Taking false branch
696 {
697 // copy blcoks from high to low
698 for (tU32 * pInCur = pResult->m_blocks + inLength, *pOutCur = pInCur + shiftBlocks;
699 pInCur >= pInBlocks;
700 --pInCur, --pOutCur)
701 {
702 *pOutCur = *pInCur;
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 )

26 |

Loop condition is false. Execution continues on line 740
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)

27 |

Loop condition is false. Execution continues on line 748
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)

28 |

Taking false branch
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 );

13 |

Calling 'BigInt_Multiply10'

16 |

Returning from 'BigInt_Multiply10'
963 BigInt_Multiply10( &scaledMarginLow );
964 if (pScaledMarginHigh != &scaledMarginLow)

17 |

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)

18 |

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)

19 |

Taking false branch
982 cutoffExponent = desiredCutoffExponent;
983 }
984 break;

20 |

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)

21 |

Assuming 'hiBlock' is >= 8

22 |

Assuming 'hiBlock' is > 429496729

23 |

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 );
1021 BigInt_ShiftLeft( &scaledValue, shift);

24 |

Calling 'BigInt_ShiftLeft'

29 |

Returning from 'BigInt_ShiftLeft'
1022 BigInt_ShiftLeft( &scaledMarginLow, shift);
1023 if (pScaledMarginHigh != &scaledMarginLow)

30 |

Taking false branch
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)

31 |

Taking false branch
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 (;;)

32 |

Loop condition is true. Entering loop body
1076 {
1077 digitExponent = digitExponent-1;
1078
1079 // divide out the scale to extract the digit
1080 outputDigit = BigInt_DivideWithRemainder_MaxQuotient9(&scaledValue, scale);

33 |

Calling 'BigInt_DivideWithRemainder_MaxQuotient9'
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}