素数相关

来源：互联网发布：知乎一个笑话笑一辈子编辑：程序博客网时间：2024/04/27 21:11

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 9const 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

0 0