P1763 埃及分数（小理解，后续补充线性方程优化）

P1763 埃及分数

1、读题：

将一个真分数表示为一堆分子为 \(1\) 的分式相加，其中我们可以简单概括为

\[\frac{a}{b} = \frac{1}{x_1} + \frac{1}{x_2} + \dots \frac{1}{x_k} \]

对于题目要求是求解找到合适的最优解法之一。

满足的最优性质按照优先级排列如下：

• \(x\) 数组的数量最少，也就是 \(k\) 最小。
• 如果 \(k\) 相同的，分母最大值越小越好，为了方便比较，并且按照题目输出，默认 \(x_1 < x_2 \dots < x_k\) 。

2、思路：

2.1、最优性的分析

基础的思路实现是将整个答案拆解为一堆数字之和，也就是每次给定一个可选的数字，可以直接采用搜索的方式解决。

但是本题有一个优先级关系，输出的数组长度越小越好，即个数最少，

此时刚好满足迭代加深搜索的基本要求

• 答案个数存在单调优先关系，越少越好
• 找到合适的答案，就可以表示这个答案个数可以。

所以我们可以限定答案的个数（也就是搜索深度），采用迭代加深搜索，只要找到合适的答案，就可以停止继续加深搜索深度。

概括如下：

while (!dfs(2, a, b, 1))
{
    max_depth++;
}

2.2、搜索的实现和细节处理

搜索的过程主要分为三步：

• 截止条件，到达最大深度，停止搜索， \(u\) 表示上一个分母的选择， \(a, b\) 分别表示剩余的分子和分母， \(depth\) 表示当前的深度

  • bool dfs(ll u, ll a, ll b, int depth) 
          if (depth == max_depth) // 此时说明找到答案了
    

• 答案比较, 其中 \(res\) 数组存储真正的最优答案， \(tmp\) 表示的是当前这个搜索过程之中得到的一组解。

  • if (a == 1)  // 此时说明只剩下一个分数了，可以直接求解答案
    {
        // 如果分母大于了10^7，则与题目给出的数据范围不符，直接返回false
        if (b > 1e7) return false;
        tmp[depth - 1] = b;
        // 如果是第一次算出解，或者已经算出的解不够最优，则更新答案
        if (!res[0] || b < res[depth - 1]) 
        {
            for (int i = 0; i < max_depth; i++) res[i] = tmp[i];
        }
        return true;
    }
    return false;
    

  • 按照上述过程，需要关注的几个重点，只剩下一个分子为 \(1\) 的分数， 此时的分母不用枚举了，可以直接得到最终的分母。

  • 更新最优答案有两个选择，如果之前没有存储过答案，也就是 \(res[0] = 0\) ，直接将 \(tmp\) 先给到最优解存储，

    如果之前存储过，只需要比较最大的分母即可，因为此时长度是一样的。

• 搜索下一个的选项范围

  • 附带性的需要重点关注，因为此时范围把控范围不对，第一容易出现负数情况，第二容易超时。

  • 此时需要给出一定的不等式关系，帮助理解，如下所示, 假设当前枚举的分母是 \(i\) , \(a', b'\) 是剩余的分子和分母， 需要满足：

    • 条件 \(1\) : \(\frac{a'}{b'} > \frac{1}{i}\) ，因为最后一项单独判断，所以每次枚举的都不能超过剩余的分数，推导不等式可以得到

      \(i > \frac{a'}{b'}\) ，此时的答案是范围的下界。

    • 条件 \(2\) ：假设当前分子是 \(x_1 = i\), 后面的所有分母都是 \(i\)， 我们可以直接知道

      \[\frac{1}{x_1}+\frac{1}{x_2}+\dots+\frac{1}{x_{maxdep-dep + 1}} < \frac{1}{i}+\frac{1}{i}+\dots+\frac{1}{i} = (maxdep-dep+1)\times \frac{1}{i} \]

      如果这个数字按照最好情况都比 \(\frac{a'}{b'}\) 还小，那么这个数值一定不可能。

      也就是

      \[\frac{a'}{b'} < (maxdep-dep+1)\times \frac{1}{i} \\ i < (maxdep-dep+1)\times \frac{b'}{a'} \]

2.3、综合性代码的展示

搜索代码实现如下：

int max_depth = 1;// 迭代加深搜索 
// u表示本层应该从哪里开始枚举分母
bool dfs(ll u, ll a, ll b, int depth) 
{
    if (depth == max_depth) 
	{
        if (a == 1) 
		{
        	// 如果分母大于了10^7，则与题目给出的数据范围不符，直接返回false
            if (b > 1e7) return false;
            tmp[depth - 1] = b;
            // 如果是第一次算出解，或者已经算出的解不够最优，则更新答案
            if (!res[0] || b < res[depth - 1]) 
			{
                 for (int i = 0; i < max_depth; i++) res[i] = tmp[i];
            }
            return true;
        }
        return false;
    }
    bool flag = false;
    for (ll i = max(u, (b / a) + 1); i <= 1ll * b / a * (max_depth - depth + 1); ++i) 
	{
        ll nx = a * i - b, ny = b * i;
        ll g = gcd(nx, ny);
        nx /= g, ny /= g;
        tmp[depth - 1] = i;
        if (dfs(i + 1, nx, ny, depth + 1)) 
        // 另外还需要关注的是，当前深度还没有找到最优解之前，不能直接范围 
        {
        	flag = true;
		}
    }
    return flag;
}

主函数实现如下：

int main() 
{
    ll a,b;
    cin>>a>>b;
    ll Gcd = gcd(a,b);
    a /= Gcd;
    b /= Gcd;
    while (!dfs(2, a, b, 1))
    {
    	max_depth++;
	}
    for (int i = 0; i < max_depth; i++)
    {
    	cout<<res[i]<<' ';
	}
    return 0;
}

2.4、额外实现参考

之前求解的时候，尝试过bfs的写法过程，即使知道一定会超过内存限制，可以拿到虚假的 90pts

仅供参考

#include <bits/stdc++.h>
using namespace std;
#define ll long long  
int max_depth = 6;// 最大深度 
const int maxn = 1000;	
int tar_dep;// 目标深度 
ll a,b;
struct P{
	int num;// 当前分母的个数
	ll num_a;// 剩余的分子 
	ll num_b;// 剩余的分母 
	ll st;// 下一次可以选择的最小数字 
	vector<ll>res;// 答案数组 
};
ll ans[maxn]; 
ll gcd(ll a,ll b)// 求两个数字的最大公约数 
{
	return (b==0)?a:gcd(b,a%b);
}
void bfs()
{
	queue<P> Q;
	Q.push({1,a,b,2});// 初始状态 
	while (!Q.empty())
	{
		P u = Q.front();
		Q.pop();
		if (tar_dep!=0 && u.num>tar_dep)// 说明在这一层已经出现答案了，不需要继续往后一层找 
		{
		     break;
		}
		if(u.num_a==1)// 说明这就是最后一个分母 ,也就是这一层就是正确答案 
		{ 
	        if (u.num_b <= u.res[u.res.size()-1]) continue;// 比最后一个枚举的还小 
			u.res.push_back(u.num_b);
			if (tar_dep==0)// 说明第一次出现答案
			{
				tar_dep = u.num;//  更新目标深度 
				for (int i = 0; i < u.res.size(); i++) 
				{
                 	ans[i] = u.res[i];
				}
			} // 说明有更优的答案 
			if (u.num_b < ans[tar_dep-1]) 
			{
                 for (int i = 0; i < u.res.size(); i++) 
				 {
                 	ans[i] = u.res[i];
				 }
            }
		} 
		if (tar_dep!=0 && u.num==tar_dep)// 说明在这一层已经出现答案了，不需要继续更新 
		{
		     continue;	
		}
		for (ll i =  max(u.st, (u.num_b / u.num_a) + 1); i <= 1ll * u.num_b / u.num_a * (max_depth - u.num + 1); ++i) // 枚举优化
		{
	        ll nx = u.num_a * i - u.num_b;
			ll ny = u.num_b * i;
	        ll g = gcd(nx, ny);
	        nx /= g, ny /= g;
	        u.res.push_back(i);// 添加i作为答案
	        Q.push({u.num+1, nx, ny, i + 1,u.res});
	        u.res.pop_back();// 下一个枚举对象也要使用，记得清空数据
	    }
	}
	return ;
}
int main()
{
    
   
	cin>>a>>b;
    ll Gcd = gcd(a,b);
    a /= Gcd;
    b /= Gcd;
    if (a==1)// 特殊判断 
    {
    	cout<<b<<endl;
    	return 0;
	}
    bfs();
    for (int i=0;i<tar_dep;++i)
    {
    	cout<<ans[i]<<' ';
	}
	return 0;
}

3、实现优化：

3.1、基础优化 \(1\) ：直接将最后两个答案单独计算。

按照最新的数据来说，上述的时间复杂度过高，考虑优化。

上述实现过程之中，枚举的数值会越来越大，所以我们考虑最终求解的时候，

假设前面的答案搜索已经完成，

之前的截止条件是

\[\frac{a'}{b'} = \frac{1}{i} \]

此时只需要判断， \(a'\) 是否等于 \(1\) ，

那么当然也可以直接考虑如果只剩下两个分数的时候，是否可以直接求解，这样就可以少一重最大的递归枚举了

假设最终剩下的数字满足这两个关系

\[\frac{a'}{b'} = \frac{1}{x} + \frac{1}{y} \]

进行通分之后可以得到关系式

\[\frac{a'}{b'} = \frac{x + y}{xy} \]

此时又因为得到的 \(a'\) 和 \(b'\) 都是最简的分数的形式，也就是呈互质关系，可以得到表示关系如下：

\[\left\{\begin{matrix} a' \times p = x + y\\ b' \times p = x \times y \end{matrix}\right. \]

其中上述的 \(p\) 是公倍数的关系。

因为存在两个表达式，两个未知数，显然上式可以得到解，

将其中一个式子，代入到另外一个式子，可以得到

\[b' \times p = x \times (a' \times p - x) \]

所以上述过程就是一个关于 \(x\) 的一元二次方程组，可以利用数学知识点解决

这里我们直接补充， 对于一元二次方程组是 \(ax^2 + bx + c = 0\) 的形式

得到的解一共有两个，分别是

\[\left\{\begin{matrix} x_1 = \frac{-b - \sqrt{b^2 - 4ac}}{2a}\\ x_2 = \frac{-b + \sqrt{b^2 - 4ac}}{2a} \end{matrix}\right. \]

其中，我们通常将公共部分表示为一个符号，叫做 \(Delta\), 也就是 $\sqrt{b^2 - 4ac} = \Delta $

所以首先将我们得到的二次函数进行基本化

\[x^2 - a'\times p \times x + b' \times p = 0 \]

此时按照上述的公式，直接求解就可以得到两个解，分别是

\[\left\{\begin{matrix} x_1 = \frac{-a'\times p - \Delta}{2}\\ x_2 = \frac{-a'\times p + \Delta}{2} \end{matrix}\right. \]

此时的 \(x_1 = x, x_2 = y\) 。

其中首先需要判断 \(\Delta\) 是否是整数答案， \(\Delta = \sqrt{(a'\times p)^2 - 4 \times b' \times p}\)

另外还需要判断 \(x_1,x_2\) 是否也是整数答案

分步骤表示为， 下方的 \(i\) 代表的就是此时的 \(p\) 。

ll delta = a * a * i * i - ((b * i) << 2);
ll Sqrt = sqrt(delta);
if (Sqrt * Sqrt != delta || ((a * i - Sqrt) & 1))
{
    continue;
}
tmp[depth - 1] = (a * i - Sqrt) >> 1; // x_1
tmp[depth] = (a * i + Sqrt) >> 1; // x_2
if (!res[0] || tmp[depth] < res[depth])
{
    for (int i = 0; i < max_depth; i++)
    {
        res[i] = tmp[i];
    }
    return true;
}

3.2、基础优化 \(2\) ： 上下界优化

关键的一点是整个过程之中，还没有处理 \(p\) 的枚举范围，

此时我们可以考虑到整个答案过程要求的分母不能超过 \(10 ^ 7\) ，我们直接设置为

\(INF = 10 ^ 7\) ，

此时的话利用好之前的一些式子，可以得到最终的上下界边界

边界 \(1\) ：

\(\Delta = \sqrt{(a'\times p)^2 - 4 \times b' \times p} > 0\)

推导后可以得到

\(p > \frac{4\times b'}{{a'}^2}\) 。

边界 \(2\) ：

\[\left\{\begin{matrix} a' \times p = x + y = x_1 + x_2\\ b' \times p = x \times y = x_1 \times x_2 \end{matrix}\right. \]

我们可以知道

\[x_1 + x_2 < 2 \times INF \\ x_1 \times x_2 < INF \times INF \]

所以可以直接知道

\(p < min(\frac{2 \times INF}{a'}, \frac{INF \times INF}{b'})\) 。

最终的枚举范围可以确定，但是此时还没有剪枝通过。

可以考虑最终剩余三个数字的情况，但是此时的关系式只有两个，优化起来过于麻烦。

关键代码实现：

#define ll long long
const int maxn = 1e7;
ll INF = 1e7;
ll res[maxn], tmp[maxn];
ll gcd(ll a, ll b)
{
    return (b == 0) ? a : gcd(b, a % b);
}
int max_depth = 1; // 迭代加深搜索
// u表示本层应该从哪里开始枚举分母
bool dfs(ll u, ll a, ll b, int depth)
// 当前选择 分子 分母  深度
{
    if (depth == max_depth)
    {
        if (a == 1)
        {
            // 如果分母大于了10^7，则与题目给出的数据范围不符，直接返回false
            if (b > 1e7)
            {
                return false;
            }
            tmp[depth - 1] = b;
            // 如果是第一次算出解，或者已经算出的解不够最优，则更新答案
            if (!res[0] || b < res[depth - 1])
            {
                for (int i = 0; i < max_depth; i++)
                {
                    res[i] = tmp[i];
                }
            }
            return true;
        }
        return false;
    }
    if (depth == max_depth - 1) // 说明到达倒数第二层
    {
        ll l = ((b << 2ll) / a / a) + 1;
        ll r = min(((INF<< 1ll)) / a, (INF - 1)* (INF / b));
        for (ll i = l; i <= r; i++)
        {
            ll delta = a * a * i * i - ((b * i) << 2);
            ll Sqrt = sqrt(delta);
            /*
            delta不为完全平方数或者为0的答案需要去除 
            */
            if (Sqrt * Sqrt != delta || ((a * i - Sqrt) & 1))
            {
                continue;
            }
            tmp[depth - 1] = (a * i - Sqrt) >> 1;
            tmp[depth] = (a * i + Sqrt) >> 1;
            if (!res[0] || tmp[depth] < res[depth])
            {
                for (int i = 0; i < max_depth; i++)
                {
                    res[i] = tmp[i];
                }
                return true;
            }
        }
        return false;
    }
    bool flag = false;
    // 1/x <= b/a  x>= (b/a)+1    maxdep-dep+1 * (1/x) >= a/b
    // 枚举当前层的分数的分母
    for (ll i = max(u, (b / a) + 1); i <= 1ll * b / a * (max_depth - depth + 1); ++i)
    {
        ll nx = a * i - b, ny = b * i;
        ll g = gcd(nx, ny);
        nx /= g, ny /= g;
        tmp[depth - 1] = i;
        if (dfs(i + 1, nx, ny, depth + 1))
        {
            flag = true;
        }
    }
    return flag;
}

3.3、附加优化：

思考一个实现过程，如果修改 \(INF\) 的数值，可以在一个比较小的时候找到答案，一定比在更大的时候优秀。

所以提前将初始的 \(INF\) 设置为 \(1e6\)，如果可以找到答案，就直接退出，否则修改为 \(1e7\) 。

4、小结：

• 利用迭代加深进行搜索的方式是必学项，广搜虽然可以较快的处理，但是实际的空间浪费的较多，另外广搜的实现中，必须找到当前这一层的最优解才能退出，这对于初学者显然不是非常友好。
• 基本的搜索剪枝也是必学项，最优性的关系可以处理上界（范围的最大值），比上一个答案大可以处理下界（范围的最小值）。
• 至于最后单独处理两个答案，本质上是数学的方法，当然剩余三个数，四个数，可以采用求解线性方程组的秩提前处理，但是作为基本选项，难度考察上暂时没有到达这个要求，学会理解二次函数即可，如果还没有学过解二次函数，建议先理解透彻前面的剪枝优化，数学学到这个时候，理解这个会更加透彻，不然欲速则不达。