表为平方和初级版

今天我们研究的问题叫做表为平方和问题，简单来说就是对于丢番图方程

来说，求满足条件的整数解。当然，根据

的范围不同，所使用的方法也可能会不一样，接下来将会以几个题为例进行深度剖析。

题目：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;
            else
                node[_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;
}