Big Integer Arithmetic

Most of computer applications operates on small integers which can be represented by 32 bit integers. However, there are times we’ll need to deal with extremely big integers which can only be represented by hundreds of bits. One of such application is cryptography.

This post gives two simple implementations of the big integer arithmetic. Both implementations use char array to store the big integers. The first one stores the digits of the integer as char, and the second one stores the digits as single digit integer values.

String Representation
We define a big integer as the structure below,

typedef struct {

   int sign;

   int length;

   char dg[MAXDG];

} BN;

For example, given two input integers “111222333444555” and “12345”, we can store them into two BNs a and b, with a.length = 15, a.dg[0] = ‘1’, a.dg[1]=’1’ …, a.dg[14]=’5’, a.sign = 1; b.length=5, b.dg[0]=’1’, b.dg[1]=’2’…a.dg[4]=’5’, b.sign=1.

Addition: We add the two integers digit by digit from right to left. If the sum of the two digits is bigger than ‘9’. Then a carry of 1 is produced for the next position. And the carry is added to the sum of the digits at the next position. We store the digits of the result from left to right and then reverse it after it’s done.

Subtraction: Subtraction is similar to addition. We minus the digits from the two integers from right to left. If the resulted digit is less than 0, we borrow 1 from the next position.

Multiplication: A naive approach is to use addition to handle multiplication. But if both inputs are big, then the algorithm is very slow.

A better approach is to handle a single digit multiplication using addition, then shift the multiplicand by 1 (which is equivalent to x10) and handle the next digit multiplication.

For example, if the two inputs are 12345 and 57, we treat the multiplication as 12345 x 7 + 123450 x 5, then each of the two multiplications can be handled using additions. In this way, the number of additions needed are reduced to the sum of all digits of the multiplier.

Division: Again the naive approach is to use subtraction. A better alternative is trying to subtract the divisor from part of the dividend from right to left. For example, to compute 234 / 11, we start by computing 23 / 11 using subtraction. The quotient is 2 and the remainder is 1.  We then compute (1*10 + 4) / 11 using subtraction, the quotient will be 1. So the result is 21.

Below is a C implementation of the algorithms described above,

/**

array implementation of bignum

**/

#include <stdio.h>

#include <string.h>

 

#define MAXDG 500

#define PLUS 1

#define MINUS -1

 

typedef struct {

    int sign;

    int length;

    char dg[MAXDG];

} BN;

 

int cmp(BN a, BN b);

void sub(BN a, BN b, BN *c);

void add(BN a, BN b, BN *c);

void reverse(BN *c);

void rm_tzeros(BN *c);

void rm_lzeros(BN *c);

void init(BN *a);

 

int cmp(BN a, BN b) {

    if (a.sign == PLUS && b.sign == MINUS) {

        return 1;

    } else if (a.sign == MINUS && b.sign == PLUS) {

        return -1;

    } else {

        //a and b have same sign

        if (a.length > b.length) {

            return 1*a.sign;

        } else if (a.length < b.length) {

            return -1*a.sign;

        } else {

            //a and b have same sign, same length

            int i;

            for (i = 0; i < a.length; ++i) {

                if (a.dg[i] > b.dg[i]) {

                    return 1*a.sign;

                } else if (a.dg[i] < b.dg[i]) {

                    return -1*a.sign;

                } 

            }

            return 0; //equal

        }

    }

}

 

void reverse(BN *c) {

    char tmp;

    int i;

    for (i = 0; i < c->length/2; ++i) {

        tmp = c->dg[i];

        c->dg[i] = c->dg[c->length-i-1];

        c->dg[c->length-i-1] = tmp;

    }  

}

 

void rm_lzeros(BN *c) {

    int i, j;

    for (i = 0; i < c->length; ++i) {

        if (c->dg[i] != '0') {

            break;

        }

    }

    for (j = 0; i < c->length; ++i, ++j) {

        c->dg[j] = c->dg[i];

    }

    c->length = j == 0?1:j;

    c->dg[c->length] = '';

}

 

void rm_tzeros(BN *c) {

    int i, j = 0;

    for (i = c->length-1; i > 0; --i) {

        if (c->dg[i] == '0') {

            j++;

        } else {

            break;

        }

    }

    c->length -= j;

    c->dg[c->length] = '';

}

 

void sub(BN a, BN b, BN *c) {

   int i, j, k;

   int borrow = 0;

   if (a.sign == MINUS) {

       if (b.sign == MINUS) {

           sub(b, a, c);

           c->sign *= -1;

       } else {

           b.sign = MINUS;

           add(a, b, c);

       }

       return;

   } else {

       if (b.sign == MINUS) {

           b.sign = PLUS;

           add(a, b, c);

           return;

       } 

   }

   //a.sign == PLUS, b.sign == PLUS

   if (cmp(a, b) < 0) {

      sub(b, a, c);

      c->sign = MINUS;

      return; 

   } else if (cmp(a, b) == 0) {

      c->dg[0] = '0';

      c->dg[1] = '';

      c->length = 1;

      c->sign = PLUS;

      return;

   } 

   //a.sign == PLUS, b.sign == PLUS, a > b

   c->sign = PLUS;

   for (i = a.length-1, j = b.length-1, k = 0; i >= 0 && j >= 0; --i, --j, ++k) {

      c->dg[k] = a.dg[i] - borrow - b.dg[j];

      if (c->dg[k] < 0) {

          c->dg[k] += 10;

          borrow = 1;

      } else {

          borrow = 0;

      }

      c->dg[k] += '0';

   }

   for (; i >= 0; --i, ++k) {

      c->dg[k] = a.dg[i] - borrow;

      if (c->dg[k] - '0' < 0) {

          c->dg[k] += 10;

          borrow = 1; 

      } else {

          borrow = 0;

      } 

   }

   c->dg[k] = '';

   c->length = k;

   //printf("%sn", c->dg);

   rm_tzeros(c);

   //printf("%sn", c->dg);

   reverse(c);

}

 

void add(BN a, BN b, BN *c) {

    int i, j, k;

    int carry = 0;

//    printf("%d:%dn", a.length, b.length);

    if (a.sign == b.sign) {

        c->sign = a.sign;

    } else {

        if (a.sign == MINUS) {

            a.sign = PLUS;    //tmp

            sub(b, a, c);

        } else {

            b.sign = PLUS;    //tmp

            sub(a, b, c);

        }

        return;

    }

    for (i = a.length-1, j = b.length-1, k = 0; i >= 0 && j >= 0; --i, --j, ++k) {

        c->dg[k] = a.dg[i] -'0' + b.dg[j] - '0' + carry;

        //printf("%d %d %c %cn",k, c->dg[k], a.dg[i], b.dg[j]);

        carry = c->dg[k]/10;

        c->dg[k] = c->dg[k]%10 + '0';

//        printf("%d %cn", k, (*c->dg[k]);

    }

    for (; i >= 0; --i, ++k){

        if (carry != 0) {

            c->dg[k] = carry + a.dg[i] - '0';

            carry = c->dg[k]/10;

            c->dg[k] = c->dg[k]%10 + '0';

        } else {

            c->dg[k] = a.dg[i];

        }

    }

    for (; j >= 0; --j, ++k){

        if (carry != 0) {

            c->dg[k] = carry + b.dg[j] - '0';

            carry = c->dg[k]/10;

            c->dg[k] = c->dg[k]%10 + '0';

        } else {

            c->dg[k] = b.dg[j];

        }

    }

    if (carry != 0) {

        c->dg[k++] = carry + '0';

    }

    c->dg[k] = '';

    c->length = k;

    //printf("%d %sn", c->length, c->dg);

    reverse(c);

    //printf("%d %sn", c->length, c->dg);

}

 

void lshift(BN *a, int num) {

    int i;

    if ((a->length == 1) && (a->dg[0] == '0')) {

        return;

    }

    for (i = 0; i < num; ++i) {

        a->dg[i+a->length] = '0';

    }

    a->length += num;

    a->dg[a->length] = '';

}

 

void mul(BN a, BN b, BN *c) {

   int i, j;

   

   init(c);

   for (i = a.length - 1; i >= 0; --i) {

       for (j = 0; j < a.dg[i]-'0'; ++j) {

          //printf("add %s to %sn", b.dg, c->dg);

          add(*c, b, c); 

       }

       lshift(&b, 1);

       //printf("%s %sn", c->dg, b.dg);

   }

   c->sign = a.sign*b.sign; 

}

 

void div(BN a, BN b, BN *c) {

    int i;

    BN tmp;

    init(c);

    init(&tmp);

    for (i = 0; i < a.length; ++i) {

        //printf("bf: %sn", tmp.dg);

        if (i > 0) {

           lshift(&tmp, 1);

        } else {

           tmp.length = 1;

        }

        //printf("af: %sn", tmp.dg);

        tmp.dg[tmp.length-1] = a.dg[i];

        c->dg[i] = '0';

        //printf("%d b: %d:%d:%s %d:%d:%sn", i, tmp.sign, tmp.length, tmp.dg, b.sign, b.length, b.dg);

        rm_lzeros(&tmp);

        while (cmp(tmp, b) >= 0) {

            //printf("sub: %s %sn", tmp.dg, b.dg);

            sub(tmp, b, &tmp);

            c->dg[i]++;

        } 

    }

    c->dg[i] = '';

    c->length = a.length;

    //printf("%d %sn", i, c->dg);

    rm_lzeros(c);

    //printf("%sn", c->dg);

}

 

void init(BN *a) {

    a->sign = PLUS;

    a->length = 1;

    a->dg[0] = '0';

    a->dg[1] = '';

}

 

int main(int argc, char** argv) {

    BN m, n, rv;

    if (argv[1][0] == '-') {

        m.sign = MINUS;

        sprintf(m.dg, "%s", &(argv[1][1]));

    } else {

        m.sign = PLUS;

        sprintf(m.dg, "%s", argv[1]);

    }

    m.length = strlen(m.dg);

    if (argv[2][0] == '-') {

        n.sign = MINUS;

        sprintf(n.dg, "%s", &(argv[2][1]));

    } else {

        n.sign = PLUS;

        sprintf(n.dg, "%s", argv[2]);

    }

    n.length = strlen(n.dg);

    printf("m:%c%sn", m.sign>0?' ':'-', m.dg);

    printf("n:%c%sn", n.sign>0?' ':'-', n.dg);

    add(m, n, &rv);

    printf("m + n = %c%sn", rv.sign>0?' ':'-', rv.dg);

    sub(m, n, &rv);

    printf("m - n = %c%sn", rv.sign>0?' ':'-', rv.dg);

    lshift(&rv, 1);

    printf("m << 1 = %c%sn", rv.sign>0?' ':'-', rv.dg);

    mul(m, n, &rv);

    printf("m * n = %c%sn", rv.sign>0?' ':'-', rv.dg);

    div(m, n, &rv);

    printf("m / n = %c%sn", rv.sign>0?' ':'-', rv.dg);

}

Save the code to bignumarray.c, and compile the code using the command below,

gcc -o bignumarray bignumarray.c

And here are some simple tests,

./bignumarray 34123432143214321 1342

m: 34123432143214321

n: 1342

m + n =  34123432143215663

m - n =  34123432143212979

m << 1 =  341234321432129790

m * n =  45793645936193618782

m / n =  25427296678997

 

./bignumarray 332 13424312432

m: 332

n: 13424312432

m + n =  13424312764

m - n = -13424312100

m << 1 = -134243121000

m * n =  4456871727424

m / n =  0

Digit Array Representation

The big integer is still defined as below,

typedef struct {

   int sign;

   int length;

   char dg[MAXDG];

} BN;

But now the dg holds the values for each digit. For example, given two inputs 12345 and 678, then the two big integers are a, b, with a.length = 5, a.sign = 1, and a.dg[0]=5, a.dg[1]=4, a.dg[2]=3, a.dg[3]=2, a.dg[4]=1; b.length=3, b.sign = 1, and b.dg[0]=8, b.dg[1]=7, b.dg[2]=6. Note that the digits are stored in a reverse manner because it is more convenient to manipulate the end of an array than the start of an array (e.g. truncate the array by 1).

The arithmetic operations are similar to the previous implementation. Below is the sample C code.

#include <stdio.h>

#include <string.h>

 

#define MAXDG 500

#define PLUS 1

#define MINUS -1

 

typedef struct {

    int sign;

    int length;

    char dg[MAXDG];

} BN;

 

 

void add(BN a, BN b, BN *c);

void mul(BN a, BN b, BN *c);

void div(BN a, BN b, BN *c);

void sub(BN a, BN b, BN *c);

void digit_shift(BN *a, int d);

int cmp(BN a, BN b);

void zero_justify(BN* a);

void init(BN *a);

void printbn(BN a);

 

void printbn(BN a) {

    int i;

    //printf("length: %d    ", a.length);

    if (a.sign == MINUS) {

        printf("-");

    }

    for (i = a.length-1; i >= 0; --i) {

        printf("%c", '0' + a.dg[i]);

    }

    printf("n");

}

 

void init(BN *a) {

    a->length = 1;

    memset(a->dg, 0x00, sizeof(char)*MAXDG);

    a->sign = PLUS;

}

 

void zero_justify(BN* a) {

    while ((a->length > 1) && (a->dg[a->length-1] == 0)) {

        (a->length)--;

    }

    //if -0, make it +0

    if ((a->length == 1) && (a->dg[0] == 0)) {

        a->sign = PLUS;

    }

}

 

int cmp(BN a, BN b) {

    if (a.sign == PLUS && b.sign == MINUS) {

        return 1;

    } else if (a.sign == MINUS && b.sign == PLUS) {

        return -1;

    } else {

        //a and b have same sign

        if (a.length > b.length) {

            return 1*a.sign;

        } else if (a.length < b.length) {

            return -1*a.sign;

        } else {

            //a and b have same sign, same length

            int i;

            for (i = a.length-1; i >= 0; --i) {

                if (a.dg[i] > b.dg[i]) {

                    return 1*a.sign;

                } else if (a.dg[i] < b.dg[i]) {

                    return -1*a.sign;

                } 

            }

            return 0; //equal

        }

    }

}

 

void digit_shift(BN *a, int d) {

    int i;

    if ((a->length == 1) && (a->dg[0] == 0)) {

        return;

    }

    for (i = a->length - 1; i >= 0; --i) {

        a->dg[i+d] = a->dg[i];

    }

    for (i = 0; i < d; ++i) {

        a->dg[i] = 0;

    }

    a->length = a->length + d;

}

 

void add(BN a, BN b, BN *c) {

    int i;

    int carry = 0;

    init(c);

    if (a.sign == b.sign) {

        c->sign = a.sign;

    } else {

        if (a.sign == MINUS) {

            a.sign = PLUS;

            sub(b, a, c);

        } else {

            b.sign = PLUS;

            sub(a, b, c);

        }

        return;

    }

    c->length = (a.length > b.length)?a.length+1:b.length+1;

    for (i = 0; i < c->length; ++i) {

        c->dg[i] = (char)(a.dg[i] + b.dg[i] + carry)%10;

        carry = (a.dg[i] + b.dg[i] + carry)/10;

    }

    zero_justify(c); 

}

 

void sub(BN a, BN b, BN *c) {

    int borrow;

    int v;

    int i;

    init(c);

    if ((a.sign == MINUS) || (b.sign == MINUS)) {

        b.sign = -1*b.sign;

        add(a, b, c);

        return;

    }

    if (cmp(a, b) < 0) {

        sub(b, a, c);

        c->sign = MINUS;

        return;

    }

    //a > b, and both a and b are +

    c->length = a.length > b.length?a.length:b.length;

    borrow = 0;

    for (i = 0; i < c->length; ++i) {

        c->dg[i] = (a.dg[i] - borrow - b.dg[i]);

        if (c->dg[i] < 0) {

            c->dg[i] += 10;

            borrow = 1;

        } else {

            borrow = 0;

        }

    }

    zero_justify(c);

}

 

void mul(BN a, BN b, BN *c) {

    BN row, tmp;

    int i,j;

    init(c);

    for (i = 0; i < b.length; ++i) {

        for (j = 0; j < b.dg[i]; ++j) {

            add(*c, a, c);

        }

        digit_shift(&a, 1);

    }

    zero_justify(c);

}

 

void div(BN a, BN b, BN *c) {

    int i;

    BN tmp;

    init(c);

    init(&tmp);

    c->sign = a.sign * b.sign;

    a.sign = PLUS;

    b.sign = PLUS;

    c->length = a.length;

    for (i = a.length - 1; i >= 0; --i) {

        digit_shift(&tmp, 1);

        //printf("%d:%dn", tmp.length, b.length);

        //printbn(tmp);

        tmp.dg[0] = a.dg[i];

        c->dg[i] = 0;

        while (cmp(tmp, b) >= 0) {

            (c->dg[i])++;

            //printbn(tmp);

            sub(tmp, b, &tmp);

            //printbn(tmp);

        }

        //printf("div: %dn", c->dg[i]);

    }

    zero_justify(c);

}

 

void readbn(char *s, BN *a) {

    int i, j;

    init(a);

//    printf("readbn: %dn", strlen(s));

    for (i = strlen(s)-1, j = 0; i >= 0; --i) {

        a->dg[j++] = s[i] - '0';

    }

    a->length = j;

}

 

int main(int argc, char** argv) {

    BN m, n, rv;

    if (argv[1][0] == '-') {

        m.sign = MINUS;

        readbn(&argv[1][1], &m);

    } else {

        m.sign = PLUS;

        readbn(argv[1], &m);

    }

    if (argv[2][0] == '-') {

        n.sign = MINUS;

        readbn(&argv[2][1], &n);

    } else {

        n.sign = PLUS;

        readbn(argv[2], &n);

    }

    printf("m:");

    printbn(m);

    printf("n:");

    printbn(n);

    add(m, n, &rv);

    printf("m + n = ");

    printbn(rv);

    sub(m, n, &rv);

    printf("m - n = ");

    printbn(rv);

    digit_shift(&rv, 1);

    printf("m << 1 = ");

    printbn(rv);

    mul(m, n, &rv);

    printf("m * n = ");

    printbn(rv);

    div(m, n, &rv);

    printf("m / n = ");

    printbn(rv);

}

Save the code to bignumarray2.c, and compile the code using the command below,

gcc -o bignumarray2 bignumarray2.c

And here are some simple tests,

./bignumarray2 34123432143214321 1342

m:34123432143214321

n:1342

m + n = 34123432143215663

m - n = 34123432143212979

m << 1 = 341234321432129790

m * n = 45793645936193618782

m / n = 25427296678997

 

./bignumarray2 332 13424312432

m:332

n:13424312432

m + n = 13424312764

m - n = -13424312100

m << 1 = -134243121000

m * n = 4456871727424

m / n = 0

Additional Notes

1. Dynamic array is necessary in order to support arbitrary length of big integers.

2. Using a single byte to represent a single digit is a huge waste. There’re better representations.

References:
1. Programming Challenges

Primality Testing — Part 4. Fermat’s Little Theorem and A Probabilistic Approach

Previous posts on primality testing are all trying to factor the input number. For a very large number with few hundreds of bits, it is difficult to compute in that manner.
This post introduce a probabilistic approach based on Fermat’s primality test.

Fermat’s Little Theorem

Fermat’s primality test is based on Fermat’s Little Theorem,

If p is prime, then for every 1 <= a < p,

[ a ^ ( p – 1) ] mod p = 1

Fermat’s Little Theorem suggests a method to test if a number is prime or not. For a big input number p, it is not desirable to test every a. Normally, we randomly pick up an a, and if it satisfy the equation above, we claim it is prime, otherwise, it is not.

However, Fermat’s little theorem is a necessary condition for a number to be prime, but not sufficient. In other words, every prime number satisfy fermat’s little theorem, but it is also possible to find an a for a composite number n where [ a ^ (p – 1) ] mod n = 1.

It can be proven that most values of a will fail the test for a composite number. And if we pick several random a and they all pass the test, the possibility of making a wrong decision is rare. This algorithm cannot guarantee that positive primality test is 100% correct, but it is fast the possibility of wrong is low enough for many practical use cases.

Implementation

Below is a simple C implementation based on the idea above,

#include <stdio.h>

#include <stdlib.h>

 

unsigned long modexp(int x, int y, int n) {

    unsigned long i;

    if (y == 0) {

        return 1%n;

    }

    if (y == 1) {

        return x%n;

    } 

    i = modexp(x, y/2, n);

    if (y&0x01) {    //odd

        //return (x*i*i)%n;

//        printf("%d:%lu:%lun", x, i, (x%n*i*i));

        return x%n*i%n*i%n;

    } else {         //even

        return i*i%n;

    }

}

 

int main(int argc, char **argv) {

    int n;

    int a;

    int i;

    n = atoi(argv[1]);

    srand(time(NULL));

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    for (i = 0; i < 20 && i < n; ++i) {

        a = rand()%100;          

        if (a == 0) {

            a = 1;

        }

        //printf("a=%dn", a);

    if (modexp(a, n-1, n) != 1) {

            printf("%d is not prime, %dn", n, a);

            break;

    } 

    }

    if (i == 20 || i == n) {

          printf("%d is primen", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

 

}

Note that the code uses the modular exponentiation algorithm. Due to the integer overflow, the largest input the program can accept is 65536 = 2^16.

Save the code to prim5.c and compile it using the command below,

gcc -o prim5 prim5.c

Below are some simpe tests,

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 65521

65521 is prime

time: 0:69

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 2

2 is prime

time: 0:53

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 65536

65536 is not prime, 8

time: 0:53

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 63703

63703 is prime

time: 0:69

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 61493

61493 is prime

time: 0:67

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 13

13 is prime

time: 0:54

References:
1. Wikipedia Primality test: http://en.wikipedia.org/wiki/Primality_test
2. Algorithms

How to Calculate Modular Exponentiation

Modular exponentiation is a exponentiation performed over modulus. It is used in many cryptographic algorithms.

Modular exponentiation accepts three input integers, say x, y and n and computes according to the formula below,

(x ^ y) mod n

Native Approaches

A naive approach to compute modular exponentiation is to compute the exponentiation first and then the modulus.

This approach is OK for small x, y and n. But in cryptographic algorithms, the values of x, y and n are normally several hundreds of bits long and this approach takes too much memory and multiplication.

A slightly improved version is to compute x mod n repeatedly and multiply to the current result. That is,

rv = x mod n => rv = rv * (x mod n) => rv = rv *(x mod n) ….

This approach reduces the usage of memory. But it still perform poorly when y is very big because the number of the multiplications it needs to perform is very large. Suppose y = 2^100, then the number of multiplications it needs to perform is also 2^100, which is almost infeasible in practice.

A Better Approach

A better approach is to repeatedly square x mod n. That is,

x mod n => (x mod n)^2 => (x mod n)^4 => (x mod n)^8 =>…=>

(x mod n) ^ logy

In this way, the number of multiplications needed can be reduced to logy. For y = 2^100, the number of multiplications needed is just log(2^100) = 100.

In general, this approach can be expressed by the equation below,

(x ^ y) mod n = [ (x ^ ( y / 2 ) ) mod n] ^ 2                                  if y is even

(x ^ y) mod n = [ x * [( x ^ ( y / 2 ) ) mod n ] ^ 2 ] mod n              if y is odd

This approach can be easily implemented in a recursive way. Below is a C implementation for this algorithm,

/**

polymial-time computation for (x^y)%N

the idea is for mod N

(x^y) = [x^(y/2)]^2 if y is even

        = x*[x^(y/2)]^2 if y is odd

**/

#include <stdio.h>

#include <stdlib.h>

 

unsigned long modexp(int x, int y, int n) {

    unsigned long i;

    if (y == 0) {

        return 1%n;

    }

    if (y == 1) {

        return x%n;

    } 

    i = modexp(x, y/2, n);

    printf("%d:%ldn", y, i);

    if (y&0x01) {    //odd

        //return (x*i*i)%n;

        return (x%n*i*i)%n;

    } else {         //even

        return (i*i)%n;

    }

}

 

int main(int argc, char **argv) {

    int x, y, N;

    x = atoi(argv[1]);

    y = atoi(argv[2]);

    N = atoi(argv[3]);

    printf("(x^y)modN = %lun", modexp(x, y, N));

    return 0;

}

One thing to note in the code above is that when y is odd, we compute using (x%n*i*i)%n instead of (x*i*i)%n. This is to prevent overflow caused by multiplication of x*i*i when x is a very big number.

Also note that when both x and n are big, the modulo x%n is big, the x%n*i*i can still cause overflow. Solving this issue requires different representations of integers and it is not discussed in this post. For most commonly used integers (except cryptographic), the program above is good enough.

For a given input number y, the number of recursions is ~logy. Therefore the time complexity is O(logy). When y is measured in bits, say, N bits, then the number of recursions is ~N. Therefore the time complexity is O(N).

Save the code as modexp.c and compile it using the command,

gcc -o modexp modexp.c

Below are some sample tests,

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp 4 13 497

3:4

6:64

13:120

(x^y)modN = 445

 

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp 1234 5678 90

2:64

5:46

11:64

22:64

44:46

88:46

177:46

354:64

709:46

1419:64

2839:64

5678:64

(x^y)modN = 46

 

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp 3 14 30

3:3

7:27

14:27

(x^y)modN = 9

 

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp 2 16 123

2:2

4:4

8:16

16:10

(x^y)modN = 100

An Alternative Approach

An alternative approach can help us to get rid of the recursion in the program above. The idea is based on the binary representations of an integer in computer. For example, in the first test case above, y = 13 = 2^3 + 2^2 + 2^0. Therefore to compute (4^13) mod 497, we can compute a = (4^(2^3)) mod 497, b = (4^(2^2)) mod 497 and c = (4^(2^0)) mod 497, and then multiply the results together.

In general, given x and n, we can pre-compute [ x ^ (2 ^ N)  ] mod n, where N is the number of bits needed to represent the biggest value of y. If y is constrained to be less than 2^100, then N = 100.

In addition, [x ^ (2 ^ N) ] mod n = { [x ^ ( 2 ^ (N-1) ) ] mod n} ^ 2. Therefore, we can compute from N = 0, 1 to N = 100 easily.

r0 = [x ^ (2 ^ 0)] mod n = x mod n
r1 = [x ^ (2 ^ 1)] mod n = r0 ^ 2
r2 = [x ^ (2 ^ 2)] mod n = r1 ^ 2

For a given y, we then scan the bits that is set to 1 and read the results from the pre-computation and multiply the results together.

Below is a sample C implementation for this approach,

/**

polymial-time computation for (x^y)%N

non-recursive version

**/

#include <stdio.h>

#include <stdlib.h>

 

#define MAXP 1000

unsigned long m[MAXP];

 

void precompute(int x, int n) {

   int i;

   m[0] = x%n;

   for (i = 1; i < MAXP; ++i) {

       m[i] = (m[i-1]*m[i-1])%n;

//       printf("%d:%ldn", i, m[i]);

   }

}

 

unsigned long modexp(int x, int y, int n) {

    int i;

    unsigned long rv;

    int tmp;

    rv = -1;

    for (i = 0; tmp != 0; ++i) {

        tmp = y >> i;

        if (tmp & 1) {

            if (rv == -1) {

                rv = m[i]%n;

            } else {

                rv = (m[i]*rv)%n;

            }

            printf("%d:%ldn", i, rv);

        }

    }

    return (rv == -1 ? 1%n:rv);

}

 

int main(int argc, char **argv) {

    int x, y, N;

    x = atoi(argv[1]);

    y = atoi(argv[2]);

    N = atoi(argv[3]);

    precompute(x, N);

    printf("(x^y)modN = %lun", modexp(x, y, N));

    return 0;

}

Save the code to modexp2.c and compile it using the command below,

gcc -o modexp2 modexp2.c

Below are some simple tests,

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp2 4 13 497

0:4

2:30

3:445

(x^y)modN = 445

 

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp2 1234 5678 90

1:46

2:46

3:46

5:46

9:46

10:46

12:46

(x^y)modN = 46

 

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp2 3 14 30

1:9

2:9

3:9

(x^y)modN = 9

 

roman10@ra-ubuntu-1:~/Desktop/integer$ ./modexp2 2 16 123

4:100

(x^y)modN = 100

References:
1. Algorithms. Dasgupta, C.H.Papadimitriou, and U.V. Vazirani
2. Wikipedia Modular Exponentiation: http://en.wikipedia.org/wiki/Modular_exponentiation
3. Numerical Recipes: http://numericalrecipes.blogspot.com/2009/03/modular-exponentiation.html

Primality Testing — Part 3. A Further Optimized Method

The previous two posts cover the naive approaches for primality testing and a method with precomputation. This post introduces a new method which could perform better than previous approaches and it can be used together with the precomputation.

The optimization is based on the observation that all prime numbers are in the form of 6k+1 or 6k+5, with 2, 3 and 5 being the only exceptions. This is because all integers can be expressed in 6k+i for some integer k, and i  is from the set {0, 1, 2, 3, 4, 5}. And 2 divides 6k+0, 6k+2 and 6k+4, 3 divides 6k+3.

Therefore, the primality test for a given integer n can be done by checking if n is divisible by 2 and 3 first, and then checking if n is divisible for all numbers of form 6k+1 or 6k+5 and less than sqrt(n).

The idea of this optimization is actually similar to the Sieve of Eratosthenes algorithm described in previous post. We eliminated testing divisibility for numbers that is divisible by the prime factors.

In general, all prime numbers are in the form gk + i, where i < k and it represents the numbers that are coprime to g. In the case above, g = 6 and the coprime of 6 are 1 and 5. If i and g are not coprime, then i must be divisible by some prime factors of g.

We go through another example, suppose g = 2*3*5 = 30, then all integers can be expressed in the form 30*k + i where i is in the set {0, 1, 2….29}. Now we mark the numbers that divisible by 2, 3 and 5 respectively. 30*k + {0, 2, 4, 6, 8, 10, 12, …28} are divisible by 2; 30*k + {3, 6, 9, 12, 15, 18, 21, 24, 27} are divisible by 3; and 30*k + {5, 10, 15, 20, 25} are divisible by 5. The rest of the numbers {1, 7, 11, 13, 17, 19, 23, 29} are not divisible by 2, 3 and 5 and are therefore coprimes of 30.

So to test if a number n is prime number, we first test if the number is divisible by 2, 3, and 5 first, and then check if the number is divisible by numbers in the form 30*k + {1, 7, 11, 13, 17, 19, 23, 29}  up to sqrt(n).

Below is an implementation of the optimization for primality testing when g is set to 6,

/**

generate prime numbers

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int rdiv(int n, int m) {

    if (n%m ==0) {

        printf("%d is divisible by %dn", n, m);

        return 0;

    }

    return 1;

}

 

int prime(int n) {

    int i, k, ed, m;

    int cpr[2] = {1, 5};

    //check if we can divide the prime factors first

    if (n == 2 || n == 3) {return 1;}

    if (rdiv(n, 2) == 0) {

        return 0;

    } else if (rdiv(n, 3) == 0) {

        return 0;

    }

    //printf("ifprime: %d", ifprime);

    //check up to sqrt(n) of the form 6k + {coprime of 6} 

    ed = sqrt(n) + 1;

    for (k = 0; ;++k) {

        for (i = 0; i < 2; ++i) {

            m = 6*k + cpr[i];

            if (m < 3) {continue;}

            if (m > ed) {

                return 1;

            } 

            if (n%m == 0) {

                printf("%d is divisible by %dn", n, m);

                return 0;

            }

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    int i;

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is prime numbern", n);

    } else {

        printf("%d is not prime numbern", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

    return 0;

}

Save the code to prim4.c and compile the code with command below,

gcc -o prim4 prim4.c –lm

Below are some simple tests:

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2

2 is prime number

time: 0:49

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 3

3 is prime number

time: 0:62

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 5

5 is prime number

time: 0:52

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2147483647

2147483647 is prime number

time: 0:267

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2147483641

2147483641 is divisible by 2699

2147483641 is not prime number

time: 0:72

Note that the code can be further improved by adding the pre-computation we mentioned in part2.

Primality Testing — Part 2. Pre-computation with Sieve of Eratosthenes

Previous post covers two naive methods for primality testing. This post introduces an algorithm to generate prime numbers called Sieve of Eratosthenes. With this algorithm, an optimization can be applied to speed up the primality testing.

Sieve of Eratosthenes

Suppose we want to generate prime numbers up to N. The algorithm consists of two loops, and it can be described as below,

  • Create a mask for all integers up to N.
  • In the first loop, let x be 2, mark 2x, 3x, 4x for all n such that nx < N.
  • finding an integer bigger than x and not marked, let it be the new x, and repeat the step above
  • the first loop terminates when we cannot find an integer bigger than x and not marked
  • In the second loop, starting from 2, collects the numbers that are not marked, which are the prime numbers up to N.

Below is an implementation of the Sieve of Eratosthenes algorithm,

/**

generate prime numbers

**/

#include <stdio.h>

#include <stdlib.h>

 

#define MAXV 10000000

#define ININ 1000

 

unsigned char bits[MAXV/8];

 

void setbit(int pos) {

    bits[pos/8] |= (0x01 << (pos%8));

}

 

int testbit(int pos) {

    return (bits[pos/8] & (0x01 << (pos%8)));

}

 

void genprime(int n, int **pr, int *npr) {

   int i, j, k;

   for (i = 2; i < n; ++i) {

       if (!testbit(i)) {

           for (j = 2; i*j < n; ++j) {

               setbit(i*j);

           }

       }

   } 

   k = ININ;

   *pr = malloc(sizeof(int) * k);

   *npr = 0;

   for (i = 2, j = 0; i < n; ++i) {

       if (!testbit(i)) {

           (*pr)[j++] = i;

           (*npr)++;

           if (j == k) {

               k += ININ;

               *pr = realloc(*pr, sizeof(int)*k);

           }

       }

   }

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    int *pr;

    int npr;

    int i;

    struct timeval stTime, edTime;

    if (n > MAXV) {

        printf("the max n is %dn", MAXV);

        return 0;

    }

    gettimeofday(&stTime, NULL);

    genprime(n, &pr, &npr);

    printf("No. of prime numbers: %dn", npr);

    for (i = 0; i < npr; ++i) {

        printf("%d  ", pr[i]);

    }

    printf("n");

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

    if (pr != NULL) {

        free(pr);

    }

    return 0;

}

Note that the code above use a single bit to indicate if an integer is marked or not.

Save the code to prim2.c, and then compile the code with the command below,

gcc -o prim2 prim2.c

Below is a screenshot of some simple tests,


Figure 1. Simple Tests of Implementation for Sieve of Eratosthenes

Primality Testing with Sieve of Eratosthenes for Precomputation

One way to speed up primality testing is to precompute a list of prime numbers up to N. Then we first check if an input integer n is divisible by the prime numbers generated. If the number is not divisible by any of those prime numbers, we then proceed to the test mentioned in the naive approach post.

Below is an implementation of this precomputation approach,

/**

a primality test with precomputation

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

#define MAXV 10000000

#define ININ 1000

 

unsigned char bits[MAXV/8];

 

void setbit(int pos) {

    bits[pos/8] |= (0x01 << (pos%8));

}

 

int testbit(int pos) {

    return (bits[pos/8] & (0x01 << (pos%8)));

}

 

void genprime(int n, int **pr, int *npr) {

   int i, j, k;

   for (i = 2; i < n; ++i) {

       if (!testbit(i)) {

           for (j = 2; i*j < n; ++j) {

               setbit(i*j);

           }

       }

   } 

   k = ININ;

   *pr = malloc(sizeof(int) * k);

   *npr = 0;

   for (i = 2, j = 0; i < n; ++i) {

       if (!testbit(i)) {

           (*pr)[j++] = i;

           (*npr)++;

           if (j == k) {

               k += ININ;

               *pr = realloc(*pr, sizeof(int)*k);

           }

       }

   }

}

 

int prime(int n, int *pr, int npr, int x) {

    int i;

    int j;

    

    for (i = 0; i < npr; ++i) {

        if (n == pr[i]) {

            return 1;

        }

        if (n%pr[i] == 0) {

            printf("%d is divisble by %dn", n, pr[i]);

            return 0;

        }

    }

    if (n == 2) {

        return 1;

    }

    if (!(n&1)) {

        printf("%d is divisible by 2n", n);

        return 0;

    }

    j = sqrt(n);

    //to be safe, run to j+1

    i = (x+1) > 3 ? (x+1):3;

    if (!(i&1)) {

        ++i;   //make sure i starts as odd number

    }

    for (; i <= j + 1; i+=2) {

        if (n%i == 0) {

            printf("%d is divisible by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int m, n;

    int *pr;

    int npr;

    int i;

    struct timeval stTime, stTime2, edTime;

    if (argc < 3) {

        printf("./prim3 <primality test number> <prime number generation limit>n");

        return 0;

    }

    m = atoi(argv[1]);

    n = atoi(argv[2]);

    if (n > MAXV) {

        printf("the max n is %dn", MAXV);

        return 0;

    }

    gettimeofday(&stTime, NULL);

    genprime(n, &pr, &npr);

    gettimeofday(&stTime2, NULL);

    //printf("No. of prime numbers: %dn", npr);

/*    for (i = 0; i < npr; ++i) {

        printf("%dn", pr[i]);

    }*/

    if (prime(m, pr, npr, n)) {

        printf("%d is prime numbern", m);

    } else {

        printf("%d is not prime numbern", m);

    }

    gettimeofday(&edTime, NULL);

    printf("total time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

    printf("test prime time: %u:%un", (unsigned int)(edTime.tv_sec - stTime2.tv_sec), (unsigned int)(edTime.tv_usec - stTime2.tv_usec)); 

    if (pr != NULL) {

        free(pr);

    }

    return 0;

}

Save the source code to prim3.c, and then compile the code with the command below,

gcc -o prim3 prim3.c –lm

Below are some simple tests for the compiled program,

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2 10

2 is prime number

total time: 0:144

test prime time: 0:48

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 3 10

3 is prime number

total time: 0:60

test prime time: 0:22

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483647 1000

2147483647 is prime number

total time: 0:469

test prime time: 0:302

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483643 1000

2147483643 is divisble by 3

2147483643 is not prime number

total time: 0:223

test prime time: 0:57

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483641 1000

2147483641 is divisible by 2699

2147483641 is not prime number

total time: 0:250

test prime time: 0:80

References:
1. Wikipedia: http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Primality Testing–Part 1. Naïve Approaches

Primality refers to determining if a given input number is prime or not. This post covers two naive approaches for primality testing.

Naive Approach 1

This most straightforward approach is to check if a given input integer n is divisible by integers from 2 to n-1. If a number is prime, then it should not divisible by any numbers in the range [2, n-1]. This is implemented as below,

/**

a naive method of determining if n is prime or not

divide from 2 to n-1

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int prime(int n) {

    int i;

    for (i = 2; i <= n-1; i++) {

        if (n%i == 0) {

            printf("%d is divisble by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is primen", n);

    } else {

        printf("%d is not primen", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Save the code to prim0.c. Compile the code with the command below,

gcc -o prim0 prim0.c –lm

Below are some sample testing results,

./prim0 2147483647

2147483647 is prime

time: 7:145569

 

./prim0 2147483640

2147483640 is divisble by 2

2147483640 is not prime

time: 0:46

 

./prim0 2

2 is prime

time: 0:51

 

./prim0 3

3 is prime

time: 0:49

Naive Approach 2: Improved Version of NA 1

The above algorithm can be improved. Firstly, if a number is divisible by an even integer, it must be divisible by 2. So we only need to check if the input number n is divisible by 2 and skip the rest of the even integers.

Secondly, if a number n is divisible by a number x, where x > sqrt(n). Then n = x*y, where y < sqrt(n), and y is also divisible by n. Therefore, we only need to check up to integer sqrt(n).

Below is the implementation with the two optimizations described,

/**

a naive method of determining if n is prime or not

divide from 2 to sqrt(n), skip the even numbers

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int prime(int n) {

    int i, j;

    if (n == 2) {

        return 1;

    }

    if (!(n&1)) {

        //n is even

        printf("%d is divisible by 2n", n);

        return 0;

    }

    j = sqrt(n);

    //to be safe, run to j+1

    for (i = 3; i <= j + 1; i+=2) {

        if (n%i == 0) {

            printf("%d is divisible by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is primen", n);

    } else {

        printf("%d is not primen", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Note that sqrt is a floating point operation. We check until sqrt(n) + 1, just to be safe.

Save the code to prim1.c, and compile the code using the command below,

gcc -o prim1 prim1.c –lm

Below are some testing results,

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2

2 is prime

time: 0:50

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 3

3 is prime

time: 0:60

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2147483647

2147483647 is prime

time: 0:308

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2147483641

2147483641 is divisible by 2699

2147483641 is not prime

time: 0:73

References:
1. Wikipedia Primality Test: http://en.wikipedia.org/wiki/Primality_test
2. Algorithms by Dasgupta, CH Papadimitriou and UV Vazirani

Computing Fibonacci Numbers

Fibonacci series is a number sequence of 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 44, … It can be expressed by the formula below,

1

Fibonacci series grows almost as fast as the power of 2. Fn ~= 2^(0.694n).

Computation of Fibonacci Numbers
The naive approach to compute Fibonacci numbers is to recursion according to the formula of Fabonacci numbers. A C implementation is given below,

/*

a naive approach of computing fibonacci series, the computation time also satisfy fabonacci series, since

fabonacci series grows almost as fast as power of 2, the computation is expoential time

*/

#include <stdio.h>

#include <stdlib.h>

 

unsigned long fibo(int n) {

    if (n == 0) {

       return 0;

    } else if (n == 1) {

       return 1;

    } else {

       return fibo(n-1) + fibo(n-2);

    }

}

 

int main(int argc, char **argv) {

    int n;

    unsigned long res;

    struct timeval stTime, edTime;

    n = atoi(argv[1]);

    gettimeofday(&stTime, NULL);

    res = fibo(n);

    gettimeofday(&edTime, NULL);

    printf("%lu, time takes (%u:%u)n", res, (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Suppose the time takes to compute fibo(n-1) and fibo(n-2) are O(n-1) and O(n-2) respectively, then the time complexity of computing fibo(n) is O(n) = O(n-1) + O(n-2). The time complexity itself follows the Fibonacci series. Since we know Fn~= 2^(0.694n), then O(n) ~=2^(0.694n). The computation is exponential. And since O(n)~=(2^0.694)^n ~=1.62^n, compute the Fn+1 takes around 1.6 times longer than computing Fn.

Compile the code using the command,

gcc -o fibo fibo.c

Below is a sample run of the program,

Figure 1. Execution of Recursive Version Program

It is not easy to observe that a Fibonacci number depends on its two previous numbers. The computation above computes some small Fibonacci numbers for many times. We can use an array to remember the previous computation. A C implementation is shown as below,

/*

use dynamic programming to remember previous computation results

this computation runs in linear time

*/

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

int main(int argc, char **argv) {

    int n;

    int i;

    struct timeval stTime, edTime;

    unsigned long* fibo;

    n = atoi(argv[1]);

    fibo = (unsigned long*) malloc((n+1)*sizeof(unsigned long));

    memset(fibo, 0x00, sizeof(unsigned long)*n);

    fibo[1] = 1;

    gettimeofday(&stTime, NULL);

    for (i = 2; i <= n; ++i) {

        fibo[i] = fibo[i-1] + fibo[i-2];

    }

    gettimeofday(&edTime, NULL);

    printf("%lu, time takes (%u:%u)n", fibo[n], (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

}

Compile the code using the command below,

gcc -o fibo1 fibo1.c

A sample run is as below,

Figure 2. Execution of Dynamic Programming Version

It is not difficult to tell the code above runs in linear time. However, the code takes a lot of memory space and if we’re to compute a Fibonacci number Fn for a very large n, then we may run out of memory space.
To compute Fn, we only need values for Fn-1 and Fn-2, this leads to the code below,

/*

use dynamic programming to remember previous computation results

this computation runs in linear time

*/

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

int main(int argc, char **argv) {

    int n;

    int i;

    struct timeval stTime, edTime;

    unsigned long fibo[2];

    n = atoi(argv[1]);

    fibo[0] = 0;

    fibo[1] = 1;

    gettimeofday(&stTime, NULL);

    for (i = 2; i <= n; ++i) {

        fibo[i%2] = fibo[0] + fibo[1];

    }

    gettimeofday(&edTime, NULL);

    if (n&0x01) {

        printf("%lu, time takes (%u:%u)n", fibo[1], (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

    } else {

        printf("%lu, time takes (%u:%u)n", fibo[0], (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec));

    }

}

Compile the code using the command below,

gcc -o fibo2 fibo2.c

A sample run is as below,

Figure 3. Execution of Dynamic Programming Version Improved

References:
1. Algorithms. Dasgupta, C.H. Papadimitriou, and U. V. Vazirani, 2006
2. Wikipedia Fibonacci Number: http://en.wikipedia.org/wiki/Fibonacci_number

Million Integer Problems — Part 2. Sorting a Million Integers

This is the second post for million integer problems. The first post is about generating million integers. You can find it here. This post deals with sorting of million integers.

It can be proven that the best that a comparison-based sorting can do has average run time of nlogn. However, given specific inputs, there’re algorithms that can do better than nlogn. The million integer sorting is one of such cases.

1. Sorting a million integers in the range of [0, 1 000 000]. Numbers can appear many times.

Since the comparison-based sorting algorithms can only do nlogn, we’ll see if we can sort without compare the one million elements.

The inputs are all integers, and it’s within a specific range [0, 1 000 000]. Integers can be the index of arrays. If we have an array with 1 million indices, every time an input appears, we mark the element in the array. If an integer appear N times, we mark the array element N times.

After reading all inputs, the indices corresponding to all input integers are marked. Because indices are ordered, we simply read from array element 0 to 1 million, if the array element is marked as N, we output the index N times.

The C implementation is as below,

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

#define MAXV 1000000

 

unsigned int buckets[MAXV+1];

 

int main(int argc, char **argv) {

   FILE* inputF, *outputF;

   int i;

   char aline[15];

   inputF = fopen("input1.txt", "r");

   if (inputF == NULL) {

       printf("cannot open input filen");

       return 0;

   }

   memset(buckets, 0x00, MAXV+1);

   while (fgets(aline, 15, inputF) != NULL) {

       i = atoi(aline);

//       printf("...%dn", i);

       ++buckets[i];

   }

   outputF = fopen("output1.txt", "w");

   for (i = 0; i <= MAXV; ++i) {

       while (buckets[i]--) {

           fprintf(outputF, "%dn", i);

       }

   }

   fclose(inputF);

   fclose(outputF);

}

You can use the code in part 1 generate 1 million integers to generate the input files for the above code.

2. Sorting a million integers in the range of [0, 1 000 000]. No number appears twice.

Compare with first problem, there’s one additional condition, the input integers are unique. This means the array element we used in first problem has either value 0 or 1. A single bit is enough to represent the array element.

The C implementation is as below,

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

#define MAXV 1000000

 

unsigned char buckets[(MAXV+7)/8];

 

void setb(int i) {

    int idx = i/8;

    int bpos = i%8;

    buckets[idx] |= (0x01 << bpos);

}

 

int checkb(int i) {

    int idx = i/8;

    int bpos = i%8;

    return (buckets[idx] & (0x01 << bpos));

}

 

int main(int argc, char **argv) {

   FILE* inputF, *outputF;

   int i;

   char aline[15];

   inputF = fopen("input.txt", "r");

   if (inputF == NULL) {

       printf("cannot open input filen");

       return 0;

   }

   memset(buckets, 0x00, (MAXV+7)/8);

   while (fgets(aline, 15, inputF) != NULL) {

       i = atoi(aline);

//       printf("...%dn", i);

       setb(i);

   }

   outputF = fopen("output.txt", "w");

   for (i = 0; i <= MAXV; ++i) {

       if (checkb(i)) {

           fprintf(outputF, "%dn", i);

       }

   }

   fclose(inputF);

   fclose(outputF);

}

You can use the code in part 1 generate 1 million integers to generate the input files for the above code.

2.1 Implementation from Programming Perls

The classic book Programming Perls has a implementation that sorts the input integers. The code is as below,

/* Copyright (C) 1999 Lucent Technologies */

/* From 'Programming Pearls' by Jon Bentley */

 

/* bitsort.c -- bitmap sort from Column 1

 *   Sort distinct integers in the range [0..N-1]

 */

 

#include <stdio.h>

 

#define BITSPERWORD 32

#define SHIFT 5

#define MASK 0x1F

#define N 10000000

int a[1 + N/BITSPERWORD];

 

void set(int i) {        a[i>>SHIFT] |=  (1<<(i & MASK)); }

void clr(int i) {        a[i>>SHIFT] &= ~(1<<(i & MASK)); }

int  test(int i){ return a[i>>SHIFT] &   (1<<(i & MASK)); }

 

int main()

{    int i;

    for (i = 0; i < N; i++)

        clr(i);

/*    Replace above 2 lines with below 3 for word-parallel init

    int top = 1 + N/BITSPERWORD;

    for (i = 0; i < top; i++)

        a[i] = 0;

 */

    while (scanf("%d", &i) != EOF)

        set(i);

    for (i = 0; i < N; i++)

        if (test(i))

            printf("%dn", i);

    return 0;

}

3. Related problems
3.1 There are 1 000 001 integers in the range of [0, 1 000 000]. Given a million unique integers in the range of [0, 1 000 000], find the integer that is missing.

This is actually the same as integer sorting problems. Once you sort the integers with the array, you can iterate through the array to find the array index whose element is not marked.

There are other problems. If you happens to encounter one, you’re welcome to add it in the comments. 🙂

References

1. Programming Perls. Second Edition.

2. Programming Perls Source Code website: http://www.cs.bell-labs.com/cm/cs/pearls/code.html

Million Integer Problems — Part 1. Generate 1 million integers

This is a post for a set of problems that deals with millions of integers. Part 1 deals with how to generate millions of integers.

This is normally used as a program to generate inputs for testing subsequent problems. We’ll look at two problems with a slight difference.

1. Allow user to input m, n, where m is the number of integers to generate, and n is the maximum value allowed. There can be repeated integers. m can be at any range of [0, 1,000,000].

This is simple. We can simply use a random number generate function to generate the required number of integers. To ensure them not exceed the max value, we use %(max+1).

The C code is as below,

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

 

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    int mx = atoi(argv[2]);

    int i;

    printf("%d...n", RAND_MAX);

 

    srand(time(NULL));

    //the problem is equivalent to reorder integers

    for (i = 0; i < n; ++i) {

        printf("%dn", rand()%(mx+1));

    }

    return 0;

}

The code accepts two input parameters. The first indicates the number of integers to generate and the second indicates the max value of the generated numbers. The code use the current time as the random number generator seed, and then generate n number of integers in the range of [0, mx].

Note that for the RAND_MAX defines the maximum value returned by rand(). The value of RAND_MAX differs on implementations, but it is at least 32767.

2. Same as 1, except that all integers has to be unique.

The first idea is to use a bitmap or integer of arrays to mark if an integer has been generated or not, and if it has been seen before, we try generate another integer.

This is easy and straightforward, but as the number of integers increases, soon we’ll find that we’re generating lots of integers that have been seen before. Generating a new unique integer becomes more and more difficult.

The better way is to think it as a integer shuffle problem. It’s pretty easy to generate a ordered set of integers within the range [0, mx], and what left to do is to shuffle the integers.

This can be done with the code below,

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

 

#define MAXV 2000000

 

unsigned int buckets[MAXV+1];

int main(int argc, char **argv) {

    int n = atoi(argv[1]);

    int mx = atoi(argv[2]);

    int i;

    int j;

    int tmp;

    printf("%d...n", RAND_MAX);

 

    for (i = 0; i < mx+1; ++i) {

        buckets[i] = i;

    }

 

    srand(time(NULL));

    //the problem is equivalent to reorder integers

    for (i = 0; i < n; ++i) {

        j = i + rand()%(mx+1-i);         //make sure the generated number in the range [i, mx]

        //swap the numbers at j and i

        tmp = buckets[j];

        buckets[j] = buckets[i];

        buckets[i] = tmp;

        printf("%dn", buckets[i]);

    }

    return 0;

}

The code takes two user inputs, including the number of integers to generate n, and the max possible value mx.

The code first generate an ordered set of integers within [0, mx] and store them in an array. In the second loop, each time we generate a random index and ensures it’s in the range of [i, mx], where i is the index of the integer we’re trying to output. We then swap the value of the current index i with the value of the randomly generated index. The swapped value at current index i is the integer for output.

In this way, every time we select an integer randomly from the set of integers not seem yet.

Code from Programming Perls

The classic book Programming Perls has a code that generates the integers. The code is as below,

/* Copyright (C) 1999 Lucent Technologies */

/* From 'Programming Pearls' by Jon Bentley */

 

/* bitsortgen.c -- gen $1 distinct integers from U[0,$2) */

 

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

#define MAXN 2000000

int x[MAXN];

int randint(int a, int b)

{    return a + (RAND_MAX * rand() + rand()) % (b + 1 - a);

}

 

int main(int argc, char *argv[])

{    int i, k, n, t, p;

    srand((unsigned) time(NULL));

    k = atoi(argv[1]);

    n = atoi(argv[2]);

    for (i = 0; i < n; i++)

        x[i] = i;

    for (i = 0; i < k; i++) {

        p = randint(i, n-1);

        t = x[p]; x[p] = x[i]; x[i] = t;

        printf("%dn", x[i]);

    }

    return 0;

}

The code crashes on my machine. Because the RAND_MAX value on my machine is 2147483647 and the method randint() will cause integer overflow.

References:
1. Programming Perls. Second Edition.

Euclid’s algorithm

Side Note: This is first written on Oct 17 2007.

Euclid’s algorithm is used to find the greatest common divisor for two positive integers, let’s say m and n. It’s the first algorithm introduced in Knuth’s classic book, <<The Art of Computer Programming >> . The following text are retrieved from the book,


a. [Find remainder] Divide m by n and let r be the remainder. (0<=r and r< n). b. [Is it zero] If r = 0, terminates the program; n is the answer. c. [Reduce] Set m = n, n = r, and go back to step a.

A C/C++ implementation of the above algorithm is as follows,

//please ensure that variables passed to the function is positive 
int euclid(int m, int n) {
	int r = m%n;
	while(r!=0) {
	    m = n;
	    n = r;
	    r = m%n;
        }
	return n;
}

This algorithm/program can also be used to decide whether two positive integers are relative prime(no common divisor except one) or not. All we need to do is to check the return value is 1 or not.