题解：P1763 埃及分数
WangTianJiao
·
·
题解
题面
update 2025/8/9：优化、修改了代码使之能通过 hack 数据，同步修改了分析。去掉了非正解代码。
分析题干
对于分数 \dfrac{a}{b}，要求把它拆分成 n 个分子为 1 的分数 \dfrac{1}{i_1} + \dfrac{1}{i_2} + \cdots + \dfrac{1}{i_n}，满足：
-
-
- 满足第一条时，i_n 尽量小，即 \dfrac{1}{i_n} 尽量大
思路
显然不能采用贪心策略而在一开始尽量使用大分数，否则可能使剩下的分数极小而又不能简单分解，不得不分解成很多分数，使最后解不是最优解。考虑搜索算法来找出解。
考虑解答树的结构。
显然每层的数目极大，单纯的广搜空间和时间都不能接受（尤其空间）。如果使用深搜，那么可能会陷入深度极大但对答案完全没有贡献的枝杈，从而时间完全不能接受。因此，我们需要更优秀的搜索方法。
迭代加深搜索（IDDFS）
简要介绍一下迭代加深搜索（迭代深化深度优先搜索）：
在这个搜索策略中，一个具有深度限制的深度优先搜索算法会不断重复地运行，并且同时放宽对于搜索深度的限制，直到找到目标状态。——百度百科
也就是说，我们进行正常的深度优先搜索，但是限制搜索的深度，如果搜索时深度超过限制，无论如何立即回溯，继续搜索其它状态，从而遍历解答树上深度不超过限制的所有状态。如果找不到解，那就扩大限制，再次搜索；一旦找到结果，立即停止搜索、停止继续扩大限制。
$A$：不会，因为解答树的每层状态数是成指数级增长的，前几层的状态数目相对于最后的一两层来说九牛一毛，即使反复搜索也成本不高。也就是说，即使假设我们提前知道解的深度，直接进行一次该深度的搜索，也比迭代加深搜索快不了多少。
显然，这种思路可以避免广搜的空间紧张，也可以避免深搜陷入无底洞。当然，搜索出的解一定层数最少。
迭代加深分数个数（即层数），每层不断枚举 $i$，从 $\dfrac{a}{b}$ 中减去 $\dfrac {1}{i}$，并向下传递：
- 剩下需要搜索的分数
- 上一个分母
- 当前层数
- 最大层数
以此进行搜索。
此时很容易发现一个优化：搜索到最后一层时不必再枚举，而是直接判定当前剩下的分数是否符合条件：分子为 $1$，分母大于上个分母。
发现严重超时，连 $\dfrac{7}{10}$ 这种小数字都超时，需要剪枝。
## 剪枝 #1
思考可行性剪枝：搜索中各个分数值严格递减，而每一层最多剩下几个分数是固定的。显然，如果剩下的分数全部顶满最大值，仍小于剩下的分数，可以直接剪枝。
在搜索时写一句剪枝：
```cpp
if ((maxdep - step + 1) * (1 / double(i)) < double(a) / double(b))//判断剩下的层数能否容纳剩下的分数
break;
```
好些了。
按这样思路写好提交，分数介于 $40$ 到 $100$ 之间（实际分数我自己两次写好提交都不一样），但是 hack 数据超时。（搜索的时间复杂度是 $O(玄学)$，两次可能差着一点常数）
## 剪枝 #2
事实上，最大分母的量级很大，使得每一层搜索工作量极大。在每次搜索时，由一次固定最大分母搜索，改为多次层数相同、但最大分母差异较大的搜索。如果最大分母每次 $+1$ 显然功效不大，可以由 $10^3$ 每次 $\times 10$ 到 $10^7$。
```cpp
ll maxdepth = 1;
for (; maxdepth <= 1e5; ++maxdepth)//迭代加深最大深度
{
ll maxb = 1000;
while (maxb <= 10000000)//迭代加深最大分母，最多5次，不会浪费大量时间
{
memset(vans, 0, sizeof(vans));
memset(vnow, 0, sizeof(vnow));
nowp = ansp = 0;
dfs(1, a, b, 0, maxdepth, maxb);
if (ansp)//找到答案立刻输出、退出
{
for (ll i = 1; i <= ansp; ++i)
{
printf("%lld ", vans[i]);
}
return 0;
}
maxb *= 10;
}
}
```
但是即使这样仍然不能通过 hack 数据。它层数极多，前面的剪枝不能通过，因此还需要剪枝。
## 剪枝 #3（最妙的剪枝）
考察最后两层搜索，它们占了主要的时间。
我们想求 $\dfrac {a}{b} = \dfrac{1}{x} + \dfrac{1}{y}$ 的整数解，通分得 $\dfrac {a}{b} = \dfrac{x+y}{xy} $，由于 $\gcd(a,b)=1$，所以存在正整数 $k$，使得：
$$
\begin {cases}
ka=x+y\\
kb=xy
\end {cases}
$$
根据韦达定理，$x,y$ 是下面这个方程的两个解：
$$
q^2-kaq+kb=0
$$
此时 $\Delta = k^2a^2-4kb$，为了有两个不等的实根，须有 $\Delta > 0$，化简得 $k>\dfrac{4b}{a^2}$。
简单枚举 $k$，计算这个方程的两个根 $x,y$，即可在一次循环中求得剩余两层的解而无需再搜索。
在枚举 $k$ 时应有范围限制：已经求得 $k>\dfrac{4b}{a^2}$，还需求得上界。
设最大分母为 $maxb$，有：
$$
x<y\\
ka=x+y<2y<2\times maxb\\
k < \dfrac{2maxb}{a}\\
kb=xy<maxb(maxb-1)\\
k<\dfrac{maxb(maxb-1)}{b}\\
$$
即
$$
k \in \left( \dfrac{4b}{a^2},\min \left( \dfrac{2maxb}{a}, \dfrac{maxb(maxb-1)}{b}\right) \right)
$$
还需验证 $\Delta$ 是完全平方数。
计算：
$$
x=\dfrac{ka-\sqrt{\Delta}}{2}\\
y=\dfrac{ka+\sqrt{\Delta}}{2}\\
$$
此处 $ka$ 和 $k^2a^2$ 同奇偶，$-4kb$ 不影响奇偶性，若 $\Delta$ 是完全平方数时 $ka$ 也和 $\sqrt{\Delta}$ 同奇同偶，从而分子是偶数，保证了解是整数。
此时 $x,y$ （当然不能大于 $maxb$ ）即是解，直接装入。
```cpp
for (ll k = 4 * b / (a * a); k*a<=2*maxb; ++k){
ll delta = (k * k * a * a) - (4 * k * b), x, y;//计算判别式
if (delta <= 0) continue;
ll s = sqrt(delta);
if (s * s != delta || delta <= 0)//判别式必须是完全平方数，保证x,y均为整数
continue;
y = ((k * a) + s) >> 1;
x = ((k * a) - s) >> 1;//计算x,y
if (x > maxb || y > maxb || x < 0 || y < 0) continue;
if ((ansp == 0 || y < vans[maxdep])){
if (x == vnow[step - 1]) continue;
for (ll i = 1; i < step; ++i)
vans[i] = vnow[i];
vans[step] = x;
vans[step+1] = y;
ansp = maxdep;
break;
}
}
```
此时这段代码能搜能剪，性能十分优越，因而此时可以用 $0.2s$ 左右的时间通过 hack 数据。
完整代码：
```cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll gcd(ll a, ll b){
return b ? gcd(b, a % b) : a;
}
const double eps = 1e-9;
ll vans[105], vnow[105];
ll nowp, ansp, maxdep = 1, maxb;
void dfs(ll step, ll a, ll b, ll lastb){
if (!a)//直接跳过
return;
ll g = gcd(a, b);
a /= g, b /= g;
if (step == maxdep - 1){
/*
解方程部分：
对于最后两层a / b = 1 / x + 1 / y:
显然a / b = (x+y) / xy;且a/b最简
那么x+y = ka, xy = kb
依据韦达定理构造二次方程：q^2 - kaq + kb = 0;
如果要方程有不等的两个根:delta = k*k*a*a - 4kb;
那么解得k>4b/a^2;为下界
令maxb为分母最大值
显然x + y < maxb
k < 2 * maxb / a;
x < y < maxb
xy < maxb(maxb-1)
k < maxb(maxb-1) / b;
min(2*maxb/a, maxb*(maxb-1)/b)即为上界
*/
//另：不开long long 见祖宗
for (ll k = 4 * b / (a * a); k*a<=2*maxb; ++k){
ll delta = (k * k * a * a) - (4 * k * b), x, y;//计算判别式
if (delta <= 0) continue;
ll s = sqrt(delta);
if (s * s != delta || delta <= 0)//判别式必须是完全平方数，保证x,y均为整数
continue;
y = ((k * a) + s) >> 1;
x = ((k * a) - s) >> 1;//计算x,y
if (x > maxb || y > maxb || x < 0 || y < 0) continue;
if ((ansp == 0 || y < vans[maxdep])){
if (x == vnow[step - 1]) continue;
for (ll i = 1; i < step; ++i)
vans[i] = vnow[i];
vans[step] = x;
vans[step+1] = y;
ansp = maxdep;
break;
}
}
return;
}
for (ll i = lastb + 1; i < maxb; i++){
if (a * i < b)//a/b小于1/i可直接跳过
continue;
//若剩下层数全部填满与当前分数相同仍不足当前分数，return true;
if ((maxdep - step + 1) * (1 / double(i)) < double(a) / double(b))//判断剩下的层数能否容纳剩下的分数
break;
ll newa = a * i - b;
ll newb = i * b;
vnow[step] = i;//装入答案
dfs(step + 1, newa, newb, i);
vnow[step] = 0;
}
}
int main(){
ll a, b;
cin >> a >> b;
//a = 73, b = 313;
for (maxdep = 1; maxdep <= 15; ++maxdep){//迭代加深最大深度
maxb = 1000;
while (maxb <= 10000000){//迭代加深最大分母，最多5次，不会浪费大量时间
memset(vans, 0, sizeof(vans));
memset(vnow, 0, sizeof(vnow));
nowp = ansp = 0;
dfs(1, a, b, 0);
if (ansp){//找到答案立刻输出、退出
for (ll i = 1; i <= ansp; ++i)
printf("%lld ", vans[i]);
return 0;
}
maxb *= 10;
}
}
return 0;
}
```