Number theory

gcd, co-primes

gcd is short for greatest common divisor
If a,b are co-primes, we denote as $(a,b)=1, \text{which means } gcd(a,b)=1$
We can use Euclid algorithm to calculate gcd of two numbers.

def gcd(a,b):
    while b!=0:
        a,b=b,a%b
    return a

Bezout's identity

Let a and b be integers with greatest common divisor d. Then, there exist integers x and y such that ax + by = d. More generally, the integers of the form ax + by are exactly the multiples of d.

we can use extended euclid algorithm to calculate x,y,gcd(a,b)

def xgcd(a,b):
    '''return gcd(a,b),  x,y  where  ax+by=gcd(a,b)'''
    if b==0:return a,1,0
    g,x,y = xgcd(b,a%b)
    return g,y,x-a//b*y

primality_test

Prime Sieve

class primeSieve:
    '''sieve of Eratosthenes, It will be more efficient when judging many times'''
    primes = [2,3,5,7,11,13]
    def isPrime(self,x):
        if x<=primes[-1]:
            return twoDivideFind(x,self.primes)
        while x>self.primes[-1]:
            left = self.primes[-1]
            right = (left+1)**2
            lst = []
            for i in range(left,right):
                for j in self.primes:
                    if i%j==0:break
                else:lst.append(i)
            self.primes+=lst
        return twoDivideFind(x,lst)

def twoDivideFind(x,li):
    a,b = 0, len(li)
    while a<=b:
        mid = (a+b)//2
        if li[mid]<x:a=mid+1
        elif li[mid]>x: b= mid-1
        else:return mid
    return -1

Miller-Rabin

Excerpted from wikipedia:Miller_Rabin_primality_test

Just like the Fermat and Solovay–Strassen tests, the Miller–Rabin test relies on an equality or set of equalities that hold true for prime values, then checks whether or not they hold for a number that we want to test for primality.

First, a lemma about square roots of unity in the finite field Z/pZ, where p is prime and p > 2. Certainly 1 and −1 always yield 1 when squared modulo p; call these trivial square roots of 1. There are no nontrivial square roots of 1 modulo p (a special case of the result that, in a field, a polynomial has no more zeroes than its degree). To show this, suppose that x is a square root of 1 modulo p. Then:
$x^2\equiv1\ (mod\ p)$
$(x-1)(x+1) \equiv 0\ (mod\ p)$

In other words, prime p divides the product (x − 1)(x + 1). By Euclid's lemma it divides one of the factors x − 1 or x + 1, implying that x is congruent to either 1 or −1 modulo p.

Now, let n be prime, and odd, with n > 2. It follows that n − 1 is even and we can write it as 2^s·d, where s and d are positive integers and d is odd. For each a in (Z/nZ)*, either

$a^d\equiv 1\ (mod\ n)$
or
$a^{2^r*d}\equiv -1\ (mod\ n), \text{where } 0\le r<s$

To show that one of these must be true, recall Fermat's little theorem, that for a prime number n:
$a^{n-1}\equiv1\ (mod\ n)$

By the lemma above, if we keep taking square roots of aⁿ⁻¹, we will get either 1 or −1. If we get −1 then the second equality holds and it is done. If we never get −1, then when we have taken out every power of 2, we are left with the first equality.

The Miller–Rabin primality test is based on the contrapositive of the above claim. That is, if we can find an a such that
$a^d\not\equiv 1\ (mod\ n)$
and

$a^{2^r*d}\not\equiv -1\ (mod\ n), \text{where } 0\le r<s$

then n is not prime. We call a a witness for the compositeness of n (sometimes misleadingly called a strong witness, although it is a certain proof of this fact). Otherwise a is called a strong liar, and n is a strong probable prime to base a. The term "strong liar" refers to the case where n is composite but nevertheless the equations hold as they would for a prime.

Every odd composite n has many witnesses a, however, no simple way of generating such an a is known. The solution is to make the test probabilistic: we choose a non-zero a in Z/nZ randomly, and check whether or not it is a witness for the compositeness of n. If n is composite, most of the choices for a will be witnesses, and the test will detect n as composite with high probability. There is, nevertheless, a small chance that we are unlucky and hit an a which is a strong liar for n. We may reduce the probability of such error by repeating the test for several independently chosen a.

For testing large numbers, it is common to choose random bases a, as, a priori, we don't know the distribution of witnesses and liars among the numbers 1, 2, ..., n − 1. In particular, Arnault ^[4] gave a 397-digit composite number for which all bases aless than 307 are strong liars. As expected this number was reported to be prime by the Maple isprime() function, which implemented the Miller–Rabin test by checking the specific bases 2,3,5,7, and 11. However, selection of a few specific small bases can guarantee identification of composites for n less than some maximum determined by said bases. This maximum is generally quite large compared to the bases. As random bases lack such determinism for small n, specific bases are better in some circumstances.

python implementation

from random import sample
def quickMulMod(a,b,m):
    '''a*b%m,  quick'''
    ret = 0
    while b:
        if b&1:
            ret = (a+ret)%m
        b//=2
        a = (a+a)%m
    return ret

def quickPowMod(a,b,m):
    '''a^b %m, quick,  O(logn)'''
    ret =1
    while b:
        if b&1:
            ret =quickMulMod(ret,a,m)
        b//=2
        a = quickMulMod(a,a,m)
    return ret


def isPrime(n,t=5):
    '''miller rabin primality test,  a probability result
        t is the number of iteration(witness)
    '''
    if n<2:
        print('[Error]: {} can\'t be classed with prime or composite'.format(n))
        return
    if n==2: return True
    d = n-1
    r = 0
    while d%2==0:
        r+=1
        d//=2
    t = min(n-3,t)
    for a in sample(range(2,n-1),t):
        x= quickPowMod(a,d,n)
        if x==1 or x==n-1: continue  #success,
        for j in range(r-1):
            x= quickMulMod(x,x,n)
            if x==n-1:break
        else:
            return False
    return True

Factorization

Pollard's rho algorithm

Excerpted from wikipedia:Pollard's rho algorithm

Suppose we need to factorize a number
$n=pq$ , where $p$ is a non-trivial factor. A polynomial modulo $n$ , called
$g(x)=x^2+c\ mod\ n$
where c is a chosen number ,eg 1.

is used to generate a pseudo-random sequence: A starting value, say 2, is chosen, and the sequence continues as
$x_1 = g(2),x_2=g(g(2)),\ldots, x_i =g^{(i)}(2) = g(x_{i-1})$
,
The sequence is related to another sequence $\{x_k\ mod \ p\}$ . Since $p$ is not known beforehand, this sequence cannot be explicitly computed in the algorithm. Yet, in it lies the core idea of the algorithm.

Because the number of possible values for these sequences are finite, both the $\{x_n\}$ sequence, which is mod $n$ , and $\{x_n\ mod\ p\}$ sequence will eventually repeat, even though we do not know the latter. Assume that the sequences behave like random numbers. Due to the birthday paradox, the number of $x_k$ before a repetition occurs is expected to be $O(\sqrt{N})$ , where $N$ is the number of possible values. So the sequence $\{x_n\ mod\ p\}$ will likely repeat much earlier than the sequence $x_k$ . Once a sequence has a repeated value, the sequence will cycle, because each value depends only on the one before it. This structure of eventual cycling gives rise to the name "Rho algorithm", owing to similarity to the shape of the Greek character ρ when the values
$x_1\ mod\ p, x_2\ mod\ p,\ldots$ are represented as nodes in a directed graph.

This is detected by the Floyd's cycle-finding algorithm: two nodes

i,j

are kept. In each step, one moves to the next node in the sequence and the other moves to the one after the next node. After that, it is checked whether

\text{gcd}(x_i-x_j,n)\neq 1

.
If it is not 1, then this implies that there ris a repetition in the

\{x_k\ mod\ p\}

swquence

This works because if the $x_i\ mod\ p$ is the same as $x_j\ mod\ p$ , the difference between $x_i,x_j$ is necessarily a multiple of $p$ . Although this always happens eventually, the resulting GCD is a divisor of $n$ other than 1. This may be $n$ itself, since the two sequences might repeat at the same time. In this (uncommon) case the algorithm fails, and can be repeated with a different parameter.

python implementation

from random import randint
from isPrime import isPrime
from gcd import gcd

def factor(n):
    '''pollard's rho algorithm'''
    if n==1: return []
    if isPrime(n):return [n]
    fact=1
    cycle_size=2
    x = x_fixed = 2
    c = randint(1,n)
    while fact==1:
        for i in range(cycle_size):
            if fact>1:break
            x=(x*x+c)%n
            if x==x_fixed:
                c = randint(1,n)
                continue
            fact = gcd(x-x_fixed,n)
        cycle_size *=2
        x_fixed = x
    return factor(fact)+factor(n//fact)

Euler function

Euler function, denoted as $\phi(n)$ , mapping n as the number of number which is smaller than n and is the co-prime of n.

e.g.: $\phi(3)=2$ since 1,2 are coprimes of 3 and smaller than 3, $\phi(4)=2$ ,(1,3)

Euler function is a kind of productive function and has two properties as follows:

$\phi(p^k) = p^k-p^{k-1}$ , where p is a prime
$\phi(mn) = \phi(m)*\phi(n)$ where $(m,n)=1$

Thus, for every narural number n, we can evaluate $\phi(n)$ using the following method.

factorize n:
$n = \prod _{i=1}^{l} p_i^{k_i}$ , where $p_i$ is a prime and $k_i,l > 0$ .
calculate $\phi(n)$ using the two properties.

$\begin{aligned} \phi(n) &=\phi( \prod _{i=1}^{l} p_i^{k_i}) \\ &=\prod _{i=1}^{l} \phi( p_i^{k_i}) \\ &=\prod _{i=1}^{l} ( p_i^{k_i}-p_i^{{k_i}-1})\\ &=\prod _{i=1}^{l}p_i^{k_i} \prod _{i=1}^{l} ( 1-\frac{1}{p_i})\\ &=n \prod _{i=1}^{l} ( 1-\frac{1}{p_i})\\ \end{aligned}$

And , $\sigma(n)$ represents the sum of all factors of n.
e.g. : $\sigma(9) = 1+3+9 = 14$
$\begin{aligned} \sigma(n) &= \prod _{i=1}^{l} \sum_{j=0}^{k_i} p_i^j \\ &=\prod _{i=1}^{l} \frac{p_i^{k_i+1}-1}{p_i-1}\\ \end{aligned}$

A perfect number is defined as $\sigma(n) = 2n$
The following is the implementation of this two functions.

from factor import factor
from collections import Counter
from functools import reduce
from operator import mul
def phi(n):
    st  = set(factor(n))
    return round(reduce(mul,(1-1/p for p in st),n))

def sigma(n):
    ct = Counter(factor(n))
    return reduce(mul,(round((p**(ct[p]+1)-1)/(p-1)) for p in ct),1)

Modulo equation

The following codes can solve a linear, group modulo equation. More details and explanations will be supplied if I am not too busy.

Note that I use -- to represent $\equiv$ in the python codes.

import re

from gcd import xgcd
from euler import phi
from isPrime import isPrime
from factor import factor

def  ind(m,g):
    ''' mod m ,primary root g  ->  {n:indg n}'''
    return {j:i for i in range(m-1) \
            for j in range(m) if (g**i-j)%m==0}

def gs(m,num=100):
    '''return list of  m's  primary roots below num'''
    p = phi(m)
    mp = factor(p)
    checkLst = [p//i for i in mp]
    return [i for i in range(2,num) if all((i**n-1)%m !=0  for n in checkLst)]

def minG(m):
    p = phi(m)
    mp = factor(p)
    checkLst = [p//i for i in mp]
    i=2
    while  1:
        if all((i**n-1)%m !=0  for n in checkLst):return i
        i+=1

class solve:
    def __init__(self,equ=None):
        self.linearPat= re.compile(r'\s*(\d+)\s*--\s*(\d+)[\s\(]*mod\s*(\d+)')
        self.sol  = []
        #self.m = m
        #self.ind_mp = ind(m,minG(m))
    def noSol(self):
        print('equation {equ} has no solution'.format(equ=self.equ))
    def error(self):
        print("Error! The divisor m must be postive integer")
    def solveLinear(self,a,b,m):
        '''ax--b(mod m): solve linear equation with one unknown
            return  ([x1,x2,...],m)
        '''
        a,b,m = self.check(a,b,m)
        g,x,y=xgcd(a,m)
        if a*b%g!=0:
            self.noSol()
            return None
        sol=x*b//g
        m0 = m//g
        sols = [(sol+i*m0)%m for i in range(g)]
        print('{}x--{}(mod {}), solution: {} mod {}'.format(a,b,m,sols,m))
        return (sols,m)
    def check(self,a,b,m):
        if m<=0:
            self.error()
            return None
        if a<0:a,b=-a,-b  ## important
        if b<0:b+= -b//m * m
        return a,b,m

    def solveHigh(self,a,n,b,m):
        ''' ax^n -- b (mod m)  ind_mp is a dict of  m's {n: indg n}'''
        ind_mp = ind(m,minG(m))
        tmp = ind_mp[b] - ind_mp[a]
        if tmp < 0:tmp+=m
        sol = self.solveLinear(n,tmp,phi(m))
        re_mp = {j:i for i ,j in ind_mp.items()}
        sols = [re_mp[i] for i in sol[0]]
        print('{}x^{}--{}(mod {}),  solution: {} mod {}'.format(a,n,b,m,sols,m))
        return sols,m

    def solveGroup(self,tups):
        '''tups is a list of tongyu equation groups, like
            [(a1,b1,m1),(a2,b2,m2)...]
            and, m1,m2... are all primes
        '''
        mp = {}
        print('solving group of equations: ')
        for a,b,m in tups:
            print('{}x--{}(mod {})'.format(a,b,m))
            if m in mp:
                if mp[m][0]*b!=mp[m][1]*a:
                    self.noSol()
                    return
            else:mp[m] = (a,b)
        product = 1
        for i in mp.keys():
            product *=i
        sol = [0]
        for i in mp:
            xs,m = self.solveLinear(product//i*mp[i][0],1,i)
            new = []
            for x in xs:
                cur = x*product//i*mp[i][1]
                for old in sol:
                    new.append(old+cur)
            sol = new
        sol= [i%product for i in sol]
        print('final solution: {} mod {}'.format(sol,product))
        return sol,product
    def __call__(self):
        s=input('输入同余方程，用--代表同于号，形如3--5(mod 7)代表3x模7同余于5')
        li= self.linearPat.findall(s)
        li = [(int(a),int(b),int(m)) for a,b,m in li]
        print(self.solveLinear(li[0]))


if __name__ == '__main__':
    solver  = solve()
    res = solver.solveLinear(3,6,9)
    print()
    res = solver.solveHigh(1,8,3,11)
    print()
    res = solver.solveGroup([(5,11,2),(3,8,5),(4,1,7)])