prime number testing

unsigned int mulmod(unsigned int a, unsigned int b, unsigned int m);
int isPrime(unsigned int n)
{
    /* given `n' determine whether it is prime.
     * return 0 => composite. 1 => prime.
     *
     * method:
     * when n is odd > 7, determine if it is a strong pseudoprime
     * to bases 2, 3, 5 & 7. if not then we are composite.
     * otherwise we are prime because n < 25x10^9 except for one
     * exception, specially tested.
     */

    unsigned int r, s;
    int i, j;
    unsigned int a;

    if (n < 2) return 0;
    if (n == 2) return 1;

    /* only spsp < 25x10^9 with bases 2, 3, 5 & 7 */
    if (n == 3215031751) return 0;

    s = n - 1;
    r = 0;
    for (;;) {
        unsigned int c = s/2;
        if (c*2 == s) {
            ++r;
            s = c;
        }
        else break;
    }

    if (!r) return 0;
    --r;

    a = 0;
    for (;;) {
        unsigned int b;
        unsigned int y;
        unsigned int z;

        a += 2;

        if (a == 9) break;
        else if (a == 4) a = 3;        
        if (n == a) break;

        b = 1;
        z = a;
        y = s;
        for (;;) {
            unsigned int c = y/2;
            if (c*2 != y) {
                b = mulmod(z, b, n);
            }
            if (!c) break;
            y = c;
            z = mulmod(z, z, n);
        }

        if (b == 1 || b == n-1) continue;

        for (j = 0; j < r; ++j) {
            b = mulmod(b, b, n);
            if (b == n-1) break;
        }
        if (j == r) return 0;
    }
    return 1;
}

unsigned int mulmod(unsigned int a, unsigned int b, unsigned int m)
{
    /* return a*b (mod m).
     * take special care of unsigned quantities and overflow.
     */

    unsigned int z;
    unsigned int y = a;
    int t;
    z = 0;

    for (;;) {
        unsigned int c = b/2;
        if (c*2 != b) {
            z += y;
            if (z < y || z >= m) z -=m;
        }

        b = c;
        if (!b) break;

        t = y;
        y <<= 1;
        if (t < 0 || y >= m) y -= m;
    }
    return z;
}