#include <NTL/ZZ.h>
#include<iostream>

using namespace std;
using namespace NTL;


//n为素数候选者,x为随机数
long witness(const ZZ& n, const ZZ& x) {
    ZZ d, y, z;
    long j, s;

    if (x == 0) return 0;

    //计算s,d，使得n-1 = 2^s * d, d是奇数:
    s = 1;
    d = n / 2;
    while (d % 2 == 0) {
        s++;
        d /= 2;
    }

    /*在NTL内置的函数中，有一个刚好满足我们的要求：
    m = n - 1;
    k = MakeOdd(m);
    */

    z = PowerMod(x, d, n); // z = x^d % n
    if (z == 1) return 0;

    j = 0;
    do {
        y = z;
        z = (y * y) % n;
        j++;
    } while (j < s && z != 1);

    return z != 1 || y != n - 1;
}

//n为待检测素数，t为检测次数
long PrimeTest(const ZZ& n, long t)
{
    if (n <= 1) return 0;

    //用2000以内的素数对n进行初筛
    PrimeSeq s;  // 生成一个素数数列
    long p;
    p = s.next();  // first prime is always 2
    while (p && p < 2000) {
        if ((n % p) == 0) return (n == p);
        p = s.next();
    }

    //Miller-Rabin法推演n的素性
    ZZ x;

    for (long i = 0; i < t; i++) {
        x = RandomBnd(n); // 随机数 between 0 and n-1
        if (witness(n, x))//有凭证
            return 0;
    }

    return 1;
}

int main()
{
    ZZ n;
    long t;
    cout << "请输入Miller-Rabin待检测的n: ";
    cin >> n;
    cout << "请输入Miller-Rabin检测次数t：";
    cin >> t;
    if (PrimeTest(n, t))
        cout << "n是大概率素数\n";
    else
        cout << "n是合数\n";
}