1. ## Algorithm Inquiry

Hello again everyone,

I have been promising myself that I would do some work on a library for large numbers since college. I read about them in a textbook by Binstock and Rex, but I didn't like the static nature of their implementation and I had a few other issues so I wrote my own.

Here is an excerpt from the main driver, test is currently set to 9999:

Code:
```   a = bigintread("9");

for(i = 0; i < test; i++, bigintcpy(c, b))
{
bigintmult(a, b, c);
}

printf( "------------------------\n" );
printf( "for( i = 0; i < test; i++ ) c = c*a\n");
bigintprint(stdout, a, 0);
printf( "\n" );
bigintprint(stdout, b, 0);
printf( "\n" );
bigintprint(stdout, c, 0);
printf( "\n" );```
The program kept up with the system's graphical calculator (xcalc) until about 9^35 when the calculator produced an error message. The loop above produced 9^10000, it is in the file that I have attached. It has about the right number of digits 9494, so I'm rather confident with it.

Does anyone have that figure handy to verify it?

The program was running in the debugger for about seven hours before producing that figure, one factor of nine at a time lol. I think maybe it would be prudent to build a database for previously performed calculations. Then use a separate thread to find them if any while waiting for a calculation thread to finish, and using the first that came back successful of course.

The code for multiplication is here:

Code:
```/*** bigintmult.c **********************************************
*
*  Multiply two big integers a and b.  Place the result into c.
*
*/

#include "bigint.h"

bigint *bigintmult( const bigint *a, const bigint *b, bigint *c )
{
bigint *r;          // pointer to the current register (product)
bigint *swp;        //
bigint *tmp;

short i;
short j;
short k;
short h;

int carry;

// verify/arrange for the register
if( NULL == c )
{
r = bigintcrt();
}
else
{
r = bigintzero(c);
}

// test that the operands are sane
if(not_a_number == a->s || not_a_number == b->s)
{
// NaN * b = NaN ^ a * NaN = NaN ^ NaN * NaN = NaN
r->s = not_a_number;

return r;
}

if(infinite_positive == a->s || infinite_negative == a->s)
{
// infinite * b = infinite
r->s = a->s;

return r;
}

if(infinite_positive == b->s || infinite_negative == b->s)
{
// a * infinite = infinite
r->s = b->s;

return r;
}

if(zero == a->s || zero == b->s)
{
return r;
}

// a != 0 ^ b != 0, a,b non-ideal integers

// clear and resize the register in case of overflow
tmp = bigintcrt();
bigintrealloc(tmp, 2*MAX(a->u, b->u) + 1);

swp = bigintcrt();

// perform multiplication
for(i = 0; i < b->mlen; ++i)
{
if(ZERO == b->m[i])
{
continue;
}

for(j = 0; j < a->mlen; j++)
{
tmp->m[i+j] = a->m[j]*b->m[i];

if(zero == tmp->s && ZERO != tmp->m[j+i])
{
tmp->s = positive;
}
}

// -- normalize tmp for addition to the result register --
for(k = 0; k < tmp->u*ARBITRARYMATH_ALLOCATION_UNIT; ++k)
{
if(9 < tmp->m[k])
{
carry = tmp->m[k]/10;
tmp->m[k] = tmp->m[k];

for(h = 1;
0 < carry
&& (k+h) < tmp->u*ARBITRARYMATH_ALLOCATION_UNIT;
++h, carry /= 10) // high_order_truncation will occur here
{
tmp->m[k+h] += (carry);

if(h+k+1 == tmp->u*ARBITRARYMATH_ALLOCATION_UNIT
&& (tmp->u+1)*ARBITRARYMATH_ALLOCATION_UNIT <= SHRT_MAX)
{
tmp->mlen = tmp->u*ARBITRARYMATH_ALLOCATION_UNIT;
bigintrealloc(tmp, tmp->u+1);
}
}
}
}

for(tmp->mlen = tmp->u*ARBITRARYMATH_ALLOCATION_UNIT;
0 < tmp->mlen; --tmp->mlen)
{
if(ZERO != tmp->m[tmp->mlen-1])
{
break;
}
}

bigintcpy(swp, r);

memset(tmp->m, ZERO, tmp->u*ARBITRARYMATH_ALLOCATION_UNIT*sizeof(bigdigit));
tmp->mlen = 1;
tmp->s = zero;
}

// select the correct sign
if( a->s != b->s )
{
r->s = negative;
}
else
{
r->s = positive;
}

// clean up
bigintdes(tmp);
bigintdes(swp);

bigintlength(r);

if(ceil(((float)r->mlen)/ARBITRARYMATH_ALLOCATION_UNIT) < r->u)
{
bigintrealloc(r, ceil(((float)r->mlen)/ARBITRARYMATH_ALLOCATION_UNIT));
}

return r;
}```
Other alternative implementations also came to mind. Like using a queue rather than a table format for the number datastructure, but it made more sense to me that it grow a 'page' at a time rather than double in size as a queue for a stream might according to the suggestion in the algorithm analysis textbook. If a carry had more than an allocation unit of figures, the instrinsic data type would have already overflowed is my justification. What sort of advantages might I have missed from using a queue for representation?

Also, in order to wring out the last bit of speed from this Pentium 4, I think that the point of parallelism is at the intermediate products before summing. That doesn't help with the situation when multiplication occurs one unit at a time. Perhaps partitioning the figure by place value is more appropriate? Tips here would be great too. I had a great parallel algorithm analysis textbook while I was in college, but a friend indeed needed one.

Thank you for your time and effort and it is great to post with you again.

Best Regards,

New Ink -- Henry
Enclosure: cprogramming_attachment.txt 2. > The program kept up with the system's graphical calculator (xcalc) until about 9^35 when the calculator produced an error message
dc(1): arbitrary precision calculator - Linux man page
dc is a reverse-polish desk calculator which supports unlimited precision arithmetic. It also allows you to define and call macros. Normally dc reads from the standard input; if any command arguments are given to it, they are filenames, and dc reads and executes the contents of the files before reading from standard input. All normal output is to standard output; all error output is to standard error.
Since both your program and dc use the console stdio, you might be able to drive both together in some script, and verify the answers automatically. 3. Thank you. I hadn't caught dc yet. I was still using bc. 4. According to dc's output I'm off a little less than a hundred digits. It should be 9543 digits long. Back to the drawing board... 5. I guess I missed something there about line 120...when tmp is still a zero. Thank you for the tip Salem. 6. It looks like line 94 is missing a modulus by ten too. 7. Hello again every one,

I think I was being a little too careful there about line 103. If k = tmp->u*ARBITRARYMATH_ALLOCATION_UNIT - 1

Then h + k + 1 is tmp->u*ARBITRARYMATH_ALLOCATION_UNIT + 1, and the number suffers a high order overflow, but never checks to see if there is more space.

Code:
```    // perform multiplication
for(i = 0; i < small->mlen; ++i)
{
if(ZERO == small->m[i])
{
continue;
}

for(j = 0; j < large->mlen; j++)
{
tmp->m[i+j] = large->m[j]*small->m[i];

if(zero == tmp->s && ZERO != tmp->m[j+i])
{
tmp->s = positive;
}
}

// -- normalize tmp for addition to the result register --
for(k = 0; k < tmp->u*ARBITRARYMATH_ALLOCATION_UNIT; ++k)
{
if(9 < tmp->m[k])
{
carry = tmp->m[k]/10;
tmp->m[k] = tmp->m[k]%10;

for(h = 1; 0 < carry; ++h, carry /= 10)
// high_order_truncation will occur here
{
if(h+k == tmp->u*ARBITRARYMATH_ALLOCATION_UNIT
&& (tmp->u+1)*ARBITRARYMATH_ALLOCATION_UNIT <= SHRT_MAX)
{
tmp->mlen = tmp->u*ARBITRARYMATH_ALLOCATION_UNIT;
bigintrealloc(tmp, tmp->u+1);
}
tmp->m[k+h] += (carry%10);
}
}
}

bigintlen(tmp); Popular pages Recent additions 