素数相关

ACM模版

判断小于MAXN的数是不是素数

/* * 素数筛选，判断小于MAXN的数是不是素数 * notprime是一张表，false表示是素数，true表示不是 */

const int MAXN = 1000010;
bool notprime[MAXN];

void init()
{
    memset(notprime, false, sizeof(notprime));
    notprime[0] = notprime[1] = true;
    for (int i = 2; i < MAXN; i++)
    {
        if (!notprime[i])
        {
            if (i > MAXN / i)   // 阻止后边i * i溢出（或者i,j用long long)
            {
                continue;
            }
            // 直接从i * i开始就可以，小于i倍的已经筛选过了
            for (int j = i * i; j < MAXN; j += i)
            {
                notprime[j] = true;
            }
        }
    }
}

查找出小于等于MAXN的素数（生成连续素数表）

/* * 素数筛选，查找出小于等于MAXN的素数 * prime[0]存素数的个数 */

const int MAXN = 100000;
int prime[MAXN + 1];

void getPrime()
{
    memset(prime, 0, sizeof(prime));
    for (int i = 2; i <= MAXN; i++)
    {
        if (!prime[i])
        {
            prime[++prime[0]] = i;
        }
        for (int j = 1; j <= prime[0] && prime[j] <= MAXN / i; j++)
        {
            prime[prime[j] * i] = 1;
            if (i % prime[j] == 0)
            {
                break;
            }
        }
    }
}

随机素数测试

/* * 随机素数测试（伪素数原理） * CALL: bool res = miller(n); * 快速测试n是否满足素数的“必要”条件，出错概率极低 * 对于任意奇数n > 2和正整数s，算法出错概率≤2^(-s) */

int witness(int a, int n)
{
    int x, d = 1;
    int i = ceil(log(n - 1.0) / log(2.0)) - 1;
    for (; i >= 0; i--)
    {
        x = d;
        d = (d * d) % n;
        if (d == 1 && x != 1 && x != n - 1)
        {
            return 1;
        }
        if (((n - 1) & (1 << i)) > 0)
        {
            d = (d * a) % n;
        }
    }
    return (d == 1 ? 0 : 1);
}

int miller(int n, int s = 50)
{
    if (n == 2)     // 质数返回1
    {
        return 1;
    }
    if (n % 2 == 0) // 偶数返回0
    {
        return 0;
    }
    int j, a;
    for (j = 0; j < a; j++)
    {
        a = rand() * (n - 2) / RAND_MAX + 1;
        // rand()只能随机产生[0, RAND_MAX)内的整数
        // 而且这个RAND_MAX只有32768直接%n的话是永远
        // 也产生不了[RAND_MAX, n)之间的数
        if (witness(a, n))
        {
            return 0;
        }
    }
    return 1;
}

大数素数测试

#define MAXL 4
#define M10 1000000000
#define Z10 9

const int zero[MAXL - 1] = {
  0};

struct bnum
{
    int data[MAXL]; // 断成每截9个长度

    // 读取字符串并转存
    void read()
    {
        memset(data, 0, sizeof(data));
        char buf[32];
        scanf("%s", buf);
        int len = (int)strlen(buf);
        int i = 0, k;
        while (len >= Z10)
        {
            for (k = len - Z10; k < len; ++k)
            {
                data[i] = data[i] * 10 + buf[k] - '0';
            }
            ++i;
            len -= Z10;
        }
        if (len > 0)
        {
            for (k = 0; k < len; ++k)
            {
                data[i] = data[i] * 10 + buf[k] - '0';
            }
        }
    }

    bool operator == (const bnum &x)
    {
        return memcmp(data, x.data, sizeof(data)) == 0;
    }

    bnum & operator = (const int x)
    {
        memset(data, 0, sizeof(data));
        data[0] = x;
        return *this;
    }

    bnum operator + (const bnum &x)
    {
        int i, carry = 0;
        bnum ans;
        for (i = 0; i < MAXL; ++i)
        {
            ans.data[i] = data[i] + x.data[i] + carry;
            carry = ans.data[i] / M10;
            ans.data[i] %= M10;
        }
        return  ans;
    }

    bnum operator - (const bnum &x)
    {
        int i, carry = 0;
        bnum ans;
        for (i = 0; i < MAXL; ++i)
        {
            ans.data[i] = data[i] - x.data[i] - carry;
            if (ans.data[i] < 0)
            {
                ans.data[i] += M10;
                carry = 1;
            }
            else
            {
                carry = 0;
            }
        }
        return ans;
    }

    // assume *this < x * 2
    bnum operator % (const bnum &x)
    {
        int i;
        for (i = MAXL - 1; i >= 0; --i)
        {
            if (data[i] < x.data[i])
            {
                return *this;
            }
            else if (data[i] > x.data[i])
            {
                break;
            }
        }
        return ((*this) - x);
    }

    bnum & div2()
    {
        int  i, carry = 0, tmp;
        for (i = MAXL - 1; i >= 0; --i)
        {
            tmp = data[i] & 1;
            data[i] = (data[i] + carry) >> 1;
            carry = tmp * M10;
        }
        return *this;
    }

    bool is_odd()
    {
        return (data[0] & 1) == 1;
    }

    bool is_zero()
    {
        for (int i = 0; i < MAXL; ++i)
        {
            if (data[i])
            {
                return false;
            }
        }
        return true;
    }
};

void mulmod(bnum &a0, bnum &b0, bnum &p, bnum &ans)
{
    bnum tmp = a0, b = b0;
    ans = 0;
    while (!b.is_zero())
    {
        if (b.is_odd())
        {
            ans = (ans + tmp) % p;
        }
        tmp = (tmp + tmp) % p;
        b.div2();
    }
}

void powmod(bnum &a0, bnum &b0, bnum &p, bnum &ans)
{
    bnum tmp = a0, b = b0;
    ans = 1;
    while (!b.is_zero())
    {
        if (b.is_odd())
        {
            mulmod(ans, tmp, p, ans);
        }
        mulmod(tmp, tmp, p, tmp);
        b.div2();
    }
}

bool MillerRabinTest(bnum &p, int iter)
{
    int i, small = 0, j, d = 0;
    for (i = 1; i < MAXL; ++i)
    {
        if (p.data[i])
        {
            break;
        }
    }
    if (i == MAXL)
    {
        // small integer test
        if (p.data[0] < 2)
        {
            return  false;
        }
        if (p.data[0] == 2)
        {
            return  true;
        }
        small = 1;
    }
    if (!p.is_odd())
    {
        return false;   // even number
    }
    bnum a, s, m, one, pd1;
    one = 1;
    s = pd1 = p - one;
    while (!s.is_odd())
    {
        s.div2();
        ++d;
    }

    for (i = 0; i < iter; ++i)
    {
        a = rand();
        if (small)
        {
            a.data[0] = a.data[0] % (p.data[0] - 1) + 1;
        }
        else
        {
            a.data[1] = a.data[0] / M10;
            a.data[0] %= M10;
        }
        if (a == one)
        {
            continue;
        }

        powmod(a, s, p, m);

        for (j = 0; j < d && !(m == one) && !(m == pd1); ++j)
        {
            mulmod(m, m, p, m);
        }
        if (!(m == pd1) && j > 0)
        {
            return false;
        }
    }
    return true;
}

int main()
{
    bnum x;

    x.read();
    puts(MillerRabinTest(x, 5) ? "Yes" : "No");

    return 0;
}

应用例题

POJ 2689 Prime Distance
POJ 2689 Prime Test