关于方程x^2+y^2=n的解的问题

来源：互联网发布：辽宁联通网络测速编辑：程序博客网时间：2024/04/29 22:34

关于方程x^2+y^2=n的解的问题

分类：数论2013-01-20 16:46 75人阅读评论(0) 收藏举报

一.有解的条件：

定理1 若正整数n可以表示成两个整数的平方之和，则在它的标准分解式n=(p1^a1)*(p2^a2)*...*(pr^ar)中，形如4k+3的素因数的指数是偶数。

定理2 正整数n 能表示成二平方数之和的充要条件，是在它的标准分解式中，形如4k+3的素因数的幂指数是偶数.

二.解的个数：

对于这个问题有人可能会想到枚举，但是枚举不是最优的，可以用下面的结论求解的个数

将正整数表示为平方数之和

分类：数论2013-01-20 23:48 96人阅读评论(0) 收藏举报

目录(?)[+]

Timus Online Judge 网站上有这么一道题目：1073. Square Country。这道题目的输入是一个不大于 60,000 的正整数，要求计算出该正整数最少能够使

用多少个正整数的平方和来表示。这道题目的时间限制是 1 秒。

四平方和定理：每个正整数都是四个平方数之和。在这个定理中，平方数是指整数(包括零)的平方。且有下述不等式成立：

所以，我们有以下C语言程序：

[cpp] view plaincopy
// http://acm.timus.ru/problem.aspx?space=1&num=1073   
#include <stdio.h>   
#include <math.h>   
    
int compute(int n)   
{   
  int i, j, k, m = 4;   
  int i0 = n / 4, i2 = n, j2, k2;   
  for (i = sqrt(n); i2 > i0; i--)   
    if ((j2 = n - (i2 = i * i)) == 0) return 1;   
    else for (j = sqrt(j2); j > 0; j--)   
      if ((k2 = n - i2 - j * j) == 0) return 2;   
      else if (k = sqrt(k2), k * k == k2 && m > 3) m = 3;   
  return m;   
}   
    
int main(void)   
{   
  int n;   
  scanf("%d", &n);   
  printf("%d", compute(n));   
  return 0;   
}   

上述程序中：

第 7 行设置 m 的初值为 4，代表一个正整数最多只需要四个平方数就可以表示了。
第 9 行开始的主循环决定第一个平方数，如果 n 刚好是平方数(第 10 行)，就直接返回 1。
第 11 行开始的内循环决定第二个平方数，如果这两个数加起来刚好等于 n (第 12 行)，就直接返回 2。
第 13 行检查 n 是否可以表示为三个平方数的和，如果是的话，就更新 m 的值为 3 。注意，此时不能直接返回 3，因为可能在后面的循环中发现 n 可以用两个平方数表示。
第 14 行返回 m 值(只可能是 3 或者 4)作为最后的答案。

上述程序在 Timus Online Judge 网站的运行时间是 0.015 秒。

上述题目有一个进一步的版本：1593. Square Country. Version 2，输入改为不大于 10¹⁵ 的正整数，时间限制还是 1 秒。上一节的程序做以下改动：

第 5 行的第 2 个 int 改为 long long
第 8 和 19 行的 int 改为 long long
第 20 行的 %d 改为 %lld

就可以适用于这道题目，但是运行结果是“Time limit exceeded”。此时，需要更好的算法。我们有以下 C 语言程序(1593.c)：

[cpp] view plaincopy
// http://acm.timus.ru/problem.aspx?space=1&num=1593   
#include <stdio.h>   
#include <math.h>   
    
int compute(long long n)   
{   
  int i, k;   
  long long i2;   
  while ((n & 3) == 0) n >>= 2;   
  if ((n & 7) == 7) return 4;   
  for (i = 8, i2 = 9; i2 <= n; i2 += i += 8)   
    while (n % i2 == 0) n /= i2;   
  if (n == 1) return 1;   
  if ((n & 1) == 0) n >>= 1;   
  if ((n & 3) == 3) return 3;   
  for (k = sqrt(n), i = 3; i <= k && n % i; i += 4) ;   
  return (i > k) ? 2 : 3;   
}   
    
int main(void)   
{   
  long long n;   
  scanf("%lld", &n);   
  printf("%d", compute(n));   
  return 0;   
}   

在上述程序中：

第 9 行消去 n 的所有值为 4 因数。
第 10 行检测 n 是否为 8m + 7 的形式，如是，直接返回 4 (请参见下节)。
第 11、12 行消去 n 的所有素因子的偶次幂(素因子 2 的偶次幂已经在第 9 行消去了)。
第 11 行中 i2 依次为：3²、5²、7²、...、t²，这是因为 (t + 1)² - (t - 1)² = 4t，每次循环 t 增加 2，所以 i 增加 4 * 2 = 8。
第 13 行，如果 n 等于 1，说明输入是个完全平方数，直接返回 1。
此时，n 的标准分解式中所有的素因子都是一次幂了。
第 14 行消去 n 的素因子 2 (如果有的话)。
第 16 行的循环中 i 从 3 开始，每次递增 4，以检查 n 是否有 4m + 3 形式的因子。
第 15 行和第 17 行根据定理 366 决定答案是两个还是三个平方之和。

这个程序在 Timus Online Judge 网站的运行时间是 0.828 秒。这道题目的最佳运行时间是 0.031 秒，不知道使用什么算法可以这么快。

上述算法的原理

定理：n ≠ (4^a)(8m + 7) 是 n 可以用三个平方数表示的一个充分必要条件

定理: 一个数 n 是两个平方之和，当且仅当在 n 的标准分解式中，它的所有形如 4m + 3 的素因子都有偶次幂

我们还有以下定理：

形如 4m + 3 的整数有形如 4m + 3 的素因子

列出平方数:

前面的 1593.c 程序只能给出答案是几个平方数之和，而对这些平方数是什么一无所知。而 1073.c 程序倒是中规中矩地想要求解这些平方数是什么，

但是从 Lagrange 定理得知最多只要四个平方数就够了，所以该程序只求解到三个平方数的情况，其余情况下答案肯定是 4 了。因此，我们将

1073.c 稍做修改，得到 1073b.c 用于列出这些平方数，如下所示：

[cpp] view plaincopy
#include <stdio.h>   
#include <stdlib.h>   
#include <math.h>   
    
static int a[5];   
    
int compute(int n)   
{   
  int i, j, k, l, m = 5;   
  int i0 = n / 4, i2 = n, j2, k2, l2;   
  for (i = sqrt(n); i2 > i0; i--)   
    if ((j2 = n - (i2 = i * i)) == 0) return a[0] = i, 1;   
    else for (j = sqrt(j2); j > 0; j--)   
      if ((k2 = n - i2 - (j2 = j * j)) == 0) return a[0] = i, a[1] = j, 2;   
      else for (k = sqrt(k2); k > 0; k--)   
        if ((l2 = n - i2 - j2 - (k2 = k * k)) == 0 && m > 3)   
          a[0] = i, a[1] =  j, a[2] = k, m = 3;   
        else if (l = sqrt(l2), l * l == l2 && m > 4)   
          a[0] = i, a[1] =  j, a[2] = k, a[3] = l, m = 4;   
  return m;   
}   
    
int main(int args, char* argv[])   
{   
  int i, n, start = 1, count = 16, k;   
  if (args > 1) start = atoi(argv[1]);   
  if (args > 2) count = atoi(argv[2]);   
  for (n = start; n < start + count; n++)   
  {   
    k = compute(n);   
    printf("%d:%6d:", k, n);   
    for (i = 0; i < k; i++) printf(" %d", a[i]);   
    puts(k > 4 ? " Error!" : "");   
  }   
  return 0;   
}   

上述程序中：

第 5 行的全局静态数组用于记录所求的平方数，数组大小为 5，而不是 4，是为了防止程序有 bug 时造成数组下标越界(第 32 行)。
第 9 行将 m 的初值从 4 改为 5，用以检测程序是否有 bug。
第 9、10 行增加了变量 l 和 l2 用于计算第四个平方数，并相应增加一层循环(第 15 行)。
第 12、14、17 和 19 行相应记录这些平方数于数组 a 中。
第 33 行在输出时检查程序是否有 bug。如果 k > 4 程序肯定有问题，违反了 Lagrange 定理。当然，k <= 4 并不意味着程序就没有问题了。:)

如果不知道 Lagrange 定理，也就是说，假设我们不知道要多少个平方数之和才够的话，这道题目看来只好用动态规划算法来求解了。

[cpp] view plaincopy
// http://acm.timus.ru/problem.aspx?space=1&num=1073   
#include <stdio.h>   
     
typedef int bool;   
     
const bool true = 1;   
const bool false = 0;   
     
bool isSquare(int n, int v, int k)   
{   
  return (n < v) ? false : (n == v) ? true : isSquare(n, v + k + 2, k + 2);   
}   
     
bool isSquareSum(int n, int m, int v, int k)   
{   
  if (n < v) return false;   
  if (m == 1) return isSquare(n, v, k);   
  return isSquareSum(n - v, m - 1, v, k) ? true : isSquareSum(n, m, v + k + 2, k + 2);   
}   
     
int compute(int n, int m)   
{   
  return isSquareSum(n, m, 1, 1) ? m : compute(n, m + 1);   
}   
     
int main(void)   
{   
  int n;   
  scanf("%d", &n);   
  printf("%d", compute(n, 1));   
  return 0;   
}   

这个程序本质上和键盘农夫园友的程序是没有区别的。分析如下：

第 9 到 12 行的 isSquare 函数判断 n 是否是不小于 v 的完全平方数。其中 k 是用于计算平方数的辅助变量。
第 14 到 19 行的 isSquareSum 函数判断 n 是否是 m 个不小于 v 的平方数之和。其中 k 是用于计算平方数的辅助变量。
第 21 到 24 行的 compute 函数计算正整数 n 最少可以表示为多少个平方数之和。

上述程序在 Timus Online Judge 网站的运行时间是 0.031 秒，而第一小节中的 1073.c 的运行时间是 0.015 秒。

如果将上述程序作如下改动：

第 9 行的前两个 int 改为 long long
第 14 行的第 1 个和第 3 个 int 改为 long long
第 21 行的第 2 个 int 改为 long long
第 28 行的 int 改为 long long
第 29 行的 %d 改为 %lld

就可以适用于“1593. Square Country. Version 2”，但是运行结果是“Crash (stack overflow)”。