51nod1186(Miller-Rabin)

题目链接：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1186

题意：中文题目诶～

思路：miller_rabin模板 (详情可参见: http://blog.csdn.net/s031302306/article/details/51606018)

注意这里的 n 为 1e30，需要用字符串处理

代码：

#include <string.h>
#include <cstdlib>

#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){//计算a0*b0%p
    ans = 0;
    bnum tmp = a0, b = b0;
    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){//计算a0^b0%p
    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(void){
    bnum x;
    x.read();
    puts(MillerRabinTest(x, 5) ? "Yes" : "No");
    return 0;
}

对于long long范围：

#include <stdio.h>
#include <iostream>
#include <cstdlib>
#define ll long long
using namespace std;

//利用二进制计算a*b%mod
ll multi_mod(ll a, ll b, ll mod){
    ll ans = 0LL;
    a %= mod;
    while(b){
        if (b & 1LL) ans = (ans+a)%mod;
        b >>= 1LL;
        a = (a+a)%mod;
    }
    return ans;
}

//计算a^b%mod
ll power_mod(ll a, ll b, ll mod){
    ll ans = 1LL;
    a %= mod;
    while(b){
        if(b & 1LL) ans = multi_mod(ans, a, mod);
        b >>= 1LL;
        a = multi_mod(a, a, mod);
    }
    return ans;
}

//Miller-Rabin测试,测试n是否为素数
bool miller_rabin(ll n, int repeat){
    if (2LL == n || 3LL == n) return true;
    if (n%2==0 || n%3==0 || n%5==0) return false;

    //将n分解为2^s*d
    ll d = n - 1LL;
    int s = 0;
    while(!(d & 1LL)){
        s++;
        d >>= 1LL;
    }
    // srand((unsigned)time(0));
    for(int i=0; i<repeat; ++i){//重复repeat次
        ll a = rand() % (n - 3) + 2;//取一个随机数,[2,n-1)
        ll x = power_mod(a, d, n);
        ll y = 0LL;
        for(int j=0; j<s; ++j){
            y = multi_mod(x, x, n);
            if (1LL == y && 1LL != x && n-1LL != x) return false;
            x = y;
        }
        if (1LL != y) return false;
    }
    return true;
}

int main(void){
    ll n;
    cin >> n;
    if(miller_rabin(n, 5)) cout << "Yes" << endl;
    else cout << "No" << endl;
    return 0;
}