P1763 埃及分数 题解

题目传送门

【-1】前言

这题的剪枝真的太妙了，很难想象巨佬是怎么独立想出来这所有的剪枝的。

本题解没有包含所有的剪枝，只选了我认为最好理解的几条剪枝。

想学习所有的剪枝的右转巨佬的题解。

【1】本题大框架：迭代加深搜索（IDDFS）

看到 $1<a<b<1000$，可以猜测分数的个数不会很多，考虑搜索。

那么怎么搜？因为我们不能确定最少分数的个数（这是我们首要要求的），可以考虑枚举分数的个数再进行搜索。

这相当于每次限制了搜索树的深度，当确定一个深度发现了解，深度就不会再加深。

这种每次限制了深度的搜索就叫迭代加深搜索（IDDFS）（又称 IDS）。

我们明显感觉到，每次限制深度为 $dep$，那么搜索深度为 $dep+1$ 时，搜索深度 $dep$ 对应的搜索树又会被重新搜一遍。浪费了一些时间。

所以这种搜索有什么优势呢？

如果我们不限定深度，那么我们可能花大力气搜出了一个答案但远没达到最优解（分数的个数超过正解），此时再及时回溯已经很难了（搜索树深度增加一层增加的节点数是指数级的），浪费的时间会更多。

const int N = 11,INF = 1e7;//N为层数上限，INF为分母最大值
void dfs(int a,int b,int x){//(a/b)表示剩余分数大小，x表示当前搜索深度
	if(x > dep)//超过了限定深度，退出
		return;
	if(a == 1){//找到了解（此时x一定等于限定的深度dep，否则在dep更小的时候就会找到解）
		if(b > tmp[x - 1]){//如果当前的分母比上
			st[x] = b;
			if(!flag || tmp[x] < ans[x])//找到了更优解
				for(int i = 1;i <= dep;i++)
					ans[i] = tmp[i];
		}
		flag = 1;//标记为已经找到了答案
		return;
	}
	int l = max((b + a - 1) / a,tmp[x - 1] + 1),r = min((dep - x + 1) * b / a,INF);//下一个分数分母的上下界
	if(flag && r >= ans[dep])
		r = ans[dep] - 1;
	for(int i = l;i <= r;i++){
		tmp[x] = i;
        //请自行模拟分数加减法
		int A = a * i - b,B = b * i;
		int gcd = GCD(A,B);
		dfs(A / gcd,B / gcd,x + 1);
	}
}
signed main(){
	a = read();b = read();
	c = GCD(a,b);
	a /= c;b /= c;
	tmp[0] = 1;
	for(dep = 1;dep <= N - 1;dep++){//枚举深度
		dfs(a,b,1);
		if(flag){找到了答案
			for(int i = 1;i <= dep;i++)
				printf("%lld ",ans[i]);
			return 0;
		}
	}
	return 0;
}

【2】初步剪枝

其实在上方的代码中已经给出了。

就是这一段：

int l = max((b + a - 1) / a/*(b + a - 1) / a = ceil(b / a)*/,tmp[x - 1] + 1);
int r = min((dep - x + 1) * b / a,INF);//下一个分数分母的上下界
if(flag && r >= ans[dep])
	r = ans[dep] - 1;

考虑上下界剪枝。

先考虑下界 $l$。这段代码中，$tm{p}_{x-1}$ 为上一个分数的分母，按照题目要求，下一个分数的分母要更大，所以下界 $l$ 至少为 $tm{p}_{x-1}+1$。又因为当前填的分数不能超过当前剩余的分数 $\frac{a}{b}$，所以有 $\frac{1}{tm{p}_x}\le \frac{a}{b}$，移项得 $tm{p}_x\ge \frac{b}{a}$。又因为 $tm{p}_x$ 为正整数，所以 $tm{p}_x\ge \lceil \frac{b}{a}\rceil$。综上，下界 $l=max(tm{p}_{x-1}+1,\lceil \frac{b}{a}\rceil )$。

上界 $r$ 呢？由题意，$tm{p}_{x+1},tm{p}_{x+2},\dots ,tm{p}_{dep}\ge tm{p}_x$，所以这些数的倒数小于 $\frac{1}{tm{p}_x}$，那么这些数倒数的和（包括 $\frac{1}{tm{p}_x}$）小于这些数的个数乘 $\frac{1}{tm{p}_x}$，即小于 $\frac{dep-x+1}{tm{p}_x}$。但这些数的和刚好等于当前的 $\frac{a}{b}$，所以 $\frac{dep-x+1}{tm{p}_x}>\frac{a}{b}$，移项得 $tm{p}_x<\frac{b(dep-x+1)}{a}$，又因为 $tm{p}_x$ 为正整数，所以 $tm{p}_x\le \lfloor \frac{b(dep-x+1)}{a}\rfloor$。题目中还说了最小的分数 $\ge \frac{1}{{10}^7}$，所以上界 $r$ 不会超过 ${10}^7$，综合有 $r=min(\lfloor \frac{b(dep-x+1)}{a}\rfloor ,{10}^7)$。

这里还有个小优化：当发现已经找到了解，那么如果当前的分数小于当前找到的最优解中最小的分数，那么接下来搜索到的解不会更优。

【3】神奇的剪枝

加了上下界剪枝依然无法通过本题，因为本题有一个很强的数据：

570 877

所以还要继续剪枝。

【3.1】神奇的剪枝 1：

我们发现，限制深度可以大大优化搜索，那么是否可以每次限制一些其它参数，进一步优化搜索呢？

题目还要求最小的分数越大越好，即最小的分数分母越小越好，考虑每次限制分母的最大值 $A$。

我的枚举是这样的：

限制深度为 $dep$ 时，先设 $A={10}^5$（这个可以自己调），因为规定了最小的分数分母小于 ${10}^7$，所以 $A\le {10}^7$。

如果每次将 $A$ 加 $1$，那么这个剪枝就没有任何意义了。因为这要求 $A$ 的初值很接近 ${10}^7$ 才能避免超时，而两个相差 $1$ 的 $A$ 对应的搜索树基本一致，相当于重复搜了很多次。

所以，我们每次将 $A$ 乘 $10$，有了这个限制，就可以对搜索进行优化。

因为题目首要满足的是分数个数尽可能小，所以枚举深度的循环在外层。

for(dep = 1;dep <= N - 1;dep++){//枚举深度
	ans[dep] = tmp[dep] = INF + 1;
	for(max_a = 100000;max_a <= INF;max_a *= 10){//即原文的A
		dfs(a,b,1);
		if(flag){//找到了答案
			for(int i = 1;i <= dep;i++)
				cout << ans[i] << ' ';
			return 0;
		}
	}
}

【3.2】神奇的剪枝 2：

这是本题最难想到但也是最妙的地方。

每次层数最深的地方搜索树的节点个数巨大，考虑对最深的两层优化。

根据题意，在 $x=dep-1$（倒数第二层）时，当前的 $\frac{a}{b}=\frac{1}{p}+\frac{1}{q}=\frac{p+q}{pq}(p<q,gcd(a,b)=1)$，因为 $a,b$ 互质，所以有：

$$\{\begin{array}{l}p+q=ak\\ pq=bk\end{array}$$

消元，得 $q=ak-p$，代入二式得 $p(ak-p)=bk$，整理得 $p^2-akp+bk=0$。

这是一个关于 $p$ 的一元二次方程，要求满足有正整数解，考虑枚举参数 $k$。

因为是一元二次方程，所以有 $\mathrm{\Delta }=a^2{k}^2-4bk\ge 0$，解这个关于 $k$ 的一元二次不等式得 $k\ge \frac{4b}{a^2}$ 或 $k\le 0$。又因为 $k$ 为正整数，所以 $k\ge \lceil \frac{4b}{a^2}\rceil$。

继续解方程，得 $p_1=\frac{ak+\sqrt{\mathrm{\Delta }}}{2},p_2=\frac{ak-\sqrt{\mathrm{\Delta }}}{2}$，要求为正整数，则 $\mathrm{\Delta }$ 为完全平方数且 $2|ak+\sqrt{\mathrm{\Delta }}$，所以当不满时直接跳过。否则，令 $tm{p}_{dep-1}=p_1,tm{p}_{dep}=p_2$ 即求得一组解。

$k$ 的下界 $l$ 确定了，接下来就确定上界 $r$。由上方方程的解，得 $ak=p_1+p_2$。因为我们限定了分母最大值 $A$，所以 $p_1<p_2\le A$，所以 $p_1+p_2<2A$，即 $ak<2A$，所以 $k<\frac{2A}{a}$。又因为 $k$ 为正整数，所以 $k\le \lfloor \frac{2A}{a}\rfloor$。

这还没完，将 $ak=p_1+p_2$ 代入原方程 $p^2-akp+bk=0$ 得 $bk=p_1{p}_2$，还是由 $p_1<p_2\le A$，得 $p_1{p}_2\le A(A-1)$，所以 $bk\le A(A-1)$ 得 $k\le \lfloor \frac{A(A-1)}{b}\rfloor$。

综上，$k$ 的上界 $r=min(\lfloor \frac{2A}{a}\rfloor ,\lfloor \frac{A(A-1)}{b}\rfloor )$。

if(x == dep - 1){
	const int l = ((b << 2) / (a * a)) + 1,r = min(((max_a << 1)) / a,max_a * (max_a - 1) / b);
	for(int i = l;i <= r;i++){
		int delta = a * a * i * i - ((b * i) << 2);
		int Sqrt = sqrt(delta);
		if(Sqrt * Sqrt != delta/*delta不为完全平方数*/ || ((a * i - Sqrt) & 1))
			continue;
		tmp[x] = (a * i - Sqrt) >> 1;tmp[x + 1] = (a * i + Sqrt) >> 1;
		if(!flag || tmp[dep] < ans[dep]){
			for(int i = 1;i <= dep;i++)
				ans[i] = tmp[i];
			flag = 1;
		}
	}
	return;
}

AC code

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N = 11;
const int INF = 1e7;
int dep;int max_a;
int tmp[N],ans[N];
int a,b,c;
bool flag;
int GCD(const int x,const int y){
	return y ? GCD(y, x % y) : x;
}
void dfs(const int a,const int b,const int x){
	if(x > dep)
		return;
	if(a == 1){
		tmp[x] = b;
		if(!flag || tmp[dep] < ans[dep])
			for(int i = 1;i <= dep;i++)
				ans[i] = tmp[i];
		flag = 1;
		return;
	}
	if(x == dep - 1){
		const int l = ((b << 2) / (a * a)) + 1,r = min(((max_a << 1)) / a,max_a * max_a / b);
		for(int i = l;i <= r;i++){
			int delta = a * a * i * i - ((b * i) << 2);
			int Sqrt = sqrt(delta);
			if(Sqrt * Sqrt != delta || ((a * i - Sqrt) & 1))
				continue;
			tmp[x] = (a * i - Sqrt) >> 1;tmp[x + 1] = (a * i + Sqrt) >> 1;
			if(!flag || tmp[dep] < ans[dep]){
				for(int i = 1;i <= dep;i++)
					ans[i] = tmp[i];
				flag = 1;
			}
		}
		return;
	}
	int l = max((b + a - 1) / a,tmp[x - 1] + 1);
	int r = min((((dep - x + 1) * b + a - 1) / a),max_a);
	if(flag && r >= ans[dep])
		r = ans[dep] - 1;
	for(int i = l;i <= r;i++){
		tmp[x] = i;
		const int A = a * i - b,B = b * i;
		const int gcd = GCD(A,B);
		dfs(A / gcd,B / gcd,x + 1);
	}
}
signed main(){
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	cin >> a >> b;
	c = GCD(a,b);
	a /= c;b /= c;
	tmp[0] = 1;
	for(dep = 1;dep <= N - 1;dep++){
		ans[dep] = tmp[dep] = INF + 1;
		for(max_a = 100000;max_a <= INF;max_a *= 10){
			dfs(a,b,1);
			if(flag){
				for(int i = 1;i <= dep;i++)
					cout << ans[i] << ' ';
				return 0;
			}
		}
	}
	return 0;
}