Sàng Eratosthenes số nguyên tố trong mảng một chiều

Chào các anh chị/bạn trong forum, e là sinh viên năm nhất và có một bài tập về nhà là duyệt một mảng một chiều, hiển thị ra các phần tử của mảng là số nguyên tố.
e cũng đã đọc thuật toán trong giáo trình( là thuật toán đếm ước của A[i] từ 2-> căn A[i]), e cảm thấy nó chưa tôi ưu lắm nên đã tham khảo phương pháp sàng Eratosthenes. Nhưng các ví dụ của phương pháp này đều là tìm các số nguyên tố trong đoạn [2,n] liên tục, chứ không có nhiều ví dụ về việc ứng dụng nó trong sàng các phần tử của mảng. Nên e đã biến đổi nó một chút thành thuật toán:
Duyệt từ đầu đến cuối mảng, nếu nó là số nguyên tố thì tìm và xóa các phần tử còn lại là bội của nó, còn nếu tại vị trí đó không phải là số nguyên tố thì tìm các phần từ bằng nó còn lại trong mảng và xóa. Sau cùng thì được một mảng chứa toàn bộ là số nguyên tố.
Mong anh/chị hay các bạn nào có phương pháp tối ưu hơn thì có thể chia sẻ để e phát triển bài toán hay hơn ạ. hì hì
Dưới đây là code C của e.

#include <stdio.h>
#include <math.h>
//Define kieu du lieu bool
#define bool int
#define TRUE 1
#define FALSE 0
// Max cua mang A
#define MAX 200

void delete(int A[], int kichThuocMang, int viTriMuonXoa){
    int i;
    for(i = viTriMuonXoa; i < kichThuocMang -1; i++ ){
        A[i] = A[i+1];
    }
}

bool laSoNguyenTo(int n){
    if(n <= 1){
        return FALSE;
    }
    else if(n == 2){
        return TRUE;
    }
    else if(n%2 == 0){
        return FALSE;
    }
    else
    {
        int m = (int)sqrt(n);
        for(int i = 3; i <= m; i+= 2){
            if(n%i == 0){
                return FALSE;
            }
        }
        return TRUE;
    }
}

void nhapMang(int A[], int n){
    for(int i = 0; i < n; i++){
        printf("Nhap A[%d] = ", i);
        scanf("%d", &A[i]);
    }
}

void xuatMang(int A[], int n){
    for(int i = 0; i < n; i++){
        printf("%4d", A[i]);
    }
}



int main()
{
    int A[MAX];
    int n;
    printf("Nhap kich thuoc mang: ");
    scanf("%d", &n);
    // Nhap mang A
    nhapMang(A, n);

    //Mang sau khi nhap
    printf("\nMang sau khi nhap la: \n");
    xuatMang(A, n);

    // Thuc thi thuat toan
    int dem = 0; // dem = so phan tu bi xoa
    for(int i = 0; i < n; i++){
        if( laSoNguyenTo(A[i]) ){ //Neu la so nguyen to thi tim va xoa cac phan tu khac la boi cua no
            for(int j = i+1; j < n; j++){
                if( A[j]%A[i] == 0 ){
                    delete(A, n, j);
                    dem++;
                    n--; // Moi lan xoa se tru di 1 phan tu
                }
            }
        }
        else // Khong la so nguyen to thi tim va xoa nhung phan tu trong mang = no
        {
            for(int j = i+1; j < n; j++){
                if(A[j] == A[i]){
                    delete(A, n, j);
                    dem++;
                    n--; // Moi lan xoa se tru di mot phan tu
                }
            }
        }
    }

    // Mang sau khi loc
    printf("\nCac so nguyen to co trong mang la : \n");
    xuatMang(A, n);
    printf("\nSo Phan tu trong mang bi xoa la %d : \n", dem);

    return 0;
}

Bạn chỉ đang gượng ép thôi :slight_smile:

Khái niệm sàng thực sự có nghĩa là: từ f(p), ta sẽ thao tác với f(p+k), f(p+2k), …; hay f(2p), f(3p), … (như tổng ước số)

2 Likes

Việc kiểm tra, tìm và delete giúp thời gian giải thuật của bạn nâng lên O(sqrt(n) * n^3).
Bằng khi ứng dụng sàng vào, bạn sẽ có 1 mảng đánh dấu 0 và 1 các số là nguyên tố hay ko. Từ đó việc kiểm tra nguyên tố chỉ tốn O(1). Và việc lặp qua mảng kiểm tra chỉ tốn O(n). Tuy nhiên khi này bộ nhớ có độ phức tạp là O(n).

Về sàng ép sàng thành O(n) thì tính ra cả thảy chương trình chỉ tốn O(n) (Nếu ko được thì tốn O(n*log(log(n))

Và đây là một ví dụ rất dễ hiểu cho kỹ thuật gọi là dynamic programming (tiếng việt là Quy Hoạch Động đó) :3

4 Likes

Em hiểu ý của mod là vì bước tìm và delete của e lại khiến tăng thời gian giải thuật lên phức tạp hơn ạ. Để giải quyết thì phải tạo mảng B, khi check thì nếu là nguyên tố thì lưu B[i] = 1 ngược lại thì 0.
Sau đó duyệt mảng B để đưa ra các phần tử nguyên tố.
Không biết e hiểu như vậy có đúng ko ạ, vì e chưa học sâu về giải thuật nâng cao như quy hoạch động nên e không rõ ạ.

1 Like

Gần đúng rồi đó bạn. Thay vì duyệt và đưa ra 1 mảng C gồm những số nguyên tố rồi lấy mảng A tìm trong mảng C.
Thì bạn dò từng phần tử trong mảng A với B luôn

if (B[A[i]]) -> số nguyên tố và chỉ tốn O(1)

À còn QHD thì khi nào học sẽ nghiệm ra tại mình giới thiệu vậy cho bạn biết thêm thôi :stuck_out_tongue:

1 Like

Pre-computation? :smiley:

Ừ thì mỗi truy vấn là O(1), nhưng mem xơi cũng rất khiếp mới được O(1). Chắc cũng trường hợp đặc biệt nào đó thì hay chứ tổng quát chả ai dùng sàng :slight_smile:

2 Likes

Bạn ấy nói về sàng nên mình tiện luôn.
Chứ A[i] mà cỡ vài triệu thì sàng chiếm biết bao nhiêu bộ nhớ rồi.
Lúc đó tìm cách giảm big-O kiểm tra nguyên tố xuống nhỉ. :thinking:

<(") Mà thôi Off-topic rùi, lên discord tám đê rogp10

Ngoài cách bình thường và cách Eratosthenes để kiếm số nguyên tố, mình còn biết thêm 1 cách nữa. Cách này được gọi là Fermat's Test of Primality
Nếu như bạn đã học qua number theory hay toán rời rạc thì bạn có thể đã nghe qua 1 theorem gọi là Fermat’s Little Theorem. Nó có phương trình như sau:

với điều kiện n là một số nguyên tố và a là số nguyên tố cùng nhau với n, i.e gcd(n,a) = 1.

Nếu như n đúng là 1 số nguyên tố, thì phương trình phía trên luôn luôn đúng cho tất cả các a (nguyên tố cùng nhau với n of course :smile:). Ví dụ:

5 là số nguyên tố, cho a = 2, 3, 4, 6

2^(5-1) % 5 = 16 % 5 = 1 
3^(5-1) % 5 = 81 % 5 = 1 
4^(5-1) % 5 = 256 % 5 = 1
6^(5-1) % 5 = 1296 % 5 = 1

Còn nếu như n không phải là nguyên tố (hợp số) thì phương trình có thể sẽ đúng, có thể sẽ sai :sweat_smile:. Ví dụ:

15 không phải là hợp số, cho a = 4

4^(15-1) % 15 = 268435456 % 15 = 1

Bạn có thể đã nhìn vô a = 4 và nghĩ rằng, ối cái phương trình Fermat đúng rồi, vậy số 15 là nguyên tố sao? :thinking:. Hmmm… vậy bây giờ mình sẽ bỏ vô thử a = 5 và a = 6

5^(15-1) % 15 = 6103515625 % 15 = 10
6^(15-1) % 15 = 78364164096 % 15 = 6

Oh wow, bây giờ phương trình đã sai rồi. Bây giờ bạn mới kết luận 15 là một hợp số.

Vậy rốt cuộc bây giờ là sao? Theo những gì chúng mình đã thấy phía trên, một con số n, nếu như nó nguyên tố thì phương trình Fermat sẽ luôn luôn đúng cho tất cả các số a nào mà nguyên tố cùng nhau nới nó (bằng 1). Nếu như n có con số a nào mà làm đáp án nó sai (không bằng 1), thì lập tức, n là hợp số.

Thuật toán:

Cho một con số n, bạn cần tạo con số a nguyên tố cùng nhau với n, rồi bỏ n và a vô phương trình Fermat.

  • Nếu như phương trình không trả về 1, n là số hỗn hợp, exit.
  • Nếu như phương trình trả về 1, n có thể là số nguyên tố.
    • Tạo ra một số a nguyên tố cùng nhau với n khác, rồi test nó lại lần nữa. Bạn làm càng nhiều a, thì khả năng n là số nguyên tố càng cao hơn
    • Khi nào bạn thấy đủ rồi (bạn đã loop tới 10-15-25 lần rồi mà phương trình cứ đưa về 1 hoài). Exit và gọi n là số nguyên tố

Pseudocode

1. Chọn 1 con số a ngẫu nhiên từ [2,3, ... , n-1]
2. Tính GCD(a,n)
   * Nếu như GCD(a,n) ! = 1, n là số hỗn hợp, exit
3. Else, tính a^(n-1) % n (Fermat)
   * Nếu như Fermat != 1, n là số hỗn hợp, exit
   * Nếu như Fermat == 1, n có thể là số nguyên tố
     - Quay lại step 1
            or
     - Nếu như bạn thấy loop đã nhiều rồi, tuyên bố n là số nguyên tố, exit
4 Likes

Đã dùng Fermat thì sao không dùng Rabin-Miller?

Số Carmichael có thể ăn được hết test Fermat, nhưng sẽ không ăn được hết RM => RM tốt hơn. RM khác hai chỗ:

  • Bước nhân mũ bt bằng +/-1 là xong.
  • Bước bình phương: nếu gặp -1 là số nguyên tố (vì x^2 = 1 mod số nguyên tố chỉ có hai nghiệm +/-1) (xem note, bđ Euclid) còn khác -1 mà bay lên 1, nghĩa là không phải. Vì như vậy là x^2 = 1 (mod p) có nghiệm khác -1

Test RM thực ra vẫn là a^(n-1) ?=? 1 (mod n), nhưng logic chi tiết hơn.

  1. Chia thử số n.
  2. Tách n-1 = 2^m * s với s lẻ.
  3. Nếu a^s = +/-1 (mod n) trả về true.
  4. t := a^s mod n
  5. for i=1 to m-1 do t = t^2 mod n, trả về “toạch” if t=1 else trả về “true” if t=n-1
  6. Nếu t^2 mod n = 1 trả về true ngược lại là false.
note

x^2 = 1 (mod p) <=> p | (x^2 - 1) <=> p | (x-1)(x+1) <=> p | (x-1) || p | (x+1) (BĐ Euclid) <=> x = +/-1 (mod p).

4 Likes

Thú thật mình chưa hiểu rõ về Rabin-Miller để giải thích nó một cách dễ dàng cho người ta được. Fermat là algo dễ nhất, nên mình cũng muốn giới thiệu nó đầu tiên vậy :smile:. Ở đây mình chỉ mong số Carmichael không xuất hiện hay con số bạn OP test đều tương đối nhỏ hết.

Khi bạn ấy nắm vững Fermat rồi, thì bạn ấy có thể lên luôn R-M để test mấy con số lớn hơn.

2 Likes

Miller-Rabin cũng từ Fermat’s little theorem mà suy ra thôi, có điều a test ở đây thuộc đoạn [2, n-2] là số nào trong đoạn đó cũng được,

từ Fermat’s little theorem ta có: ap-1 ≡ 1 (mod p) với p là số nguyên tố

vì p là số nguyên tố, p là số lẻ (p = 2 loại ra dễ dàng), nên p-1 là số chẵn. Vậy ap-1 = a(p-1)/2 * a(p-1)/2.

Ta có: ap-1 ≡ 1 (mod p) nên a(p-1)/2 ≡ 1 (mod p) hoặc a(p-1)/2 ≡ -1 (mod p). Chứng minh dễ dàng:

ap-1 = kp + 1,
ap-1 - 1 = kp
(a(p-1)/2 - 1)(a(p-1)/2 + 1) = kp

hay nói cách khác a(p-1)/2 - 1 chia hết cho p hoặc a(p-1)/2 + 1 chia hết cho p, nên a(p-1)/2 ≡ 1 (mod p) hoặc a(p-1)/2 ≡ -1 (mod p).

tổng quát hóa nếu x2 ≡ 1 (mod p) thì x ≡ 1 (mod p) hoặc x ≡ -1 (mod p)

đi sâu hơn nữa: a(p-1)/2 ≡ -1 (mod p) thì ta ko thể phân tích x thêm được nữa, nhưng a(p-1)/2 ≡ 1 (mod p) có thể phân tích thêm được, nếu (p-1)/2 là số chẵn. Khi đó ta lại có a(p-1)/4 ≡ 1 (mod p) hoặc a(p-1)/4 ≡ -1 (mod p). Tiếp tục nếu (p-1)/4 chẵn, ta lại có a(p-1)/8 ≡ 1 (mod p) hoặc a(p-1)/8 ≡ -1 (mod p) v.v… tới khi nào (p-1)/2s là số lẻ thì dừng. Đặt d = (p-1)/2s, hay p-1 = 2s * d. Khi đó ta có:

  • ad ≡ 1 (mod p), hoặc
  • ad ≡ -1 (mod p), hoặc
  • a2d ≡ -1 (mod p), hoặc
  • a4d ≡ -1 (mod p), hoặc
  • a8d ≡ -1 (mod p), hoặc
  • a(p-1)/2 ≡ -1 (mod p)
    Trong đó (p-1)/2 = 2s-1 * d

đây là căn bản của Miller-Rabin test. Nếu n là số nguyên tố thì với mọi 1 < a < n-1 hay 2 ≤ a ≤ n-2 phải thỏa 1 trong các trường hợp ở trên. Nếu có số a nào ko thỏa tất cả các trường hợp đó thì n ko phải là số nguyên tố.

ví dụ n = 217:
n-1 = 216 = 23 * 27. d = 27, s = 3. Với số a thuộc [2, 215], ta cần kiểm tra

  • a27 ≡ 1 (mod 217), hoặc
  • a2^0 * 27 = a27 ≡ -1 (mod 217), hoặc
  • a2^1 * 27 = a54 ≡ -1 (mod 217), hoặc
  • a2^2 * 27 = a108 ≡ -1 (mod 217)
    với a = 211, ta có 21127 ≡ 1 (mod 217) nên 211 “chém gió” 217 là số nguyên tố
    với a = 2, kết quả lần lượt là 190, 190, 78, 8 nên 2 khẳng định 217 ko phải số nguyên tố. 211 là tên lừa đảo

ví dụ n = 13 (số to quá viết ko nổi :V)
n-1 = 12 = 22 * 3. d = 3, s = 2. Với số a thuộc [2, 11], ta cần kiểm tra

  • a3 ≡ 1 (mod 13), hoặc
  • a3 ≡ -1 (mod 13), hoặc
  • a2^1 * 3 = a6 ≡ -1 (mod 13)
    với a = 2, ta có 26 = 64 = 5 * 13 - 1 nên 26 ≡ -1 (mod 13), 2 “chém gió” 13 là số nguyên tố
    với a = 3, ta có 33 = 27 = 2 * 13 + 1 nên 33 ≡ 1 (mod 13), 3 cũng “chém gió” 13 là số nguyên tố
    v.v… với mọi a thuộc [2,11] đều chém 13 là số nguyên tố, vậy chắc chắn nó là nguyên tố :V

nếu n lớn thì kiểm tra ko hết a nổi, phải dựa vào tỉ lệ phần trăm mấy đứa a chém gió rồi tính xác suất n là số nguyên tố. Ở đây tỉ lệ mấy đứa chém gió nếu chọn a ngẫu nhiên lên tới 75% (ko phải 25%, nhầm :V), hay nếu có 20 số a chọn ngẫu nhiên chém gió n là số nguyên tố thì chỉ có 2-40 khả năng n ko phải là số nguyên tố.

4 Likes

Bổ đề Euclid :slight_smile: còn 1/4 là upper bound error, tức là ít nhất 75% số base (a) sẽ kick hợp số ra, và tăng dần khi hợp số lớn hơn. Số lớn (cho RSA key) có thể quay 10 vòng.

2 Likes

quay 10 vòng ít quá, 2^-20 là 1 phần triệu, mà key RSA xài cho cả tỉ trang web mỗi năm, 1 phần triệu thành ra mỗi năm có 1000 key ko phải nguyên tố à, sao được

20 vòng ít nhất :V

kiểm tra lẹ lắm, số 1024 bit dưới 1 giây là ra rồi.

2 Likes

edit: ehehe lười đọc lười cập nhật thì phán 20 lần test, lạc hậu quá rồi :flushed:

“Average Case Error Estimates for the Strong Probable Prime Test” từ trang 177 tới 194 khẳng định khả năng với số 1024 bit chỉ cần 6 lần chọn a ngẫu nhiên mà pass M-R test là tỉ lệ < 2-128 rồi :V

p(k, t) < k3/2 2t t-1/2 42 - sqrt(tk) với t = 2, k ≥ 88 hoặc 3 ≤ t ≤ k/9, k ≥ 21

thế số vô k = 1024, t = 6 ta có p(1024, 6) < 10243/2 26 6-1/2 42 - sqrt(1024*6) = 8.81x10-41 < 2-128

FIPS 186-4 table C.2, C.3 cũng chỉ rõ số round tối thiểu cho k=1024 chỉ cần là 4 hoặc 5. Mấy ông NIST tính theo cái thuật toán gì phức tạp ở appendix F.1, đọc hoa cả mắt =)

3 Likes
83% thành viên diễn đàn không hỏi bài tập, còn bạn thì sao?