P1763 题解

王熙文 · 2023-07-16 11:42:44 · 题解

前言

偶然做到这题，看到讨论区有一个 hack 将所有的题解都叉掉了，因此就想用我的代码通过这个 hack。经过了两天的优化，代码终于过了，因此来发一篇题解。

upd 2023.7.17: 修改了一些细节，重新交了一下题解。

思路

首先发现题目要求单位分数的个数最少，且直接 dfs 可能会进入一个错解出不来，bfs 会空间超限，所以使用迭代加深搜索，从小到大枚举选的分数个数 cnt，并 dfs。

注意到题目也要求最大的分母最小，且如果限制了分母的最大值有利于剪枝，所以也将最大的分母迭代加深，设限制为 ax，ax 从 10^3 到 10^7 每次乘 10。

可以参考下面这个代码理解一下。代码中 ch 表示当前选的数，ans 表示答案选的数。

void dfs(int now,int a,int b)
{
    if(now==cnt+1)
    {
        if(a!=0) return;
        if(ans[1]==0 || ch[cnt]<ans[cnt])
        {
            for(int i=1; i<=cnt; ++i) ans[i]=ch[i];
        }
        return;
    }
    for(int i=ch[now-1]+1; i<=ax; ++i)
    {
        int gta=a*i-b,gtb=b*i,g=gcd(gta,gtb);
        gta/=g,gtb/=g;
        ch[now]=i,dfs(now+1,gta,gtb);
    }
}

这样并不能通过，考虑优化代码。以下是我加的几个优化：

• 设接下来选的分数为 \dfrac{1}{i}，则 \dfrac{a}{b}\ge\dfrac{1}{i}，所以 i \ge \dfrac{b}{a}。

• 因为后面选的分数一定小于当前，所以必须满足 (cnt-now+1)\cdot \dfrac{1}{i} \ge \dfrac{a}{b} 才有可能满足。因此 i \le \dfrac{(cnt-now+1)\cdot b}{a}。

• 当 now \le cnt-1 时，需要满足下一次 i 存在（下一次 i 的限制区间左端点小于等于右端点）。因为 i 的限制有关 a,b 的都是两个数的比例（\dfrac{a}{b} 或 \dfrac{b}{a}），所以在考虑下一次的时候 a,b 不需要约分。设下一次的 a,b 为 nxta,nxtb。提出来 \dfrac{nxtb}{nxta} \le i 和 i \le ax 这两个条件，得 \dfrac{nxtb}{nxta} \le ax，即 \dfrac{bi}{ai-b} \le ax。所以 i \ge \dfrac{axb}{axa-b}。

• 用类似的方法推 now \le cnt-2（\dfrac{axbi}{ax(ai-b)-bi} \le ax）可以得到 i \ge \dfrac{axb}{axa-2b}。因此猜测对于任意的 now 都有 i \ge \dfrac{axb}{axa-(cnt-now)b}。使用类似的方法可以证明。

• 假设 cnt=3，那么 \dfrac{a}{b}=\dfrac{1}{x}+\dfrac{1}{y}+\dfrac{1}{z}=\dfrac{xy+xz+yz}{xyz}。又因为 x,y,z \le ax，所以若 a>cnt\cdot ax^2 或 b > ax^3，则无解。类似的，若在过程中 a > (cnt-now+1)\cdot ax^{cnt-now} 或 b > ax^{cnt-now+1} 则无解。

• 当 now=cnt 时，最后一个分数一定为 \dfrac{a}{b}，所以直接检查是否合法即可，不需要递归到下一层。

• 当 now=cnt-1 时，还剩下两个分数。设 \dfrac{a}{b}=\dfrac{1}{x}+\dfrac{1}{y}=\dfrac{x+y}{xy}。所以 a,b 同时乘上某个数（设其为 z）后分别等于 x+y,xy。因为 x+y \le 2ax，所以 z \le \dfrac{2ax}{a}，并不会太大。因此当 z 小于等于某个阈值时，可以直接枚举 z，并解方程组 \begin{cases}az=x+y\\bz=xy\end{cases} 即可得到最后两个分数，检查是否合法即可。

• 在上面优化的过程中，方程有不相等的两解当且仅当 \Delta=(az)^2-4bz > 0 即 z >\dfrac{4b}{a^2}，因此枚举的时候下界为 \lfloor \dfrac{4b}{a^2}\rfloor+1。

使用这些优化后，便可以快速通过讨论区的 hack。

代码

代码中使用了 __int128，因为不确定是否会爆 long long。

#include<bits/stdc++.h>
using namespace std;
#define int __int128
int ceil(int a,int b) { return (a+b-1)/b; }
int gcd(int a,int b) { return b==0?a:gcd(b,a%b); }
int a,b;
int cnt,ax;
int ch[1000010],ans[1000010];
int ppow[1000010];
int sum=0;
void dfs(int now,int a,int b)
{
    if(now==cnt-1 && min(ax*2/a,ax*ax/b)<=1000)
    {
        int tmp=min(ax*2/a,ax*ax/b);
        for(int i=(4*b/a/a)+1; i<=tmp; ++i)
        {
            int delta=a*i*a*i-4*b*i;
            int sqr=sqrtl((long double)delta)+2;
            while(sqr*sqr>delta) --sqr;
            if(sqr*sqr!=delta || (a*i+sqr)%2!=0) continue;
            ch[cnt-1]=(a*i-sqr)/2,ch[cnt]=(a*i+sqr)/2;
            if(ch[cnt-1]>ch[cnt-2] && ch[cnt]>ch[cnt-1] && ch[cnt]<ax)
            {
                ax=ch[cnt];
                for(int i=1; i<=cnt; ++i) ans[i]=ch[i];
                for(int i=1; i<=cnt; ++i) ppow[i]=min(ppow[i-1]*ax,(int)1e29);
            }
        }
        return;
    }
    if(now==cnt)
    {
        if(b%a!=0 || b/a<=ch[cnt-1]) return;
        ch[cnt]=b/a;
        if(ch[cnt]<ax)
        {
            ax=ch[cnt];
            for(int i=1; i<=cnt; ++i) ans[i]=ch[i];
            for(int i=1; i<=cnt; ++i) ppow[i]=min(ppow[i-1]*ax,(int)1e29);
        }
        return;
    }
    if(a>(cnt-now+1)*ppow[cnt-now] || b>ppow[cnt-now+1]) return;
    int l=max(ch[now-1]+1,max(ceil(b,a),ceil(b*ax,a*ax-(cnt-now)*b))),r=min(ax,(cnt-now+1)*b/a);
    if(now==cnt-1)
    {
        for(int i=l; i<=r; ++i) ch[now]=i,dfs(now+1,a*i-b,b*i);
    }
    else
    {
        for(int i=l; i<=r; ++i)
        {
            ch[now]=i;
            int g=gcd(a*i-b,b*i);
            dfs(now+1,(a*i-b)/g,b*i/g);
        }
    }
}
signed main()
{
    ppow[0]=1;
    long long re; cin>>re; a=re; cin>>re; b=re;
    int g=gcd(a,b); a/=g,b/=g;
    if(a==1) return cout<<(long long)b,0;
    cnt=2;
    while(1)
    {
        ans[cnt]=1e9;
        ax=1e3;
        while(ax<=1e7)
        {
            for(int i=1; i<=cnt; ++i) ppow[i]=min(ppow[i-1]*ax,(int)1e29);
            dfs(1,a,b);
            if(ans[1]!=0 && ans[1]!=1e9)
            {
                for(int i=1; i<=cnt; ++i) cout<<(long long)ans[i]<<' ';
                return 0;
            }
            ax*=10;
        }
        ++cnt;
    }
    return 0;
}