Modular multiplicative inverse – Tính nghịch đảo theo modulo

Nghịch đảo của a theo modulo m (Modular Multiplicative Inverse – MMI) là một số nguyên x sao cho:
a^{-1}\equiv x \enspace (\textrm{mod} \, m)

Tương đương với:
ax \equiv aa^{-1}\equiv 1 \enspace (\textrm{mod} \, m)

Nghịch đảo của a theo module m (là x) tồn tại chỉ khi a và m nguyên tố cùng nhau, tức là:
gcd(a,m)=1

Hãy xem xét một vài cách để tính MMI

1. Brute Force
Chúng ta có thể duyệt tìm số x sao cho:
ax \equiv 1 \enspace (\textrm{mod} \, m)

int divMod(int a, int m) {
    a %= m;
    for(int x = 1; x < m; x++)
        if ((a*x) % m == 1) return x;
}

Độ phức tạp: O(M)

2. Sử dụng giải thuật Euclid mở rộng (Extended Euclidean Algorithm – EEA)
Chúng ta cần tìm số x sao cho a·x = 1 (mod m). Biểu thức này có thể viết lại thành a·x = 1 + m·y, tương đương với a·x – m·y = 1. Bởi vì x và y không nhất thiết phải là số dương, vì vậy chúng ta viết lại biểu thức ở dạng chuẩn:
a·x + m·y = 1

Trong số học, bổ đề Bézout phát biểu rằng: Nếu d là ước số chung lớn nhất của 2 số nguyên a, b; thì tồn tại 2 số nguyên x và y sao cho:
a·x + b·y = d
Và để tìm x, y; chúng ta có thể sử dụng EEA.

EEA là một sự mở rộng từ thuật toán Euclid. Bên cạnh việc tìm ước chung lớn nhất của 2 số a và b, EEA còn tìm ra 2 số x và y (thông thường 1 trong 2 số này sẽ âm) thoả mãn bổ đề Bézout: a·x + b·y = gcd(a, b).

EEA đặc biệt hữu dụng khi a và b là 2 số nguyên tố cùng nhau, khi đó x là nghịch đảo của a theo modulo b, và y là nghịch đảo của b theo modulo a (quy về bài toán hiện tại của chúng ta)

EEA kết hợp quá trình tìm GCD(a, b) trong thuật toán Euclid; đồng thời tìm thêm cặp số x, y thoả mãn bổ đề Lemma sử dụng công thức truy hồi:
1. Nếu a = 0, thuật toán kết thúc, trả về kết quả x = 0, y = 1 (Vì lúc này gcd(a, b) = b, vậy ta có ax + by = 0·0 + b·1 = b)
2. Ngược lại:
+ Tính thương q = b / a và số dư r = b % a
+ Đệ quy tìm cặp số x1, y1 sao cho a·x1 + r·y1 chia hết bởi a và r
+ Cuối cùng thuật toán trả về kết quả x = y1 – q · x1 và y = x1 (*)
(*)
Tại hàm đệ quy (a, b), ta cần tìm x, y sao cho a·x + b·y = gcd(a, b) = d.
Hàm này lại gọi đến hàm (b%a, a)
Tại hàm đệ quy (b%a, a), ta cần tìm x1, y1 sao cho:
(b%a)·x1 + a·y1 = gcd(b%a, a) = d.
Vậy ax + by = (b%a)x1 + ay1
Nếu chọn y = x1 (vì sao ? Chọn y khác thì sao ?) , tính toán ta có x = y1 – (b/a)·x1

Code C++:

int64_t extGcd(int64_t a, int64_t b, int64_t& x, int64_t& y) {
    if (!a) {
        x = 0;
        y = 1;
        return b;
    }
    int64_t x1, y1;
    int64_t d = extGcd(b % a, a, x1, y1);
    x = y1 - (b / a) * x1;
    y = x1;
    return d;
}

Độ phức tạp: O(log(max(a, b))

Như vậy để tính nghịch đảo của a theo modulo m, ta sử dụng hàm extGcd như sau:

int64_t x, y;
int64_t g = extGcd(a, m, x, y);
if (g != 1) // bad gcd
else // x là nghịch đảo của a theo modulo m

Tham khảo:
[1] http://vuontoanblog.blogspot.com/2012/11/Euclidean-algorithm1.html
[2] https://comeoncodeon.wordpress.com/2011/10/09/modular-multiplicative-inverse/
[3] http://vnoi.info/wiki/translate/he/Number-Theory

Happy coding,

[NTUCoder] UOCLE – Ước lẻ

Đề bài: http://laptrinh.ntu.edu.vn/Problem/Details/3267

 Lời giải:
+ Đầu tiên ta chứng minh: Số tự nhiên N có số ước là lẻ <=> N là số chính phương.
Thực vậy, phân tích N thành thừa số nguyên tố: N = p1^a1 * p2^a2 … * pm^am, như vậy số ước của N = (a1 + 1) * (a2 + 1) * .. * (am + 1). Muốn tích này lẻ thì mỗi thừa số (ai + 1) phải lẻ => ai chẳn => N có thể được viết dưới dạng N = [p1^(a1/2) * .. * pm^(am/2)] ^ 2 => N chính phương
+ Vậy Số số có ước lẻ nhỏ hơn bằng N = Số số chính phương nhỏ hơn bằng N = sqrt(N) (Các số chính phương này có được bằng cách bình phương các số từ 1..sqrt(N))

Pillai’s arithmetical function – Thuật toán tính hàm Pillai

Trong lý thuyết số, hàm số Pillai được định nghĩa theo công thức sau:

P(n)=\sum_{k=1}^n\gcd(k,n)

Ta có thể viết lại như sau:
P(n) = \sum_{d\mid n} d cnt(d), với d là ước số của n và cnt(d) là số gặp (i,n) thoả mãn gcd(i,n) = d.

Theo tính chất của hàm gcd, nếu ta có gcd(a,b) = d thì gcd(a/d, b/d) = 1. Do đó công thức của hàm Pillai sẽ là:
P(n) = \sum_{d\mid n} d \phi(n/d), với \phiHàm phi Euler

Có thể tính được giá trị hàm Pillai của các số trong khoảng [1..N] bằng cách sàng Sieve tính trước Phi hàm Euler, sau đó sàng thêm một lần nữa để tính P(n).
Code
Độ phức tạp O(NlogN)

Nguồn tham khảo:

[1] https://en.wikipedia.org/wiki/Pillai’s_arithmetical_function

Happy coding

 

Euler’s totient function – Thuật toán tính Phi hàm Euler

Hàm số Euler của một số nguyên dương n (ký hiệu \phi(n)) là số các số i nguyên dương nhỏ hơn hoặc bằng n nguyên tố cùng nhau với n (GCD(i, n) = 1).

Một số công thức tính \phi(n):

  • \phi(n)=n\prod_{p|n}\left(1-\frac{1}{p}\right), với p là ước nguyên tố của n
  • \phi(n)=\sum_{d\mid n} d \cdot \mu(n/d), với d là ước của n và \mu là hàm Mobius

Công thức đầu tiên cho ta ý tưởng về cách tính Phi hàm Euler bằng sàng nguyên tố Sieve với độ phức tạp O(NlogN):
Code

Các nguồn tham khảo:

[1] https://vi.wikipedia.org/wiki/Phi_hàm_Euler

Happy coding

Convert from decimal to base-n number system and vice versa

Convert from decimal to base-n number system
We can implement this with a simple, easy to understand program:

void g( int x, int n ) {
    if ( x ) {
        g( x / n, n );
        std::cout << x % n;
    }
}


Convert from decimal to base-n number system
Let’s  consider an example converting a binary number to decimal number:
10101_(2) = 2^4 * 1 + 2^3 * 0 + 2^2 * 1 + 2^1 * 0 + 2^0 * 1
We also can write it like this:
10101_(2) = 1 + 2 * ( 0 + 2 * ( 1 + 2 * ( 0 + 2 * 1 ) ) )

Expanding for base-n number system, if we have a number in base-n number system: c_{0}c_{1}...c_{m-1}( n )
Then:
c_{0}c_{1}...c_{m-1}( n ) = c_{m-1} + n * ( c_{m-2} + n * ( ... n * c_{0} ) ..
Finally we have a program without the need of pow function:

int h( char * src, int n ) {
    int decimal = 0;
    for( char * p = src; p - src < strlen( src ); p++ ) {
        decimal = decimal * n + *p - '0';
    }

    return decimal;
}