uvaoj 138 Street Numbers 佩尔方程

uvaoj 138 Street Numbers 佩尔方程

这道题可以用暴力搜索的方法来做，算法效率是O(n3)，但事实上是可以用数论的方法来生成所有解的，下面将用最简单的数学语言进行介绍。

设全部的房子数为n，她住的为k，那么按照题目的要求n和k必须满足：
1 + 2 + ... + (k – 1) = (k + 1) + (k + 2) + ... + (n)    (1)
用等差数列求和公式Sn=n(a1+an)/2对(1)进行化简，得：
(k – 1)[1 + (k – 1)] / 2 = (n – k)[(k + 1) + n] / 2, 两边同乘以2
(k – 1)[1 + (k – 1)] = (n – k)[(k + 1) + n]    (2)
将(2)式再次化简，得：
2k^2 = n^2 + n    (3)
题目就转化为求出满足(3)式的最小的前10对正整数：n和k。先将上式两边配方，移项得：
(2n + 1)^2 – 2(2k)^2 = 1    (4)
令：
x = 2n + 1，y = 2k  (5)
将(5)代入(4)式，可得：
x^2 – 2y^2 = 1    (6)
不定方程(6)是一个典型的佩尔（Pell）方程，要推导求解该方程需要数论的一些基本概念，篇幅所限不再赘述。为了方便程序实现，应使用递归公式进行求解。易知方程(6)的最小一组正整数解为(x0=3, y0=2)，将题目所给的最小一组解(n=8，k=6)代入(5)式可得方程的解为(17, 12)。经验算可知（也可由下面的公式对(x0, y0)进行迭代求出），(17, 12)恰好为(6)式的第二组正整数解。注意，由3,2生成的不是(3）式的解。

至此初始值（前两组(x, y)解）计算完毕，后面的(x, y)可用迭代公式求出：
xi+1 = xix0 + nyiy0
yi+1 = xiy0 + yix0    (7)
用(7)迭代出前10组(x, y)，然后用(5)式进行反算：
n = (x – 1) / 2, k = y / 2    (8)
即可求得前10组n和k。算法实现中，注意保留前一次的x和y值。若在计算下一个y值前已完成x的迭代，那么y值的迭代将出现错误。此外由于本题没有输入，您也可以预先将答案生成，在程序中直接按格式输出结果，这样速度最快（不过有点无耻哈~）。

代码如下：

/*************************************************************************
	> File Name: 138.cpp
	> Author: gwq
	> Mail: gwq5210@qq.com 
	> Created Time: 2014年12月23日 星期二 23时24分10秒
 ************************************************************************/

#include <cmath>
#include <ctime>
#include <cctype>
#include <climits>
#include <cstdio>
#include <cstdlib>
#include <cstring>

#include <map>
#include <set>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <sstream>
#include <iostream>
#include <algorithm>

#define INF (INT_MAX / 10)
#define clr(arr, val) memset(arr, val, sizeof(arr))
#define pb push_back
#define sz(a) ((int)(a).size())

using namespace std;
typedef set<int> si;
typedef vector<int> vi;
typedef map<int, int> mii;
typedef long long ll;

const double esp = 1e-5;

int main(int argc, char *argv[])
{
	int t = 10;
	int m = 3;
	int n = 2;
	int a = 3;
	int b = 2;
	for (int i = 0; i < t; ++i) {
		int t1 = a * m + 2 * b * n;
		int t2 = a * n + b * m;
		m = t1;
		n = t2;
		int x = (m - 1) / 2;
		int y = n / 2;
		printf("%10d%10d\n", y, x);
	}

	return 0;
}

参考：

1）http://www.cnblogs.com/devymex/archive/2010/09/07/1818983.html