数论

#pragma warning (disable: 4786)
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<deque>
#include<map>
#include<algorithm>
using namespace std;

#ifdef _WIN32
typedef unsigned __int64 UINT64;
#elif
typedef unsigned long long UINT64;
#endif

//1000以内素数表
UINT64 gPrimeTable1000[] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
    101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,
    211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,
    337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,
    461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,
    601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,
    739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,
    881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997};

/******************************************************************************
                         数论类
******************************************************************************/
class Number
{
public:
    Number(UINT64 n=0) : num(n) {} //构造函数
    void print (ostream& oo) const; //打印
    static Number rand(); //随机数

    //基本运算
    bool    operator == (const Number& n) const {return num == n.num;}
    bool    operator != (const Number& n) const {return num != n.num;}
    bool    operator >  (const Number& n) const {return num > n.num;}
    bool    operator <  (const Number& n) const {return num < n.num;}
    bool    operator >= (const Number& n) const {return num >= n.num;}
    bool    operator <= (const Number& n) const {return num <= n.num;}
    Number  operator +  (const Number& n) const {return num + n.num;}
    Number  operator -  (const Number& n) const {return num - n.num;}
    Number  operator *  (const Number& n) const {return num * n.num;}
    Number  operator /  (const Number& n) const {return num / n.num;}
    Number  operator %  (const Number& n) const {return num % n.num;}
    Number& operator += (const Number& n) {num += n.num; return *this;}
    Number& operator -= (const Number& n) {num -= n.num; return *this;}
    Number& operator *= (const Number& n) {num *= n.num; return *this;}
    Number& operator /= (const Number& n) {num /= n.num; return *this;}
    Number& operator %= (const Number& n) {num %= n.num; return *this;}

    bool isPrime() const; //是否素数
    Number power(unsigned int index) const; //幂
    Number powerModule(Number index, Number mod) const; //幂余
    void factorization(map<Number,int>& factor) const; //因子分解
    Number primeNumber(Number n) const;  //n以内素数个数
    Number eulerFunc() const; //Euler函数φ
    Number mobiusFunc() const; //M?bius函数μ
    int legendre(Number p) const; //Legendre符号, 二次剩余, p必须是素数
    static Number gcd(Number a, Number b); //最大公约数
    static Number gcd(Number a[], int n); //一组数的最大公约数
    static void primeTable(int n, deque<Number>& table); //素数表
    static Number fabonacci(int n); //Fabonacci数
    static Number combinatory(unsigned int n, unsigned int m);
    static Number permutation(unsigned int n, unsigned int m);

    Number RSA_encryption(Number r, Number e);
    Number RSA_decryption(Number r, Number d);
private:
    Number powerForIsPrime(Number a, Number p) const; //isPrime的辅助函数,幂余的时候加二次探测
    Number pollard() const; //Pollard法求一个因子
public:
    UINT64 num;
};

ostream& operator << (ostream& oo,Number& n) //重载插入运算符
{
    n.print(oo);
    return oo;
}

void Number::print(ostream& oo) const
{
#ifdef _WIN32
    printf("%I64u",num);
#elif
    printf("%llu",num);
#endif
}

Number Number::rand()
{
    return Number(UINT64(::rand()) * ::rand() * ::rand() * ::rand());
}

Number Number::power(unsigned int index) const
{
    Number ans = 1;
    Number pow = num;

    while (index)
    {
        if (index & 1)
            ans.num *= pow.num;
        pow.num *= pow.num;
        index >>= 1;
    }

    return ans;
}

Number Number::powerModule(Number index, Number mod) const
{
    Number r = 1;
    Number pow = num;

    while (index.num)
    {
        if (index.num & 1)
            r.num = r.num * pow.num % mod.num;
        pow.num = pow.num * pow.num % mod.num;
        index.num >>= 1;
    }

    return r;
}

Number Number::powerForIsPrime(Number a, Number p) const
{
    if (p == 0)
        return 1;
    Number x = powerForIsPrime(a, p/2);
    Number result = x*x % num; //二次探测
    if (result == 1 && x != 1 && x != num-1) //若num为素数, 则方程x^2≡1的解为x=1或p-1
        return 0; //是合数
    if (p%2 == 1)
        result = result * a % num;
    return result;
}

bool Number::isPrime() const
{
    if (num <= 1000)
        return find(gPrimeTable1000,gPrimeTable1000+168,num) != gPrimeTable1000+168;

    for (int time=1; time<10; time++)
    {
        Number a = rand()%(num-3) + 2;
        Number result = powerForIsPrime(a, num-1); //Fermart小定理
        if (result == 0 || result != 1)
            return false;
    }
    return true;
}

Number Number::pollard() const
{
    for (int j=0; j<168; j++)
        if (num % gPrimeTable1000[j] == 0)
            return gPrimeTable1000[j];

    Number x = rand() % num;
    Number y = x;
    int k = 2;
    for (int i=2; i<1024; i++)
    {
        x = (x*x - 1) % num;
        Number d = gcd(y-x, num);
        if (d > 1 && d < num)
            return d;
        if (i == k)
        {
            y = x;
            k *= 2;
        }
    }
    return 1;
}

void Number::factorization(map<Number,int>& factor) const
{
    if (isPrime())
    {
        factor[num]++;
        return;
    }

    Number d;
    for (int t=1; t<=10; t++)
    {
        d = pollard();
        if (d.num != 1)
            break;
    }
    if (d.num == 1)
        return;

    d.factorization(factor);
    (*this/d).factorization(factor);
}

Number Number::eulerFunc() const
{
    map<Number,int> factor;
    factorization(factor);
    Number ans = num;
    for (map<Number,int>::const_iterator it=factor.begin(); it!=factor.end(); ++it)
        ans = ans * (it->first-1) / it->first;
    return ans;
}

Number Number::mobiusFunc() const
{
    if (num == 1)
        return 1;

    map<Number,int> factor;
    factorization(factor);
    for (map<Number,int>::const_iterator it=factor.begin(); it!=factor.end(); ++it)
        if (it->second > 1)
            return 0;

    if(factor.size() % 2 == 0)
        return 1;
    else
        return -1;
}

int Number::legendre(Number p) const
{
    if (num % p.num == 0)
        return 0;
    else if (powerModule((p-1)/2, p) == 1)
        return 1;
    else
        return -1;
}

Number Number::gcd(Number a, Number b)
{
    UINT64 r;
    while (b != 0)
    {
        r = a.num % b.num;
        a.num = b.num;
        b.num  = r;
    }
    return a;
}

Number Number::gcd(Number a[], int n)
{
    if (n < 2)
        return 1;
    Number d = gcd(a[0], a[1]);
    for (int i=2; i<n; i++)
        d = gcd(d, a[i]);
    return d;
}

void Number::primeTable(int n, deque<Number>& table)
{
    for (int i=3; i<n; i+=2)
    {
        int j;
        for (j=0; j<table.size(); j++)
            if (i % table[j].num == 0)
                break;
        if (j == table.size())
            table.push_back(i);
    }
    table.push_front(2);
}

Number Number::combinatory(unsigned int n, unsigned int m)
{
    if (m > n)
        return 0;
    UINT64 ans = 1;
    for (int i=0; i<m; i++)
        ans = (ans * (n-i)) / (i+1);
    return ans;
}

Number Number::permutation(unsigned int n, unsigned int m)
{
    if (m > n)
        return 0;
    UINT64 ans = 1;
    for (int i=0; i<m; i++)
        ans *= (n-i);
    return ans;
}

//RSA加密
//r=p*q为两个大素数的乘积, e与(p-1)*(q-1)互素为明文
Number Number::RSA_encryption(Number r, Number e)
{
    return powerModule(e,r);
}

//RSA解密
//r=p*q为两个大素数的乘积, d为密文, d*e≡1 (mod (p-1)*(q-1))
Number Number::RSA_decryption(Number r, Number d)
{
    return powerModule(d,r);
}