#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); }