表为平方和初级版
今天我們研究的問題叫做表為平方和問題，簡單來說就是對于丟番圖方程來說，求滿足條件的整數
解。當然，根據的范圍不同，所使用的方法也可能會不一樣，接下來將會以幾個題為例進行深度剖析。
?
?
題目：http://www.lydsy.com/JudgeOnline/problem.php?id=1041
?
題意：給定一個圓的方程，其中有，求在這個圓上有多少個整數點。
?
分析：這個問題比較簡單，因為對于丟番圖方程來說，整數解的個數可以通過如下方法計算
?
???? ?
?
?????其中，滿足
?
?????（1）當為偶數時，有
???? （2）當為奇數時，有
?
?????那么，本題的方法就很明確了，先對進行大數分解，然后搜出所有因子，再枚舉所有因子判斷累加即可。由
???? 于比較簡單，省略代碼。
?
?
題目：http://acm.timus.ru/problem.aspx?num=1593
?
題意：給定一個正整數，其中，求該正整數最少能夠使用多少個正整數的平方和來表示。
?
分析：先來認識幾個重要的定理，如下
?
????（1）費馬平方和定理
???????奇質數能表示為兩個平方數之和的充分必要條件是該素數被4除余1
?
????（2）費馬平方和定理的拓展定理
???????正整數能表示為兩平方數之和的充要條件是在它的標準分解式中，形如素因子的指數是偶數
?
????（3）Brahmagupta–Fibonacci identity
????????如果兩個整數都能表示為兩個平方數之和，則它們的積也能表示為兩個平方數之和。公式及拓展公式為
?
?????? ?
?
??????? 從這個定理可以看出：如果不能表示為三個數的平方和，那么也就不能表示為兩個數的平方和。
?
????（4）四平方和定理
????????每個正整數都可以表示成四個整數的平方數之和
?
????（5）表為3個數的平方和條件
????????正整數能表示為三個數的平方和的充要條件是不能表示成的形式，其中和為非負
????????整數。
?
?????通過上面的幾個重要的定理，就可以來做這題了。
?
代碼：
#include <iostream> #include <string.h> #include <stdio.h> #include <math.h>using namespace std; typedef long long LL;int Work(LL n) {while(n % 4 == 0) n >>= 2;if(n % 8 == 7) return 4;LL i = 8, t = 9;while(t <= n){while(n % t == 0) n /= t;i += 8;t += i;}if(n == 1) return 1;if(n % 2 == 0) n >>= 1;if(n % 4 == 3) return 3;LL k = 3;while(k * k <= n){if(n % k == 0) return 3;k += 4;}return 2; }int main() {LL n;while(cin>>n)cout<<Work(n)<<endl;return 0; }
上面第13行至第19行的作用是消去中所有素數的偶次冪，這里為，每次增加2，而每次相鄰兩項滿
足，所以每次增加8。
?
接下來，開始進入我們今天的重點，求丟番圖方程的解。這里很大，范圍是。
?
題目：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1171
?
分析：通過Brahmagupta–Fibonacci identity知道，如果兩個整數都能表示為兩個平方數之和，則它們的積也
?????能表示為兩個平方數之和，即
?
?????
?
?????那么，我們可以先對素因子分解，采用Pollard-rho分解法，然后對于每個素數，求出的
?????解，然后進行合并即可。而求解丟番圖方程（其中為素數）的方法如下
?
???? （1）首先判斷奇素數是否滿足條件，如果是才有解，并且解唯一，否則無解。
???? （2）找出二次同余方程的一個解
???? （3）對奇素數和進行歐幾里德輾轉相除運算，記錄每次的余數為（注意第一次的余數為），當滿
? ? ? ? ?足時結束，此時的就是丟番圖方程中的，進而通過得到。
?
代碼：
#include <iostream> #include <stdlib.h> #include <string.h> #include <algorithm> #include <stdio.h> #include <math.h> #include <set>const int Times = 10; const int N = 65;using namespace std; typedef long long LL;LL ct, cnt; LL fac[N], num[N];struct T {LL p, d; }; LL w;struct Pair {LL x, y;bool operator < (const Pair &a) const{if(a.x == x) return a.y > y;return a.x > x;} }; Pair node[N]; set<Pair> Set; set<Pair> Mem;int sign[2][2] = {{-1, 1}, {1, -1}};LL gcd(LL a, LL b) {return b? gcd(b, a % b) : a; }LL multi(LL a, LL b, LL m) {LL ans = 0;a %= m;while(b){if(b & 1){ans = (ans + a) % m;b--;}b >>= 1;a = (a + a) % m;}return ans; }LL quick_mod(LL a, LL b, LL m) {LL ans = 1;a %= m;while(b){if(b & 1){ans = multi(ans, a, m);b--;}b >>= 1;a = multi(a, a, m);}return ans; }bool Miller_Rabin(LL n) {if(n == 2) return true;if(n < 2 || !(n & 1)) return false;LL m = n - 1;int k = 0;while((m & 1) == 0){k++;m >>= 1;}for(int i=0; i<Times; i++){LL a = rand() % (n - 1) + 1;LL x = quick_mod(a, m, n);LL y = 0;for(int j=0; j<k; j++){y = multi(x, x, n);if(y == 1 && x != 1 && x != n - 1) return false;x = y;}if(y != 1) return false;}return true; }LL pollard_rho(LL n, LL c) {LL i = 1, k = 2;LL x = rand() % (n - 1) + 1;LL y = x;while(true){i++;x = (multi(x, x, n) + c) % n;LL d = gcd((y - x + n) % n, n);if(1 < d && d < n) return d;if(y == x) return n;if(i == k){y = x;k <<= 1;}} }void find(LL n, int c) {if(n == 1) return;if(Miller_Rabin(n)){fac[ct++] = n;return ;}LL p = n;LL k = c;while(p >= n) p = pollard_rho(p, c--);find(p, k);find(n / p, k); }void Partion(LL n) {ct = 0;find(n, 120);sort(fac, fac + ct);num[0] = 1;int k = 1;for(int i=1; i<ct; i++){if(fac[i] == fac[i-1])++num[k-1];else{num[k] = 1;fac[k++] = fac[i];}}cnt = k; }bool Filter() {for(int i = 0; i < cnt; i++){if(fac[i] % 4 == 3 && (num[i] & 1))return true;}return false; }T multi_er(T a, T b, LL m) {T ans;ans.p = (multi(a.p, b.p, m) + multi(multi(a.d, b.d, m), w, m)) % m;ans.d = (multi(a.p, b.d, m) + multi(a.d, b.p, m)) % m;return ans; }T power(T a, LL b, LL m) {T ans;ans.p = 1;ans.d = 0;while(b){if(b & 1){ans = multi_er(ans, a, m);b--;}b >>= 1;a = multi_er(a, a, m);}return ans; }LL _power(LL a, LL b) {LL ans = 1;while(b){if(b & 1){ans = ans * a;b--;}b >>= 1;a = a * a;}return ans; }LL Legendre(LL a, LL p) {return quick_mod(a, (p - 1) >> 1, p); }LL Work(LL n, LL p) {LL a = -1, t;while(1){a = rand() % p;t = a * a - n;w = (t % p + p) % p;if(Legendre(w, p) + 1 == p) break;}T tmp;tmp.p = a;tmp.d = 1;T ans = power(tmp, (p + 1) >> 1, p);return ans.p; }void Solve(LL n, LL p, LL &_x, LL &_y) {LL x = Work(n, p);if(x >= p - x)x = p - x;LL t = p;LL r = x;LL limit = (LL)sqrt(1.0 * p);while(r > limit){x = r;r = t % x;t = x;}_x = r;_y = (LL)sqrt((p - r * r));if(_x > _y) swap(_x, _y); }void getData(Pair node[], int &_cnt) {for(int i = 0; i < cnt; i++){if(fac[i] == 2){LL res = _power(fac[i], num[i]>>1);if(num[i] & 1)node[_cnt].x = res;elsenode[_cnt].x = 0;node[_cnt].y = res;_cnt++;continue;}if(fac[i] % 4 == 3){LL res = _power(fac[i], num[i]>>1);node[_cnt].x = 0;node[_cnt].y = res;_cnt++;continue;}Solve(-1, fac[i], node[_cnt].x, node[_cnt].y);_cnt++;for(int j = 1; j < num[i]; j++){node[_cnt].x = node[_cnt - 1].x;node[_cnt].y = node[_cnt - 1].y;_cnt++;}} }void dfs(int dept, int cnt, Pair ans) {if(dept == cnt - 1){Set.insert(ans);return;}for(int i = 0; i < 2; i++){Pair _ans;_ans.x = ans.x * node[dept + 1].x + sign[i][0] * ans.y * node[dept + 1].y;_ans.y = ans.x * node[dept + 1].y + sign[i][1] * ans.y * node[dept + 1].x;if(_ans.x < 0) _ans.x = -_ans.x;if(_ans.y < 0) _ans.y = -_ans.y;if(_ans.x > _ans.y) swap(_ans.x, _ans.y);if(Mem.find(_ans) != Mem.end()) return;Mem.insert(_ans);dfs(dept + 1, cnt, _ans);} }void Merge(Pair node[], int cnt) {Pair ans;ans.x = node[0].x;ans.y = node[0].y;Mem.insert(ans);dfs(0, cnt, ans); }int main() {LL n;while(cin>>n){if(n == 1){puts("0 1");continue;}Partion(n);if(Filter()){puts("No solution");continue;}Set.clear();Mem.clear();int _cnt = 0;getData(node, _cnt);Merge(node, _cnt);set<Pair>::iterator it;for(it = Set.begin(); it != Set.end(); it++)cout<<it->x<<" "<<it->y<<endl;}return 0; }
?
總結
- 上一篇： 关于丢番图方程x^2-dy^2=-1
- 下一篇： 2013年东北赛B题(数位DP)