cmathmultiplicationmodulo

Modulo multiplication (in C)


Specifically: I have two unsigned integers (a,b) and I want to calculate (a*b)%UINT_MAX (UINT_MAX is defined as the maximal unsigned int). What is the best way to do so?

Background: I need to write a module for linux that will emulate a geometric sequence, reading from it will give me the next element (modulo UINT_MAX), the only solution I found is to add the current element to itself times, while adding is done using the following logic:(that I use for the arithmetic sequence)

for(int i=0; i<b; ++i){
  if(UINT_MAX - current_value > difference) {
    current_value += difference;
  } else {
    current_value = difference - (UINT_MAX - current_value);
  }

when current_value = a in the first iteration (and is updated in every iteration, and difference = a (always). Obviously this is not a intelligent solution. How would an intelligent person achieve this?

Thanks!


Solution

  • As has been mentioned, if you have a type of twice the width available, just use that, here

    (unsigned int)(((unsigned long long)a * b) % UINT_MAX)
    

    if int is 32 bits and long long 64 (or more). If you have no larger type, you can split the factors at half the bit-width, multiply and reduce the parts, finally assemble it. Illustrated for 32-bit unsigned here:

    a_low = a & 0xFFFF;  // low 16 bits of a
    a_high = a >> 16;    // high 16 bits of a, shifted in low half
    b_low = b & 0xFFFF;
    b_high = b >> 16;
    /*
     * Now a = (a_high * 65536 + a_low), b = (b_high * 65536 + b_low)
     * Thus a*b = (a_high * b_high) * 65536 * 65536
     *          + (a_high * b_low + a_low * b_high) * 65536
     *          + a_low * b_low
     *
     * All products a_i * b_j are at most (65536 - 1) * (65536 - 1) = UINT_MAX - 2 * 65536 + 2
     * The high product reduces to
     * (a_high * b_high) * (UINT_MAX + 1) = (a_high * b_high)
     * The middle products are a bit trickier, but splitting again solves:
     * m1 = a_high * b_low;
     * m1_low = m1 & 0xFFFF;
     * m1_high = m1 >> 16;
     * Then m1 * 65536 = m1_high * (UINT_MAX + 1) + m1_low * 65536 = m1_high + m1_low * 65536
     * Similar for a_low * b_high
     * Finally, add the parts and take care of overflow
     */
    m1 = a_high * b_low;
    m2 = a_low * b_high;
    m1_low = m1 & 0xFFFF;
    m1_high = m1 >> 16;
    m2_low = m2 & 0xFFFF;
    m2_high = m2 >> 16;
    result = a_high * b_high;
    temp = result + ((m1_low << 16) | m1_high);
    if (temp < result)    // overflow
    {
        result = temp+1;
    }
    else
    {
        result = temp;
    }
    if (result == UINT_MAX)
    {
        result = 0;
    }
    // I'm too lazy to type out the rest, you get the gist, I suppose.
    

    Of course, if what you need is actually reduction modulo UINT_MAX + 1, as @Toad assumes,then that's just what multiplication of unsigned int does.