Thread: Algorithm Inquiry

  1. #1
    Registered User
    Join Date
    Jul 2010
    Location
    Oklahoma
    Posts
    107

    Question 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");
       b = 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;
                }
            }
            
            bigintadd(r, tmp, swp);
            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
    Last edited by new_ink2001; 12-26-2011 at 01:00 AM. Reason: Placement of attachment link
    Kept the text books....
    Went interdisciplinary after college....
    Still looking for a real job since 2005....

    During the interim, I may be reached at ELance, vWorker, FreeLancer, oDesk and WyzAnt.

  2. #2
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    39,660
    > 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.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  3. #3
    Registered User
    Join Date
    Jul 2010
    Location
    Oklahoma
    Posts
    107
    Thank you. I hadn't caught dc yet. I was still using bc.
    Kept the text books....
    Went interdisciplinary after college....
    Still looking for a real job since 2005....

    During the interim, I may be reached at ELance, vWorker, FreeLancer, oDesk and WyzAnt.

  4. #4
    Registered User
    Join Date
    Jul 2010
    Location
    Oklahoma
    Posts
    107
    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...
    Kept the text books....
    Went interdisciplinary after college....
    Still looking for a real job since 2005....

    During the interim, I may be reached at ELance, vWorker, FreeLancer, oDesk and WyzAnt.

  5. #5
    Registered User
    Join Date
    Jul 2010
    Location
    Oklahoma
    Posts
    107
    I guess I missed something there about line 120...when tmp is still a zero. Thank you for the tip Salem.
    Kept the text books....
    Went interdisciplinary after college....
    Still looking for a real job since 2005....

    During the interim, I may be reached at ELance, vWorker, FreeLancer, oDesk and WyzAnt.

  6. #6
    Registered User
    Join Date
    Jul 2010
    Location
    Oklahoma
    Posts
    107
    It looks like line 94 is missing a modulus by ten too.
    Kept the text books....
    Went interdisciplinary after college....
    Still looking for a real job since 2005....

    During the interim, I may be reached at ELance, vWorker, FreeLancer, oDesk and WyzAnt.

  7. #7
    Registered User
    Join Date
    Jul 2010
    Location
    Oklahoma
    Posts
    107
    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.

    It should read something like:

    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);
            bigintadd(r, tmp, swp);
            bigintcpy(swp, r);
    
            memset(tmp->m, ZERO, tmp->u*ARBITRARYMATH_ALLOCATION_UNIT*sizeof(bigdigit));
            tmp->mlen = 1;
            tmp->s = zero;
        }
    Thank you for the attention.
    Kept the text books....
    Went interdisciplinary after college....
    Still looking for a real job since 2005....

    During the interim, I may be reached at ELance, vWorker, FreeLancer, oDesk and WyzAnt.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. GNU strip inquiry
    By Chris87 in forum Tech Board
    Replies: 2
    Last Post: 07-05-2009, 08:04 PM
  2. Bitflags inquiry
    By Chris87 in forum C Programming
    Replies: 4
    Last Post: 01-03-2009, 06:22 PM
  3. Loops: The 'any' inquiry
    By chubbs1900 in forum C++ Programming
    Replies: 6
    Last Post: 12-10-2007, 10:35 AM
  4. for loop inquiry
    By cunnus88 in forum C Programming
    Replies: 36
    Last Post: 10-20-2005, 11:03 PM
  5. Prolog inquiry
    By magnum38 in forum C Programming
    Replies: 2
    Last Post: 05-25-2002, 10:05 AM