Lỗi khi cài đặt thuật toán Rabin - Miller

[rogp10]

Vì sao 93719377 = 4139 * 22643 lại ra true nhỉ mọi người, trên wiki ghi 2, 3, 5, 7, 11 là dư rồi. (gcc >= 5)

long long modPow(int base, int expo, long long modulo) {
	unsigned __int128 res=1;
	unsigned __int128 tmp=base;
	while(expo > 0) {
		if(expo%2) res = (res*tmp) % (unsigned __int128) modulo;
		tmp = (tmp*tmp) % (unsigned __int128) modulo;
		expo/=2;
	}
	return (long long) tmp;
}

int isPrime(long long n)
{
	if((n==23) || (n==3137) || (n==8389) || (n==157163)) return true;
	/* for(int i=0; i<PREFILTER_LIMIT; ++i) {
		if(n % primeList[i] == 0) return false;
	} */
	
	int base[] = {2, 3, 5, 7, 11};
	int tn = n-1; //chết vì câu này (dành cho bạn nào đến sau)
	int exp_2 = 0;
	while(tn % 2 == 0) { tn/=2; ++exp_2; }
	printf("%lld ", tn);
	
	for(int i=0; i<5; ++i) {
		unsigned __int128 tmp = modPow(base[i],tn,n);
		if(tmp == n-1 || tmp == 1) continue;
                 int j;
		for(int j=exp_2; j>=0; --j) { //vòng này chạy exp_2 + 1 lần => FAIL!
            //viết ngược lại từ j = 0 đến < exp_2 - 1 thôi. int j ở đây là sai.
			tmp = (tmp*tmp) % (unsigned __int128) n;
			if(tmp == 1) return false;
			else if(tmp == n-1) break;
		}
                if(j > 0 && tmp != 1) return false; //sai ở đây nữa, != n-1 mới đúng!
	}
	return true;
}

[noname00 (HK boy)]

else if(tmp == n-1) break;
    		}
    	}
    	return true;
    }

Dòng return true ở cuối phải là return false chứ nhỉ? Em chạy return false thì ra. Kể cả mấy trường hợp bác đặt ở trong if kia cũng đúng luôn.

Code của em chả khác mấy so với code của bác:

#define ll long long

const int base[] = {2, 3, 5, 7, 11, 13}; // bỏ 13 đi cũng ok

ll powM(ll bs, ll pwr, ll mod) {
	ll func = 1LL;
	int bfst=64 - __builtin_clzll(pwr);
	
	for (int i=bfst-1; i>=0; i--) {
		ll pp = func;
		func = (func * func) % mod;
		if (pwr & (1LL<<i)) func = (func * bs) % mod;
	}
	return func;
}

bool isPrime(ll n) {
	ll _n = n-1;
	int exp_2 = 0;
	while (_n % 2 == 0) _n /= 2, exp_2++;

	for (int i=0; i<sizeof(base)/sizeof(base[0]); i++) {
		if (powM(base[i], _n, n) == 1) return true;

		for (int r=0; r<exp_2; r++)
			if (powM(base[i], (1LL << r) * _n, n) == n - 1) return true;
	}
	return false;
}

[Dark.Hades (Ichiro)]

// Vote typedef thay cho define

[rogp10]

Mình đã sửa lại hàm mũ mod theo bạn và giờ đã tính đúng nhưng còn hàm chính thì dùng phép bình phương vẫn sai, ra toàn false với n đủ lớn.

[noname00 (HK boy)]

nhưng còn hàm chính thì dùng phép bình phương vẫn sai, ra toàn false với n đủ lớn.

Em vừa thử với 191507399987 và ra false, mặc dù nó đúng là số nguyên tố. Em phát hiện ra lúc nhân đã bị tràn số.

Thêm đoạn nhân mod ngoài là ngon ngay:

const ll m = 100000000LL;

ll mulM(ll a, ll b, ll mod) {
	if (a <= m && b <= m)
		return (a * b) % mod;
	if (a == 0 || b == 0)
		return 0LL;

	ll a1 = a / m, a2 = a % m;
	ll b1 = b / m, b2 = b % m;
	// a*b = (a1*m + a2)*(b1*m + b2) = a1*b1*m*m + m*(a1*b2 + a2*b1) + a2*b2

	return (mulM((a1*b1)%mod, (m*m)%mod, mod)
	+ mulM(m, mulM(a1, b2, mod) + mulM(a2, b1, mod), mod)
	+ (a2*b2)%mod) % mod;
}

Sau đó các đoạn nhân ở hàm pow đều sửa lại

func = mulM(func, func, mod);
if (pwr & (1LL<<i)) func = mulM(func, bs, mod);

// Vote typedef thay cho define

Tại quen tay code define

[tntxtnt]

đầu tiên chọn 1 cái thư viện số lớn đã rồi áp dụng, mấy kiểu int 32 64 bit này nhằm nhò gì @_@

[tntxtnt]

à thấy lỗi sai rồi, ở 2 dòng này:
else if(tmp == n-1) break;
if(j > 0 && tmp != 1) return false;
trong wiki ghi là
continue WitnessLoop
return composite
chỗ continue WitnessLoop ko đơn giản là break, vì break xong gặp ngay câu if(j > 0 && tmp != 1) return false; ở dưới là return false luôn rồi, ko continue witness loop được nữa.

phải viết là
if (tmp != n-1) return false;
vì điều kiện dể continue witness loop là tmp == n-1.
xóa cái int j kia đi ko cần nó nữa. Hơn nữa j lặp r - 1 lần, ghi vòng lặp thế kia là r + 1 lần rồi @_@

hơn nữa a = 2,3,5,7,11 là chưa chắc đâu, vì cái “bảo đảm” này là áp dụng cho cái thuật toán xác định deterministic ở dưới, ko phải cho bản rabin miller này. a ở đây phải được chọn ngẫu nhiên trong đoạn [2, n-2]

[rogp10]

Mình sẽ dùng flag ntn:

   flag = false;
   {
    //...
       if(tmp == n-1) { flag=true; break; }
       else if(tmp == 1) {flag=false; break; }
    //...
    }
    if(!flag) return false;

còn base cố định là do mình biết max rồi

[tntxtnt]

Chỗ tmp==1 return false luôn rồi cần gì flag nữa @@ hơn nữa flag tối nghĩa đặt cái tên nào khác như continuewitnessloop cho dễ hiểu

weidai11/cryptopp/blob/master/nbtheory.cpp#L151

openssl/openssl/blob/master/crypto/bn/bn_prime.c#L225

người ta tách cái witness test ra thành hàm riêng nè, khỏi cần flag

[rogp10]

Cuối cùng cũng đã tìm ra lỗi:
int tn = n-1;
vì n cỡ 64 bit thì không lí nào n-1 chỉ có 32 bit.

Cảm ơn mọi người

p/s: RM nhanh hơn chia trâu 4 lần.

[tntxtnt]

4 lần? Chậm vậy, int 32 bit à? Quất lên 2048 bit luôn cho nó chạy nhanh hơn tỉ tỉ tỉ tỉ… lần

[rogp10]

À 9 lần: 0.078 -> 0.009s so với chỉ chia cho số nguyên tố :v (thực ra dãy số nguyên tố này là kết quả trung gian trong một bài code nên không tổng quát)