<center><h1>Số Chính Phương </h1></center> Trong bài viết này, mình sẽ nói về số chính phương và bài toán kiểm tra một số có phải là số chính phương hay không bằng những cách từ đơn giản tới phức tạp. # 1. Lý thuyết Số chính phương là số tự nhiên có căn bậc hai là một số tự nhiên, hay nói cách khác, số chính phương bằng bình phương (lũy thừa bậc 2) của một số nguyên. :::info Số n được gọi là số chính phương khi và chỉ khi tồn tại $a \in \mathbb{Z}$ sao cho $a^2 = n$ ::: > Ví dụ như 0, 1, 4, 9, 16, ... Một số tính chất của số chính phương: - Chữ số tận cùng của một số chính phương luôn là 0, 1, 4, 5, 6, 9. - Khi phân tích một số chính phương ra thừa số nguyên tố ta được các thừa số là lũy thừa của số nguyên tố với số mũ chẵn. - Thật vậy, giả sử $n$ là một số chính phương, tức là ta có thể viết $n = a^2$, mà mọi số nguyên dương $> 1$ đều có thể phân tích thừa số nguyên tố, do đó ta có thể viết $a = p_1^{k_1}.p_2^{k_2}....p_n^{k_n}$ nên $n = (p_1^{k_1}.p_2^{k_2}....p_n^{k_n})^2 = p_1^{2k_1}.p_2^{2k_2}....p_n^{2k_n}$. Hệ quả là tất cả số mũ đều sẽ là một số chẵn (đây là một tính chất quan trọng mà chúng ta sẽ dùng) - Số chính phương chia cho 3 luôn có số dư là 0 hoặc 1; chia cho 4 luôn dư 0 hoặc 1 - Mọi số chính phương lẻ khi chia cho 8 luôn luôn dư 1. - Một số tự nhiên là số chính phương khi và chỉ khi nó có số lượng ước số là một số lẻ. # 2. Xây dựng thuật toán kiểm tra số chính phương Bài toán kiểm tra một số có phải là số chính phương hay không thì không phải là một bài toán khó, có rất nhiều cách để giải, ta sẽ đi qua từng thuật toán. ## 2.1 Thuật toán đơn giản Ta dùng các hàm có sẵn như `sqrt()` hay `math.isqrt` để tính căn bậc 2 của số nguyên dương $n$, sau đó bình phương nó lên, nếu hai giá trị đó bằng nhau thì đó là số chính phương. ```cpp= bool is_square(long long int n){ if(n < 0) return false; if(n == 0) return true; long long int a = sqrtl(n); if(a * a == n) return true; else return false; } ``` Về mặt lý thuyết, các hàm tính căn bậc 2 có độ phức tạp thuật toán là $O(1)$, nên thuật toán trên là một thuật toán tốt. Nhưng vấn đề đặt ra là nếu đầu vào là một số nguyên rất lớn, lớn đến mức ngoài vùng lưu trữ của kiểu dữ liệu `long long` luôn thì sao?, ta không thể dùng cách đó được nữa. Có một giải pháp đơn giản hơn là chuyển sang python, vì python cho phép ta lưu một số nguyên có giá trị lớn hơn nhiều so với c++, thế nhưng nên ta dùng hàm `sqrt` trong python với đầu vào là số $n$ có độ dài bit tùy ý (BigInt), độ phức tạp của nó sẽ thường là $O(n \log n)$, nên nếu giả sử ta cần phải kiểm tra một đầu vào gồm rất nhiều số lớn, cách trên sẽ không tối ưu, ta cần những phương pháp như tìm kiếm nhị phân hoặc các phương pháp tìm nghiệm xấp xỉ như Newton-Raphson, ... các phương pháp đó tuy hiệu quả, nhưng mình sẽ không bàn đến ở đây, vì nó quá nổi tiếng, mình sẽ giới thiệu 2 phương pháp rất đơn giản, nhưng hiệu quả. ## 2.2 Phương pháp xác suất (Thuật toán Monte Carlo) Phương pháp này dựa trên một tính chất quan trọng: ::: info Nếu $n$ là một số chính phương, thì với mọi số nguyên dương $m$, giá trị $n \pmod m$ bắt buộc phải là một thặng dư bậc hai (Quadratic Residue) modulo $m$. ::: **Ý tưởng của thuật toán:** Thay vì cố gắng tính căn bậc hai, chúng ta sẽ thực hiện kiểm tra tính chất thặng dư bậc hai với các số nguyên tố $p$ ngẫu nhiên. Nếu $n$ thất bại ở bất kỳ lần thử nào, nó chắc chắn không phải số chính phương. Nếu nó vượt qua nhiều lần thử, ta có thể tin tưởng nó là số chính phương với **xác suất sai số phụ thuộc vào số lần thử** - Bước 1: Chọn một số nguyên tố $p$ ngẫu nhiên. - Bước 2: Kiểm tra tính chất thặng dư bậc hai: Nếu $n^{(p-1)/2} \equiv p-1 \pmod p$, thì $n$ không phải là số chính phương (theo tiêu chuẩn Euler). - Bước 3: Lặp lại quá trình trên $k$ lần với các số $p$ khác nhau. - Xác suất để một số không phải số chính phương vượt qua bài test thành công qua $k$ lần thử là khoảng $(1/2)^k$. Với $k=20$, xác suất sai số chỉ còn khoảng $1/1.048.576$. ```py def is_square(n, k=20): if n < 0: return False if n < 2: return True for _ in range(k): p = get_random_prime() if pow(n, (p - 1) // 2, p) == p - 1: return False return True ``` Nhận xét: Thuật toán này đã rất mạnh để loại bỏ nhanh các số không phải số chính phương. Tuy nhiên, chi phí để chọn số nguyên tố ngẫu nhiên và thực hiện phép tính lũy thừa modulo (pow) vẫn còn khá lớn nếu ta cần kiểm tra hàng triệu số mỗi giây. Nếu ta muốn tối ưu thêm nữa thì sao? ## 2.3 Thật toán kiểm tra nhanh (chưa nghĩ ra tên) Từ những tính chất trên, ta có thể nhận ra rằng việc kiểm tra một số là số chính phương sẽ khó hơn việc chứng minh nó không phải chính phương, do đó thay vì hỏi câu hỏi "Trong các số của input, số nào là số chính phương" thì ta sẽ hỏi "Trong tập dữ liệu đó thì số nào không phải là số chính phương", việc thay đổi góc nhìn như thế sẽ giúp ta giải quyết bài toán này dễ hơn. Đây là hàm kiểm tra mà mình tự custom: ```py= import math _MODULI = [735, 899, 793, 583, 851] _MASK_64 = 0x0202021202030213 _M0, _M1, _M2, _M3, _M4 = _MODULI _MSK0 = 0x60008a0004402102410806084830004500022010812084030424180022804110084090422182120c00114000880420482100c10906000ca00044021024108060848301045000220108120840304241800228001100a409042018213 _MSK1 = 0x204038620200e624020129c80021cb00094632241a034000711e881088340820c4224043281224014081490e008b82180416201002620801169c90483ca120446330022830e08151cd400a825a001c439120328e080014144890e220012114122201a89060874089c80a205a121102b3 _MSK2 = 0x1613a09402483221b80d102042461190c480a40200380cc826c12010b0c40660366100059822830013000504e815400b08502843400aa05c8280032003051066800219b019808c3420120d904cc070010094048c262189081022c07611304900a41721b _MSK3 = 0x334408c009a0a432804002190728c50120306408409819014069092180440e914a1b060060158230540011022a4300e9018014768c9158084068a405a2a0660651822947002018a13 _MSK4 = 0x130040a668a040150201a201211b0a008e04c978001a6843042082624000168290a100305d02000604cc3c001830412010a02b21008740992884324e0a442210445c211a98123100a4498108478288ac00311700621000892c008a48513080226ac14403c2100e01121b def is_square(n): if n < 0: return False if n == 0: return True if not (_MASK_64 & (1 << (n & 63))): return False a = n tz = (n & -n).bit_length() - 1 if tz: if tz & 1: return False a >>= tz if (a & 7) != 1: return False if not (_MSK0 & (1 << (n % _M0))): return False if not (_MSK1 & (1 << (n % _M1))): return False if not (_MSK2 & (1 << (n % _M2))): return False if not (_MSK3 & (1 << (n % _M3))): return False if not (_MSK4 & (1 << (n % _M4))): return False root = math.isqrt(n) return root * root == n ``` Mình sẽ giải thích chi tiết từng phần Đầu tiên là ```py= if not (_MASK_64 & (1 << (n & 63))): return False ``` Lệnh này có thể được chia làm 3 bước: Bước 1: n & 63 tương đương với $n \pmod{64}$. Bước 2: 1 << (n & 63) (Tạo con trỏ bit) Sau khi có số dư $r = n \pmod{64}$ (giá trị từ 0 đến 63), lệnh này sẽ dịch bit 1 sang trái $r$ vị trí. Mục đích: Tạo ra một mask mà chỉ có duy nhất bit ở vị trí thứ $r$ được bật lên 1. Nó đóng vai trò như một con trỏ để truy vấn vào bảng lookkup của chúng ta. Bước 3: _MASK_64 & ... (Truy vấn cơ sở dữ liệu)_MASK_64 là một hằng số 64 bit (0x0202021202030213). Đây không phải là một con số ngẫu nhiên, nó chính là một Lookup Table (Bảng tra cứu) thu nhỏ, mỗi bit trong 64 bit của hằng số này đại diện cho một số dư từ 0 đến 63. - Nếu bit thứ $i$ bằng 1: Nghĩa là $i$ là một Thặng dư bậc hai mod 64. - Nếu bit thứ $i$ bằng 0: Nghĩa là không có số chính phương nào khi chia 64 lại dư $i$.2. - Ý nghĩa: Theo lý thuyết số, các số chính phương $k^2 \pmod{64}$ chỉ có thể nhận một trong 12 giá trị sau:$$\{0, 1, 4, 9, 16, 17, 25, 33, 36, 41, 49, 57\}$$ - Hằng số _MASK_64 đó có các bit tại đúng 12 vị trí này được bật lên 1.Nếu $n \pmod{64}$ rơi vào một trong 12 số này thì kết quả phép & sẽ khác 0.Nếu $n \pmod{64}$ rơi vào 52 số còn lại thì kết quả phép & bằng 0 vậy ta return luôn là False ngay lập tức. Tỷ lệ loại bỏ bằng cách này là: $1 - (12/64) = 81.25\%$ chỉ bằng 1 dòng code Phần tip theo: ```py= a = n tz = (n & -n).bit_length() - 1 if tz: if tz & 1: return False a >>= tz ``` Ta có hệ quả cơ bản mọi số nguyên dương $n$ đều có thể phân tích duy nhất dưới dạng:$$n = 2^{tz} \cdot a$$Trong đó $tz$ là số lượng bit 0 tận cùng (Trailing Zeros) và $a$ là một số lẻ. Điều kiện 1: Số bit 0 tận cùng ($tz$) phải là số chẵn. Dựa trên định lý qmà ta đã viết ở Phần 1: "Phân tích ra thừa số nguyên tố thì các số mũ phải chẵn". Nên nếu $n$ là số chính phương, thì thừa số nguyên tố $2$ trong $n$ phải có số mũ chẵn.Trong hệ nhị phân, số mũ của thừa số 2 chính là số lượng bit 0 ở cuối. ```py= tz = (n & -n).bit_length() - 1 ``` Đây là một cách "trick" trong Python để tìm số lượng bit 0 tận cùng mà không cần dùng vòng lặp, tận dụng cách máy tính lưu trữ số âm (Số bù 2). ```py= if tz & 1: return False ``` - Nếu $tz$ lẻ (ví dụ $8 = 2^3$, $32 = 2^5$), chúng ta loại ngay lập tức. Điều kiện 2: Phần lẻ $a$ phải đồng dư $1 \pmod 8$.Đây là một tính chất cực mạnh mà thường ta ít để ý. Sau khi đã loại bỏ hết các thừa số 2 bằng phép dịch bit (a >>= tz), chúng ta thu được phần lẻ $a$. Và ta có tính chất mọi số chính phương lẻ khi chia cho 8 luôn dư 1. **Chứng minh**: Một số lẻ luôn có dạng $2k+1$, bình phương của nó là $(2k+1)^2 = 4k^2 + 4k + 1 = 4k(k+1) + 1$. Vì $k(k+1)$ là tích hai số liên tiếp nên luôn chia hết cho 2, dẫn đến $4k(k+1)$ chia hết cho 8. Do đó, $(2k+1)^2 \equiv 1 \pmod 8$. ```py= if (a & 7) != 1: return False ``` - Phép & 7 tương đương với mod 8. Cuối cùng là ta sẽ sử dụng kiến thức về thặng dư bậc 2. Dù đã qua hai lớp lọc, vẫn có những con số không phải chính phương nhưng lại có cấu trúc bit rất "đẹp" (ví dụ: các số có dạng $k^2 + 128$). Để xử lý những trường hợp này, chúng ta áp dụng lý thuyết thặng dư bậc hai trên các hệ cơ số khác nhau. Tại sao lại chọn bộ số [735, 899, 793, 583, 851]? Đây không phải là những con số ngẫu nhiên. Mình đã thử brute-force so sánh nhiều tổ hợp với nhau và thu được kết quả tốt nhất là bộ này, chúng được tính toán dựa trên hai tiêu chí: - Tính nguyên tố cùng nhau : Các bộ số này không chia hết cho nhau, đảm bảo mỗi hệ mod sẽ không bị trùng lặp với hệ trước đó. - $735 = 3 \times 5 \times 7^2$ - $899 = 29 \times 31$. - $793 = 13 \times 61$. - $583 = 11 \times 53$. - $851 = 23 \times 37$. Việc chọn các số mod là tích của các số nguyên tố giúp chúng ta kiểm tra tính chính phương trên nhiều trường số cùng lúc chỉ với một phép tính chia dư. Ban đầu mình lưu trữ một danh sách các số dư và dùng vòng lặp để kiểm tra, nhưng qua trao đổi với đứa bạn học pwn của mình (@nctgaming) thì tụi mình đã tạo ra đoạn code tà đạo hơn giúp chạy nhanh hơn bằng một kỹ thuật gọi là bitmask: ```py= _MODULI = [735, 899, 793, 583, 851] _MASK_64 = 0x0202021202030213 _M0, _M1, _M2, _M3, _M4 = _MODULI _MSK0 = 0x60008a0004402102410806084830004500022010812084030424180022804110084090422182120c00114000880420482100c10906000ca00044021024108060848301045000220108120840304241800228001100a409042018213 _MSK1 = 0x204038620200e624020129c80021cb00094632241a034000711e881088340820c4224043281224014081490e008b82180416201002620801169c90483ca120446330022830e08151cd400a825a001c439120328e080014144890e220012114122201a89060874089c80a205a121102b3 _MSK2 = 0x1613a09402483221b80d102042461190c480a40200380cc826c12010b0c40660366100059822830013000504e815400b08502843400aa05c8280032003051066800219b019808c3420120d904cc070010094048c262189081022c07611304900a41721b _MSK3 = 0x334408c009a0a432804002190728c50120306408409819014069092180440e914a1b060060158230540011022a4300e9018014768c9158084068a405a2a0660651822947002018a13 _MSK4 = 0x130040a668a040150201a201211b0a008e04c978001a6843042082624000168290a100305d02000604cc3c001830412010a02b21008740992884324e0a442210445c211a98123100a4498108478288ac00311700621000892c008a48513080226ac14403c2100e01121b def is_square(n): if not (_MSK0 & (1 << (n % _M0))): return False if not (_MSK1 & (1 << (n % _M1))): return False if not (_MSK2 & (1 << (n % _M2))): return False if not (_MSK3 & (1 << (n % _M3))): return False if not (_MSK4 & (1 << (n % _M4))): return False ``` - Mỗi Modulo $m$ sẽ có một mask $2^m$ bits.Bit thứ $r$ sẽ được bật lên 1 nếu $r$ là thặng dư bậc hai của $m$, tụi mình đã tính toán và chuyển nó thành hex - `if not (_MSK0 & (1 << (n % _M0))): return False`. Phép toán này kiểm tra xem số dư $n \pmod m$ có nằm trong danh sách của chúng ta hay không. CPU thực hiện việc này chỉ trong nháy mắt. Và cuối cùng là sử dụng hàm căn ở bước cuối: ```py= root = math.isqrt(n) return root * root == n ``` Qua thực nghiệm với 10 triệu số ngẩu nhiên, tốc độ xử lý: 6.94 giây cho 10 triệu số $\approx$ 0.69 µs (micro-giây) cho mỗi phép thử, nhanh hơn dùng hàm isqrt thường rất nhiều, và tỉ lệ các số bị lọc là khoảng $99.98377%$, tức là cứ 6161 số đi vào, chỉ có 1 số duy nhất lọt qua được để đi đến hàm math.isqrt(). Nếu không có bộ lọc: ta phải tính isqrt 10 triệu lần. Nếu có bộ lọc: chỉ tính isqrt 1623 lần, ta sẽ tiết kiệm được rất nhiều thời gian, tài nguyên tính toán. --- Qua đó ta đã giải quyết được bài toán kiểm tra tính chính phương của một số nguyên một cách hiệu quả hơn bằng cách thay đổi góc nhìn, bên cạnh đó là sử dụng các phép thao tác trên bit thay vì sử dụng các phép mod để tối ưu tốc độ xử lý. Bài viết này chỉ là một cái chia sẻ của tụi mình trong lúc rảnh rỗi, có thể sẽ có rất nhiều lỗi nên mọi người tham khảo thôi nhé.