algorithm-examples/cpp-algorithm/src/math/miller_rabin.h at main · codejsha/algorithm-examples · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#ifndef CPP_ALGORITHM_MILLER_RABIN_H
#define CPP_ALGORITHM_MILLER_RABIN_H

#include <random>

namespace MillerRabin
{
    /**
     * \brief Compute the modular exponentiation.
     * \param base base
     * \param exponent exponent
     * \param mod modulus
     * \return result: (base^exponent) % mod
     */
    int ModularExponentiation(int base, int exponent, int mod);

    /**
     * \brief Test whether the number is prime.
     * Witness the Fermat's theorem: a^(p-1) = 1 (mod p)
     * \param exponent exponent
     * \param number number
     * \return whether the number is prime
     */
    bool Witness(int exponent, int number);

    /**
     * \brief Miller-Rabin primality test.
     * \param number input number
     * \param repeats repeat times to increase the probability of correctness
     * \return whether the number is prime or not
     */
    bool MillerRabinPrimalityTest(int number, int repeats);
}

// ----------------------------------------------------------------------------
inline int MillerRabin::ModularExponentiation(int base, int exponent, const int mod)
{
    int result = 1;
    while (exponent > 0)
    {
        if (exponent % 2 == 1)
        {
            result = (result * base) % mod;
        }
        exponent /= 2;
        base = (base * base) % mod;
    }
    return result;
}

// ----------------------------------------------------------------------------
inline bool MillerRabin::Witness(int exponent, const int number)
{
    std::random_device rd;
    std::mt19937 generator(rd());
    std::uniform_int_distribution distribution(2, number - 1);
    const int random = distribution(generator);
    int x = ModularExponentiation(random, exponent, number);

    if (x == 1 || x == number - 1)
    {
        return true;
    }

    while (exponent != number - 1)
    {
        x = (x * x) % number;
        exponent *= 2;

        if (x == 1)
        {
            return false;
        }
        if (x == number - 1)
        {
            return true;
        }
    }

    return false;
}

// ----------------------------------------------------------------------------
inline bool MillerRabin::MillerRabinPrimalityTest(const int number, const int repeats)
{
    if (number < 2)
    {
        return false;
    }
    if (number == 2)
    {
        return true;
    }
    if (number % 2 == 0)
    {
        return false;
    }

    int d = number - 1;
    while (d % 2 == 0)
    {
        d /= 2;
    }

    for (int i = 0; i < repeats; ++i)
    {
        if (!Witness(d, number))
        {
            return false;
        }
    }

    return true;
}

#endif