[GO] [Internal_math version (2)] Decoding the AtCoder Library ~ Implementation in Python ~

0. Introduction

AtCoder Official Algorithm Collection AtCoder Library (** ACL **) released on September 7, 2020 it was done. I thought it was a good opportunity because most of the algorithms included in ACL were new to me, so I did everything from studying algorithms to implementing them in Python.

In this article we will look at ** internal_math **.

internal_math is an assortment of number-theoretic algorithms used inside ACLs and has the following contents.

Of these, in this article

• is_prime
• primitive_root

Handles. I will not touch on constexpr (constant expression) itself.

Not covered in this article

• safe_mod
• pow_mod
• inv_gcd
• barrett

Is dealt with in the following article. Please see that if you like. [[Internal_math version ①] Decoding the AtCoder Library ~ Implementation in Python ~] qiita_internal_math_1

Target readers

• If you look at the ACL code but ** don't know what you're doing **.
• Those who do not understand C ++ and want to read it in Python.

Untargeted readers

• Those who want Python-optimized code for ACLs. → Since it was implemented to be the same as ACL as much as possible, execution speed in Python is not considered at all. There seems to be a movement to make it available directly from Cython, so it may be good to follow that.

Referenced

This is a wikipedia page related to is_prime.

• [Japanese Wikipedia --Fermat's Little Theorem] wiki_fermat
• [Japanese Wikipedia-Probable Prime] wiki_prp

This is a paper on strong pseudoprime.

• GERHARD JAESCHKE (1993) "ON STRONG PSEUDOPRIMES TO SEVERAL BASES"

Miller-Rabin This is an article by @ srtk86 about the primality test. It's easy to understand.

• [Various prime judgments-Fermat test mirror rabin prime judgment method] qiita_miller_rabin

Describes primitive roots.

• [Contributors to Wikibooks, "Primitive Number Theory / Primitive Roots and Indexes," Wikibooks, Free Textbooks] wikibooks_primitive_root

1. is_prime

Consider determining if the positive integer $ n $ is a prime number.

1.1. Definitive primality test

The definition of a prime number is "a natural number greater than 1 with only one positive divisor and itself", so divide $ n $ by a natural number from 2 to $ n -1 $, and if it is not divisible, $ n $ Can be said to be a prime number. This is the so-called ** trial split ** method. When implemented in Python, it looks like this:

def deterministicPT(n):
    if n <= 1: return False
    if n == 2: return True
    for i in range(2, n):
        if n % i == 0: return False
    return True

If you use this to determine the prime number, it will be as follows.

for n in range(1, 10):
    if deterministicPT(n):
        print(f'{n}Is a prime number')
    else:
        print(f'{n}Is not a prime number')

#1 is not a prime number
#2 is a prime number
#3 is a prime number
#4 is not a prime number
#5 is a prime number
#6 is not a prime number
#7 is a prime number
#8 is not a prime number
#9 is not a prime number

This method adheres to the definition, so you can be confident that $ n $ is a prime number. Such a method is called ** deterministic primality test **. In other words, the definitive primality test is to perform a test with the following properties.

• If it is a prime number, it will definitely pass **
• If it is not a prime number, it will definitely fail **

1.2. Probabilistic primality test

In contrast to the definitive primality test method that determines whether it is a "prime number" or "not a prime number", the algorithm that determines either "may be a prime number" or "not a prime number" is ** a probable primality test * *Is called. The tests used in the probabilistic primality test have the following properties:

• If it is a prime number, ** must pass **
• If it is not a prime number, it can be either pass or fail **
• ** Built-in natural number $ a $ **

And the natural number that passes such a test is called "** probable prime number of base $ a $ **".

Let's look at a concrete example. Consider the following test for the natural number $ n $ you want to determine.

• Fail if $ n = 1 $
• Pass if $ 2 \ leq n \ leq a $
• If $ n> a $, pass if not divisible by $ a $, fail if divisible

This test meets the three properties listed above and can be used in probabilistic primality testing. Implemented in Python, it would look like this:

class ProbabilisticPT:
    def __init__(self, a=2):
        self._a = a
    
    def change_base(self, a):
        self._a = a
    
    def test(self, n):
        if n == 1: return False
        if n <= self._a: return True
        if n % self._a != 0:
            return True
        else:
            return False

Let's make a primality test when $ a = 2 $.

a = 2
ppt = ProbabilisticPT(a)
for n in range(10):
    if ppt.test(n):
        print(f'{n}Is the bottom{a}Probable prime number of')
    else:
        print(f'{n}Is not a prime number')

#1 is not a prime number
#2 is the base 2 probable prime
#3 is the base 2 probable prime
#4 is not a prime number
#5 is the base 2 probable prime
#6 is not a prime number
#7 is the base 2 probable prime
#8 is not a prime number
#9 is the base 2 probable prime

In the probable primality test, the judgment that "it is not a prime number" can be trusted. This is because one of the properties of the test, the kinematic pair of "If it is a prime number, it will definitely pass" is "If it fails, it will not always be a prime number". Therefore, it was decided that 4, 6 and 8 are not prime numbers. However, if it is judged as a "probable prime number", you need to be careful because you cannot get much information from it.

The merit of probable primality test is probably its calculation speed. In this example, you can judge whether it is a "probable prime number" or "not a prime number" by just one division. However, the result obtained is "possible prime number", and I do not know if it is really a prime number. In fact, the composite number 9 is also judged to be a probable prime number.

1.3. To improve accuracy

The probabilistic primality test uses multiple bases to test to improve the accuracy of the test. If it is judged as "not a prime number" at one base, it is confirmed that the natural number is not a prime number, so improvement of judgment accuracy can be expected. Let's actually test for natural numbers less than 30 when the base is 2, 3, 5.

ppt = ProbabilisticPT()
p = {}
for a in [2, 3, 5]:
    p[a] = set()
    ppt.change_base(a)  #Set bottom to a
    for n in range(30):
        if ppt.test(n):
            p[a].add(n)
for k, v in p.items(): print(f'bottom{k}Probable prime number of: {v}')

#Find the intersection of a set of probable primes for each base
print('Probable primes at all bases:', p[2] & p[3] & p[5])



#Probable prime number of base 2: {2, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29}
#Probable prime number of base 3: {2, 3, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20, 22, 23, 25, 26, 28, 29}
#Probable prime number of base 5: {2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 16, 17, 18, 19, 21, 22, 23, 24, 26, 27, 28, 29}
#Probable primes at all bases: {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

A 9 that was judged to be a probable prime number when the base was 2 was confirmed to be a composite number by a test with a base of 3. In this way, the idea of the probable primality test is to improve the probability that the probable prime numbers at all the bases are prime numbers by combining multiple bases to a sufficient accuracy. In the test used this time, by combining three bases (2, 3, 5), it was possible to determine a natural number less than 30 as a prime number with 100% accuracy.

1.4. Fermat test

The one that uses ** Fermat's Little Theorem ** for testing is called ** Fermat test **. Fermat's Little Theorem is:

When $ p $ is a prime number and $ a $ is an integer that is not a multiple of $ p $ ($ a $ and $ p $ are relatively prime)

a^{p-1} \equiv 1 \pmod{p}

Is established.

In other words, what is the Fermat test?

• Pass if $ a ^ {n-1} \ equiv 1 \ pmod {n} $ is satisfied for $ n $ and the coprime natural number $ a $
• Fail if not met

That is. The Fermat test is much more powerful than the test I used earlier, and of the 11,408,012,595 odd composites up to $ 25 \ times 10 ^ 9 $, only 21,853 can pass the Fermat test with a bottom of 2. .. (From "[Japanese Wikipedia-Probable Prime] wiki_prp")

However, the Fermat test also has its drawbacks. A composite number called the ** Carmichael number ** passes any bottom Fermat test. Therefore, no matter how many bottoms are combined, the possibility of misjudgment remains in a size that cannot be ignored in the Fermat test.

1.5. Square root of 1 in the mod world

In preparation for improving the Fermat test, first consider the square root of 1 in the mod world. The square root of 1 in the mod world is $ x $ that satisfies the following congruence formula for two or more natural numbers $ n $.

x^2 \equiv 1 \pmod{n}

Obviously $ x = 1, n-1 $ satisfies this equation, so we call these ** trivial ** square roots and the others ** non-trivial ** square roots.

Consider the non-trivial square root. So consider the case where $ x $ is neither $ 1 $ nor $ n-1 $. For example, when $ n = 15 $, $ x = 4, 11 $ is a non-trivial square root.

\begin{aligned}
(4)^2 &= 16 \equiv 1 \pmod{15}\\
(11)^2 &= 121 \equiv 1 \pmod{15}
\end{aligned}

But what if $ n $ is a prime number? Actually, at this time ** there is no non-trivial square root **. Let's show this by reductio ad absurdum. Now there is a non-trivial square root modulo the prime number $ p $, which we call $ x $. That is, $ x $ is neither $ 1 $ nor $ p-1 $,

x^2 \equiv 1 \pmod{p}

Meet. At this time,

\begin{aligned}
x^2 &\equiv 1 \pmod{p}\\
x^2 - 1 &\equiv 0 \pmod{p}\\
(x + 1) (x - 1) &\equiv 0 \pmod{p}
\end{aligned}

And since $ p $ is a prime number, at least one of $ (x + 1) $ and $ (x -1) $ is divisible by $ p $. But since $ x $ is neither $ 1 $ nor $ p-1 $

\begin{aligned}
(x + 1) \not \equiv 0 \pmod{p}\\
(x - 1) \not \equiv 0 \pmod{p}
\end{aligned}

It will be. That is, both $ (x + 1) $ and $ (x -1) $ are not divisible by $ p $. Therefore, it was shown that there is no non-trivial square root modulo the prime number $ p $ because of the contradiction.

1.6. Fermat's Little Theorem for odd prime numbers

If the natural number you want to judge is 2, it is obvious that it is a prime number, so let's assume that it was processed in advance. At this time, the prime numbers are always odd, so consider Fermat's little theorem for odd prime numbers $ p $. Now, since $ p-1 $ is an even number, we use the natural numbers $ s $ and the odd number $ d $.

p - 1 = 2^s \cdot d

Can be expressed as. Therefore, Fermat's little theorem

a^{2^s \cdot d} \equiv 1 \pmod{p}

And this is

(a^{2^{s-1} \cdot d})^2 \equiv 1 \pmod{p}

You can also see that. As shown in the previous section, there is no non-trivial square root of $ 1 $ modulo the prime $ p $, so

\begin{aligned}
a^{2^{s-1} \cdot d} &\equiv \;\;\:1 \pmod{p}\\
a^{2^{s-1} \cdot d} &\equiv -1 \pmod{p}
\end{aligned}

It will be either. If it was the first, then $ s -1> 0 $

(a^{2^{s-2} \cdot d})^2 \equiv 1 \pmod{p}

It can be seen that this is also

\begin{aligned}
a^{2^{s-2} \cdot d} &\equiv \;\;\:1 \pmod{p}\\
a^{2^{s-2} \cdot d} &\equiv -1 \pmod{p}
\end{aligned}

It will be either. If you repeat this, someday

\begin{aligned}
a^{d} &\equiv \;\;\:1 \pmod{p}\\
a^{d} &\equiv -1 \pmod{p}
\end{aligned}

It will be. The summary up to this point is as follows.

When $ p $ is an odd prime number, use the natural numbers $ s $ and the odd $ d $

p - 1 = 2^s \cdot d

Can be written, and ** always satisfies one of the following congruence formulas modulo $ p $ **.

\begin{aligned}
a^{2^{s-1} \cdot d} &\equiv -1\\
a^{2^{s-2} \cdot d} &\equiv -1\\
\cdots\\
a^{2 \cdot d} &\equiv -1\\
a^{d} &\equiv -1\\
a^d &\equiv \;\;\:1
\end{aligned}

This congruence formula is about $ \ log_2 {p} $ at most, so let's check all **.

From the above, the test for the natural number $ n $ that we want to judge is as follows.

• Fail if $ n = 1 $

• Pass if $ n = 2 $

• Fail if $ n \ geq 3 $ is even

• If $ n \ geq 3 $ is odd, then $ n -1 = 2 ^ s \ cdot d $, and $ n $ as a law for $ n $ and the coprime natural number $ a $.

  \begin{aligned}
  a^{2^{s-1} \cdot d} &\equiv -1\\
  a^{2^{s-2} \cdot d} &\equiv -1\\
  \cdots\\
  a^{2 \cdot d} &\equiv -1\\
  a^{d} &\equiv -1\\
  a^d &\equiv \;\;\:1
  \end{aligned}
  

Pass if you meet at least one of them, fail if you do not meet any of them

The probable primality test using this test is called the ** Miller-Rabin primality test **.

1.7. How to implement the test

Organize the judgments for odd $ n $ of 3 or more. Since the probabilistic primality test method tests with multiple bases and judges that it is a prime number only when all pass, it is only necessary to detect the case of failure in the test process. Now let's say $ n -1 = 2 ^ s \ cdot d $ and the bottom of the test is $ a $. At this time, $ a ^ {2 ^ rd} $ is calculated for $ r $ such that $ 0 \ leq r <s $. First, when $ r = 0 $

a^d \equiv 1\; or\; -1 \pmod{n}

If so, the test for the bottom $ a $ is complete. If not, for $ r = 1, 2, \ cdots, s-1 $

a^{2^rd} \not \equiv 1\; or\; -1 \pmod{n}

Calculate as long as. Also, the termination condition $ r <s $ for $ r $ corresponds to $ 2 ^ rd <n -1 $. If $ a ^ {2 ^ rd} \ equiv -1 \ pmod {n} $, the test for the base $ a $ is complete. But what if it becomes $ a ^ {2 ^ rd} \ equiv 1 \ pmod {n} $? At this time, it has already been confirmed that $ a ^ {2 ^ {r-1} d} $ is neither $ 1 $ nor $ -1 $, so $ a ^ {2 ^ {r-1} d} $ is The non-trivial square root of $ 1 $. Therefore, $ n $ is not a prime number and will be rejected. Also, if it is $ a ^ {2 ^ rd} \ not \ equiv 1 ; or ; -1 \ pmod {n} $ until the end, it will be rejected. The above is summarized in the figure below.

1.8. How to choose the bottom

In the Miller-Rabin primality test, the number of bases is generally decided by oneself, and a natural number $ a $ such that $ a <n $ is selected ** randomly and used as the base. Since there is a trade-off between calculation speed and judgment accuracy, it is desirable to have as few bases as possible while ensuring the required accuracy. Therefore, bottom combinations that can efficiently improve accuracy have been studied. The ACL uses $ a = \ {2, 7, 61 } $ as the base. According to [Jaeschke (1993)] min_2_7_61, the minimum composite number that passes all this bottom test is $ 4759123141 ; (= 48781 \ cdot 97561> 2 ^ {32}) $. Therefore, in the range of $ n <2 ^ {32} $ (range of unsigned 4-byte integer), it can be judged with the precision of $ 100 % $.

1.9. Implementation

Let's implement it. The part that uses pow_mod in ACL is replaced by the Python built-in function pow, which is an equivalent function.

def is_prime(n):
    #Trivial part
    if (n <= 1): return False
    if (n == 2 or n == 7 or n == 61): return True
    if (n % 2 == 0): return False

    d = n - 1  #n is odd
    while (d % 2 == 0): d //= 2  #Find the odd d

    for a in (2, 7, 61):
        t = d  #Hold d because it is also used at other bottoms
        y = pow(a, t, n)

        # a^d = 1, -If it is 1, this loop will not enter
        # a^t = 1, -Repeat until 1
        while (t != n - 1 and y != 1 and y != n - 1):
            y = y * y % n
            t <<= 1
        
        # a^d = 1, -1 passes(t%2 == 0)
        # a^t = -1 passes(y != n - 1)
        if (y != n - 1 and t % 2 == 0):
            return False
    return True



print(is_prime(17))  # True
print(is_prime(1000000007))  # True
print(is_prime(121))  # False
print(is_prime(561))  # False (561 is one of the Carmichael numbers)

#bottom{2, 7, 61}Minimum number of composites that pass the test
print(is_prime(4759123141))  # True

1.10. Comparison of calculation speed

Let's compare it with the definitive primality test method of time complexity $ O (\ sqrt {n}) $. The code used for the comparison is below.

import random

#mirror-Rabin Primality Test (Probable Prime Primality Test)
def pro_is_prime(n):
    if (n <= 1): return False
    if (n == 2 or n == 7 or n == 61): return True
    if (n % 2 == 0): return False
    d = n - 1
    while (d % 2 == 0): d //= 2
    for a in (2, 7, 61):
        t = d
        y = pow(a, t, n)
        while (t != n - 1 and y != 1 and y != n - 1):
            y = y * y % n
            t <<= 1
        if (y != n - 1 and t % 2 == 0):
            return False
    return True

#Deterministic primality test
def det_is_prime(n):
    if n < 2: return False
    if n == 2: return True
    if n % 2 == 0: return False
    for i in range(3, int(n ** 0.5) + 1):
        if n % i == 0: return False
    return True


def random_1():
    l1, r1 = [3, 1 << 16]
    return random.randrange(l1, r1, 2)

def random_2():
    l2, r2 = [(1 << 16) + 1, 1 << 32]
    return random.randrange(l2, r2, 2)

def random_3():
    l3, r3 = [3, 1 << 32]
    return random.randrange(l3, r3, 2)


def main():
    """
    seed_list = [111, 222, 333, 444, 555, \
                 666, 777, 888, 999, 101010]
    """
    random.seed(111)  #Random number fixed
    loop = 10000  #Number of loops 10^4 or 10^6

    for _ in range(loop):
        n = random_1()
        #n = random_2()
        #n = random_3()

        pro_is_prime(n)  #Probabilistic
        #det_is_prime(n)  #Definitive
    

if __name__ == "__main__":
    main()

Measurement method

I used AtCoder's code tests (Python 3.8.2). Random numbers were generated from 3 types of range with 10 types of seed values, and the execution time per 10,000 primality tests was investigated. In the case of even numbers, both are processed first, so only odd numbers are used.

• random_1
• Random number range: $ [3, 2 ^ {16}] $ odd number
• Number of loops: $ 10 ^ 4, 10 ^ 6 $ (Because $ 10 ^ 4 $ may be the main overhead such as import)
• random_2
• Random number range: $ [2 ^ {16}, 2 ^ {32}] $ odd number
• Number of loops: $ 10 ^ 4 $
• random_3
• Random number range: $ [3, 2 ^ {32}] $ odd number
• Number of loops: $ 10 ^ 4 $

Measurement result

The result is shown in the figure below. The numerical value is the time per 10,000 primality tests, and the unit is ms.

If the number you want to judge is as small as $ 10 ^ 5 $, the deterministic primality test seems to be slightly faster. However, as the number increases, the superiority of the Miller-Rabin primality test becomes more pronounced.

2. primitive_root Find the smallest of the primitive roots modulo the prime number $ m $.

2.1. Order

First, let's look at the important term "order" in explaining primitive roots. The definition is this.

For two or more natural numbers $ m $ and $ m $ and relatively prime integers $ a $

a^{n} \equiv 1 \pmod{m}

The smallest natural number $ n $ is called the order ** of $ a $ modulo ** $ m $.

Let's look at a concrete example. When $ m = 5, a = 3 $

\begin{aligned}
3^1  &= 3 & \not \equiv 1 \pmod{5}\\
3^2  &= 9 & \not \equiv 1 \pmod{5}\\
3^3  &= 27 & \not \equiv 1 \pmod{5}\\
3^4  &= 81 & \equiv 1 \pmod{5}
\end{aligned}

So the order of $ 3 $ modulo $ 5 $ is $ 4 $. Also, when $ m = 12, a = 5 $,

\begin{aligned}
5^1  &= 2 & \not \equiv 1 \pmod{12}\\
5^2  &= 25 & \equiv 1 \pmod{12}
\end{aligned}

So the order of $ 5 $ modulo $ 12 $ is $ 2 $.

The nature of the order

The order $ e $ of the integer $ a $, which is relatively prime to $ m $ modulo $ m $, has the following properties for the positive integer $ n $.

a^n \equiv 1 \pmod{m} \Leftrightarrow e\;Is\;n\;Divisor of

(Proof) \Leftarrow: Since $ e $ is a divisor of $ n $, we use the natural number $ k $ and multiply by $ n = ke $. Therefore, using $ m $ as the law

\begin{aligned}
a^n &=a^{ke}\\
&= (a^e)^k\\
&\equiv 1
\end{aligned}

And established. \Rightarrow: We express $ n $ as $ n = qe + r ; (0 \ leq r <e) $ using non-negative integers $ q and r $. At this time, using $ m $ as the law

\begin{aligned}
a^r &\equiv (a^{e})^qa^r\\
&\equiv a^{qe+r}\\
&\equiv 1
\end{aligned}

Will be. If $ 0 <r <e $, then it contradicts that $ e $ is an order (the order is a ** smallest ** integer that satisfies $ a ^ e \ equiv 1 \ pmod {m} $). So $ r = 0 $, that is, $ e $ is a divisor of $ n $.

(End of proof)

Especially when $ m $ is a prime number, the order is a divisor of $ m -1 $ according to Fermat's Little Theorem.

2.2. Primitive root

Consider the case where $ m $ is a prime number. From Fermat's Little Theorem, the order of $ m $ and the relatively prime integer $ a $, modulo $ m $, is always less than or equal to $ m -1 $. Therefore, we specialize in $ a $ whose order is exactly $ m -1 $, and call it a primitive root ** modulo $ m $ **.

For example, when $ m = 7 $, $ a = 3 $

\begin{aligned}
3^1  &= 3  &\not \equiv 1 \pmod{7}\\
3^2  &= 9  &\not \equiv 1 \pmod{7}\\
3^3  &= 27  &\not \equiv 1 \pmod{7}\\
3^4  &= 81  &\not \equiv 1 \pmod{7}\\
3^5  &= 243  &\not \equiv 1 \pmod{7}\\
3^6  &= 729  &\equiv 1 \pmod{7}\\
\end{aligned}

So the order of $ 3 $ modulo $ 7 $ is $ 6 $. So $ 3 $ is a primitive root modulo $ 7 $. There is not always one primitive root. From the first example shown in Orders, we can see that $ 3 $ is a primitive root modulo $ 5 $. On the other hand, when $ a = 2 $

\begin{aligned}
2^1  &= 2 & \not \equiv 1 \pmod{5}\\
2^2  &= 4 & \not \equiv 1 \pmod{5}\\
2^3  &= 8 & \not \equiv 1 \pmod{5}\\
2^4  &= 16 & \equiv 1 \pmod{5}
\end{aligned}

So $ 2 $ is also a primitive root modulo $ 5 $.

Note that the primitive root always exists when $ m $ is a prime number. Check the proof with "Primitive Root Theorem".

Primitive root nature

Primitive root $ g $ modulo the prime number $ m $ has the following properties:

A set of powers of $ g $ modulo $ m $

\{g\pmod{m}, g^2\pmod{m}, \cdots , g^{m - 1}\pmod{m}\}

And an integer divided by $ m $, the set obtained by subtracting $ 0 $ from the set

\{1, 2, \cdots , m - 1\}

Match.

This can be paraphrased as follows: About the natural number $ k $ of $ 1 \ leq k \ leq m-1 $

• g^k \not \equiv 0
• $ g ^ k $ are all different.

The first is clear from the fact that $ g $ and $ m $ are relatively prime. Therefore, I will prove the second one.

(Proof) Now suppose there is a natural number $ k, l $ such that $ k <l \ leq m --1 $, and it is $ g ^ k \ equiv g ^ l \ pmod {m} $. At this time

\begin{aligned}
g^l - g^k &\equiv 0 \pmod{m}\\
g^k(g^{l-k} - 1) &\equiv 0 \pmod{m}\\
g^{l-k} - 1 &\equiv 0 \pmod{m}\\
g^{l-k} &\equiv 1 \pmod{m}
\end{aligned}

It will be. Since $ l --k <m --1 $, the order of the primitive root $ g $ is inconsistent with $ m-1 $. So for $ 1 \ leq k \ leq m-1 $, $ g ^ k $ is all different. (End of proof)

2.3. How to check if it is a primitive root 1

To determine if a natural number $ g $ greater than or equal to 2 is a primitive root modulo $ m $, $ g, g ^ 2, \ cdots, modulo $ m $ from the definition of the primitive root. Just make sure that all g ^ {m --2} $ are not congruent with $ 1 $. Therefore, the implementation of the algorithm for finding the smallest primitive root is as follows.

def naive_primitive_root(m):
    g = 2
    while True:
        for i in range(1, m - 1):
            if pow(g, i, m) == 1:
                g += 1
                break
        else:
            return g


print(naive_primitive_root(7))  # 3
print(naive_primitive_root(23))  # 5
#print(naive_primitive_root(1000000007))  #Very time consuming

This method looks up all powers of $ g $ up to $ m − 2 $, so it takes a lot of time as $ m $ grows.

2.4. How to check if it is a primitive root 2

In fact, it is easier to confirm whether it is a primitive root by using the properties of the primitive root. The properties used are:

Using the prime number $ p_i $ and the natural number $ e_i $, $ m-1 $ is for a prime number $ m $ over $ 3 $.

m-1 = \prod_{i = 1}^n{p_i^{e_i}}

When factored into prime factors, the following holds for natural numbers $ g $ over $ 2 $.

g\;But\;m\;Primitive root modulo
\Leftrightarrow 1 \leq i \leq n\;Against\; g^{\frac{m-1}{p_i}}\not\equiv 1 \pmod{m}

(Proof) \Rightarrow: Since $ p_i $ is a prime factor of $ m-1 $, $ \ frac {m-1} {p_i} $ is an integer and satisfies $ 1 \ leq \ frac {m-1} {p_i} <p --1 $. When $ g $ is a primitive root, $ g ^ {p-1} $ is the power that should be congruent with $ 1 $ for the first time, so against $ 1 \ leq i \ leq n ; ; g ^ {\ frac {m-1} It becomes {p_i}} \ not \ equiv 1 \ pmod {m} $.

\Leftarrow: The kinematic pair of this claim

g is not a primitive root\Rightarrow g^{\frac{m-1}{p_i}}\equiv 1 \pmod{m}There exists i

So I will show this. Now that $ g $ is not a primitive root, the order $ e $ is less than $ m-1 $. As shown in the order property, $ e $ is a divisor of $ m-1 $, so we use $ p_i $, which is the same as the prime factor of $ m-1 $.

e = \prod_{j = 1}^n{p_j^{e_j^{\prime}}}\;\;\;(e_j^{\prime} \leq e_j)

There is at least one $ j $ that can be factored into prime factors, especially $ e_j ^ {\ prime} <e_j $. Let one of these $ j $ be $ i $. At this time

\begin{aligned}
\frac{m-1}{p_i} &= \frac{1}{p_i}\prod_{j}{p_j^{e_j}}\\
&= \frac{1}{p_i}\left(\prod_{j}{p_j^{e_j - e_j^{\prime}}}\right)\left(\prod_j{p_j^{e_j^{\prime}}}\right)\\
&= p_i^{e_i - e_i^{\prime} - 1}\left(\prod_{j \ne i}{p_j^{e_j - e_j^{\prime}}}\right)\left(\prod_j{p_j^{e_j^{\prime}}}\right)
\end{aligned}

Here, $ e_i --e_i ^ {\ prime} --1 \ geq 0 $ and $ e_j --e_j ^ {\ prime} \ geq 0 $

\frac{m-1}{p_i} = (Natural number) \times e

Will be. Therefore, if $ g $ is not a primitive root

g^{\frac{m-1}{p_i}} \equiv \left(g^e\right)^{Natural number}　\equiv 1 \pmod{m}

There exists $ i $.

(End of proof)

From the above

1. Find the prime factor of $ m-1 $.
2. For each prime factor $ p_i $, check if it becomes $ g ^ {\ frac {m-1} {p_i}} \ not \ equiv 1 \ pmod {m} $.

You can determine if $ g $ is a primitive root by following the procedure.

2.5. Implementation

Let's implement it. Like is_prime, pow_mod is replaced by the built-in function pow.

# @param m must be prime
def primitive_root(m):
    if m == 2: return 1
    if m == 167772161: return 3
    if m == 469762049: return 3
    if m == 754974721: return 11
    if m == 998244353: return 3

    # m-Prime factor extraction of 1
    divs = [2]
    x = (m - 1) // 2
    while x % 2 == 0: x //= 2
    i = 3
    while i * i <= x:
        if x % i == 0:
            divs.append(i)
            while x % i == 0: x //= i
        i += 2
    if x > 1: divs.append(x)

    #Find the smallest g that is not congruent with 1 in all prime factors
    g = 2
    while True:
        if all(pow(g, (m - 1) // div, m) != 1 for div in divs): return g
        g += 1


print(primitive_root(7))  # 3
print(primitive_root(23))  # 5
print(primitive_root(1000000007))  # 5

The first few $ m $ to determine seems to be the mods used in convolution.

3. Conclusion

This time we have looked at primality tests and primitive roots. I especially found the idea of the probabilistic primality test method interesting.

Among the internal_math, the ones that I did not touch on this time are written in [internal_math edition ①] qiita_internal_math_1, so please refer to that as well.

Please let us know if you have any mistakes, bugs, or advice.