只要求出梅森素数,就可以求出一个完满数。代码如下:(以后有时间会用重定义操作符的方法让代码更简洁明了的)
#include <iostream> #include <cmath> #include <string> #define MLEN 101 #define MAX(a, b) ((a) > (b) ? (a) : (b)) using namespace std; bool isABiggerB(int a[], int b[]);//a是否比b大 void add(int a[], int b[], int ans[]);//加数,加数,和 void minus(int a[], int b[], int ans[]);//被减数,减数,差 void times(int a[], int b[], int c[]);//乘数,乘数,积 void division(int a[], int b[], int qu[], int rem[]);//被除数,除数,商,余数 void fastPower(int x, int n, int a[]);//底数,指数,幂 void fp(int &x, int p, int a[], int b[]);//fastPower的递规函数 bool isPrime(int a[]);//a是否为素数 bool theSame(int a[], int b[]);//a,b是否相等 //以下所有高精度统一从第1格开始存储数据,第0格存储该数的长度 bool isABiggerB(int a[], int b[]) { if (a[0] > b[0]) return true; if (a[0] < b[0]) return false; for (int i = a[0]; i >= 1; --i) { if (a[i] > b[i]) return true; else if (a[i] < b[i]) return false; } return true; } void add(int a[], int b[], int ans[]) { int jw = 0; int c[MLEN]; memset(c, 0, MLEN * sizeof(int)); c[0] = MAX(a[0], b[0]); for (int i = 1; i <= c[0]; ++i) { jw = jw + (i <= a[0] ? a[i] : 0) + (i <= b[0] ? b[i] : 0); c[i] = jw % 10; jw = jw / 10; } if (jw > 0) c[c[0] + 1] = jw, c[0]++; for (int i = 0; i <= MLEN; ++i) ans[i] = c[i]; } void minus(int a[], int b[], int ans[]) { int c[MLEN]; memset(c, 0, MLEN * sizeof(int)); for (int i = 1; i <= a[0]; ++i) { int bvalue = i > b[0] ? 0 : b[i]; c[i] = a[i] - bvalue; if (c[i] < 0) { c[i] += 10; --a[i + 1]; } } while (!c[a[0]] && a[0] > 1) --a[0]; c[0] = a[0]; for (int i = 0; i <= MLEN; ++i) ans[i] = c[i]; } void times(int a[], int b[], int ans[]) { int c[MLEN]; memset(c, 0, MLEN * sizeof(int)); for (int i = 1; i <= a[0]; ++i) { for (int j = 1; j <= a[0]; ++j) { c[i + j - 1] += a[i] * b[j]; c[i + j] += c[i + j - 1] / 10; c[i + j - 1] %= 10; } } c[0] = a[0] + b[0]; while (c[c[0]] == 0 && c[0] > 1) --c[0]; for (int i = 0; i <= MLEN; ++i) ans[i] = c[i]; } void division(int a[], int b[], int qu[], int rem[]) { for (int i = a[0]; i >= 1; --i) { for (int j = rem[0]; j >= 1; --j) rem[j + 1] = rem[j]; rem[1] = a[i]; if (rem[rem[0] + 1] > 0) ++rem[0]; while (isABiggerB(rem, b)) { minus(rem, b, rem); ++qu[i]; } } qu[0] = a[0]; while (qu[0] > 1 && qu[qu[0]] == 0) --qu[0]; } void fastPower(int x, int n, int a[]) { int b[MLEN]; a[0] = floor(n * log((double) x) / log((double) 10) + 1); memset(b, 0, MLEN * sizeof(int)); a[1] = 1; fp(x, n, a, b); } void fp(int &x, int p, int a[], int b[]) { if (p == 0) return; fp(x, p / 2, a, b); for (int i = 1; i <= a[0]; ++i) { for (int j = 1; j <= a[0]; ++j) { if (p % 2 == 1) b[i + j - 1] += a[i] * a[j] * x; else b[i + j - 1] += a[i] * a[j]; } } for (int i = 1; i <= a[0]; ++i) { b[i + 1] += b[i] / 10; a[i] = b[i] % 10; } memset(b, 0, MLEN * sizeof(int)); } bool theSame(int a[], int b[]) { if (a[0] != b[0]) return false; for (int i = 1; i <= a[0]; ++i) { if (a[i] != b[i]) return false; } return true; } bool isPrime(int n[]) { int zero[2], one[2];//表示0和1 int judge[MLEN];//存储余数 int i[MLEN];//计数器 int temp[MLEN];//i的平方 zero[0] = 1, zero[1] = 0; one[0] = 1, one[1] = 1; memset(judge, 0, MLEN * sizeof(int)); memset(i, 0, MLEN * sizeof(int)); memset(temp, 0, MLEN * sizeof(int)); i[0] = 1, i[1] = 2;//i = 2 do { memset(temp, 0, MLEN * sizeof(int)); memset(judge, 0, MLEN * sizeof(int)); times(i, i, temp); if (isABiggerB(temp, n)) { break; } memset(temp, 0, MLEN * sizeof(int)); division(n, i, temp, judge); if (theSame(zero, judge)) return false; add(i, one, i); } while (true);//同下面的“bool isPrime(int n)”逻辑 return true; } bool isPrime(int n) { if (n < 2) return false; int i; for (i = 2; i * i <= n; ++i) if (n % i == 0) return false; return true; } int main() { int one[2]; one[0] = 1, one[1] = 1; int m[MLEN], nm[MLEN]; memset(m, 0, MLEN * sizeof(int)); memset(nm, 0, MLEN * sizeof(int)); for (int p = 1; p < 30; ++p) { if (isPrime(p)) { fastPower(2, p, m); minus(m, one, m); if (isPrime(m)) { fastPower(2, p - 1, nm); times(m, nm, m);//由梅森素数求完满数的公式 for (int i = m[0]; i > 0; --i) cout << m[i]; cout << endl; } } memset(m, 0, MLEN * sizeof(int)); memset(nm, 0, MLEN * sizeof(int)); } system("pause"); return 0; }
如果把MAXL再调大一点的话,是可以求出接下来几个数的。
不过如果输出不出来的话,那就到网上去搜一个在线c++编译器再编译一下吧。
下面附上一个在线编译器的地址:http://codepad.org/