Prime number test

MillerRabinAlgorithm.cpp

//
//  Miller-Rabin Algorithm
//
//  Created by Hao Ling on 2014-03-19.
//  Copyright (c) 2014 Hao Ling. All rights reserved.
//

#include <iostream>
#include "math.h"

using namespace std;

/** getBitNum() ********************************************************
 * @params: long long x: input number
 * @return: the number of bits
 * @descrip: Get the number of bits of a given input number.
 ***********************************************************************/
int getBitNum
(long long x){
    // Record the number of bits
    int i = 1;
    // Count the bit number
    while ( x/2 != 0 ) {
        i++;
        x /= 2;
    }
    return i;
}

/** getBits() **********************************************************
 * @params: long long b: input number
            bool bits[]: binary representation
            int length: the number of bits
 * @return: none
 * @descrip: Convert a number to its binary representation.
 ***********************************************************************/
void getBits
(long long b, bool bits[], int length){
    // Convert a number to its binary representation
    for (int i=0; i<length; i++) {
        bits[i] = b % 2;
        b /= 2;
    }
}

/** sm() ***************************************************************
 * @params: long long x, long long b, long long n:
            parameters in pow(x,b) mod n
 * @return: the formula result
 * @descrip: Square and multiply algorithm.
 ***********************************************************************/
long long sm
(long long x, long long b, long long n){
    
    long long result = 1;
    
    // b is an l-bit binary number
    int l = getBitNum(b);
    bool bits[l];
    getBits(b, bits, l);
    
    for (int i=l-1; i>=0; i--) {
        
        // Square
        result = (result*result) % n;
        
        // Multiply
        if ( bits[i] == 1 ) {
            result *= x;
            result = result % n;
        }
    }
    
    return result;
}

/** mrTest() ***********************************************************
 * @params: long long m: the number to test
 * @return: Returns 1 if composite; returns 0 if inconclusive.
 * @descrip: Test if a number is composite or inconclusive.
 ***********************************************************************/
bool mrTest(long long m){
    
    /*
     * Step 1: Find integers k and z with k>0 and z is odd
     * so that m-1 = pow(2,k)*z.
     */
    int k = 0;
    long long z;
    long long temp = m-1;
    
    while ( temp%2 == 0 ) {
        k++;
        temp /= 2;
    }
    z = temp;
    
    /*
     * Step 2: Run the test.
     */
    long long g;
    
    // Select a random integer c so that 1<c<(m-1)
    long long c;
    c = rand() % (m-2) + 2;
    
    g = sm(c, z, m);
    
    if ( g%m == 1 ) {
        // Return inconclusive
        return 0;
    }
    
    for (int i=0; i<k; i++) {
        
        if ( g%m == (m-1) ) {
            // Return inconclusive
            return 0;
        }
        
        g = g*g;
        
    }

    // Return composite
    return 1;
}

/** main() *************************************************************
 * @descrip: The main function.
 ***********************************************************************/
int main(int argc, const char * argv[])
{
    
    // Define the number to test
    long long m = 849330323;
    // Define the maximum probability
    // that a composite number is not conclusively determinded as composite
    double maxProb = 0.0001;
    
    // Define the minimum number of tests
    int t = 1;
    double prob = 0.25;
    // Get the minimum number of tests
    while ( prob > maxProb ) {
        t++;
        prob *= 0.25;
    }
    
    cout << "Number of test needed is "<<t<<endl;
    
    for (int i=0; i<t; i++) {
        
        if ( mrTest(m) == 1 ) {
            cout << m << " is composite." << endl;
            return 1;
        }
        
    }
    
    cout << m << " is prime with the probability of " << (1-maxProb)*100 << "%." << endl;
    
    return 0;
}