题解 AT_abc332_e 【Lucky bag】

· · 题解

思路

本题要将 N 个物品放到 D 个袋子中，每个物品有一个重量 W_i，放入袋子中后，每个袋子内物品总重 x_i，最后求这些袋子内物品总重数据 x_i 的方差最小值。

方差式子可以简化：

V=\dfrac{1}{D}\sum\limits_{i=1}^D(x_i-\overline x)^2,

其中

\begin{aligned}&\sum\limits_{i=1}^D(x_i-\overline x)^2\\&=\sum\limits_{i=1}^D(x_i^2-2x_i\overline x+\overline x^2)\\&=\left(\sum\limits_{i=1}^Dx_i^2\right)-2\overline x\left(\sum\limits_{i=1}^Dx_i\right)+\left(\sum\limits_{i=1}^D\overline x^2\right)\\&=\left(\sum\limits_{i=1}^Dx_i^2\right)-2\overline xD\overline x+D\overline x^2\\&=\left(\sum\limits_{i=1}^Dx_i^2\right)-D\overline x^2\\&=\left(\sum\limits_{i=1}^Dx_i^2\right)-\dfrac{1}{D}\left(\sum\limits_{i=1}^Dx_i\right)^2,\end{aligned}

所以代回原式中

V=\dfrac{D\left(\sum\limits_{i=1}^Dx_i^2\right)-\left(\sum\limits_{i=1}^Dx_i\right)^2}{D^2}.

根据质量守恒定律，我们知道 \sum\limits_{i=1}^Dx_i=\sum\limits_{i=1}^NW_i，所以上式分子的右边那项是常数项，关键要求出 \min\sum\limits_{i=1}^Dx_i^2。

求解方法

（随机化+贪心）我们充分发挥人类智慧，按照某种顺序依次将物品加入内部重量最小的袋子里，只要我们枚举的物品随机加入顺序足够多，所有方案的最小值就很可能是

\min\sum\limits_{i=1}^Dx_i^2.

使用 std::shuffle 枚举加入顺序即可：

#include<random>
#include<vector>
#include<iomanip>
#include<iostream>
#include<algorithm>
using namespace std;
using ll=long long;
int n,d;
long long x2(const vector<int>&w)
{
    vector<ll> ds(d);
    for(auto i:w)*min_element(ds.begin(),ds.end())+=i;
    //每次放到最小的袋子
    ll ans=0;for(auto i:ds)ans+=i*i;return ans;
    //求平方和
}
int main()
{
    cin>>n>>d;ll sum=0;
    vector<int> w(n);for(auto&i:w)cin>>i,sum+=i;
    mt19937 rnd(0xab8d5f6);
    ll minx2=9e18;
    for(int cnt=5e6;cnt--;)//卡一个 1s 左右的时间就好
    {
        shuffle(w.begin(),w.end(),rnd);//随机打乱顺序
        minx2=min(minx2,x2(w));//贪心求解，并取最小值
    }
    cout<<fixed<<setprecision(15)<<(minx2*d-sum*sum)*1.0l/d/d;
    //此处需要注意精度问题，尽量在做除法时再转化成long double
    return 0;
}

也可以通过状态压缩动态规划求解：

（状态压缩动态规划）设状态 \mathrm{dp}_{i,S} 表示包含了物品集合 S 的 i 个袋子内重量平方和的最小值，转移方程为。

\mathrm{dp}_{i,S}=\min\limits_{T\subseteq S}\mathrm{dp}_{i-1,S-T}+\mathrm{dp}_{1,T},\quad i\geqslant 2.

集合的表示用全集中某个元素是否存在在该集合（0、1）表示，因而一个元素对应一个 bit，枚举子集的位运算方法则是初始化 t=s，下一个集合表示是 t=(t-1)&s，枚举到 t==0 时可以结束。(t-1)&s 表示如果 t 不为 0（T 非空），则 t-1 相当于将 t 最小的 bit 去掉，加上所有位次比它低的 bit，由于又包含在 s 中，再将结果和 s 做与运算，从而排除掉在 s 之外的 bit；如果 t 为 0（T 为空），假设 t 是带符号补码实现的整数，那么 (t-1)&s 结果为 s。

程序对应的算法时间复杂度分析：

枚举所有集合表示 s，以及其子集表示 t 的过程中，考虑到一个物品有三种情况：包含在 S 和 T 中；包含在 S 但不包含在 T 中；不包含在 S 和 T 中。（由于 T\subseteq S，所以不能只在 T 中而不在 S 中）由于有 N 个物品，可知有 3^N 种不同的方案数。

再考虑上上袋子的维数，可知对应算法状态转移的时间复杂度为 \Theta(D3^N)。

#include<cstdio>
#include<algorithm>
long long dp[16][1<<15],w[15];
int main()
{
    int n,d;scanf("%d%d",&n,&d);long long sum=0;
    for(int i=0;i<n;i++)scanf("%lld",w+i),sum+=w[i];
    __builtin_memset(dp,0x3f,sizeof dp);
    for(int s=0;s<1<<n;s++)
    {
        dp[1][s]=0;
        for(int i=0;i<n;i++)if(s&1<<i)dp[1][s]+=w[i];
        dp[1][s]*=dp[1][s];
        //预处理 dp[1][s]，即只有一个袋子时 s 状态对应的总重平方
    }
    for(int i=2;i<=d;i++)for(int s=0;s<1<<n;s++)
    {
        int t=s;
        do{
            dp[i][s]=std::min(dp[i][s],dp[i-1][s^t]+dp[1][t]);
        }while((t=(t-1)&s)!=s);
        //for(int t=s;t;t=(t-1)&s)中 t 取遍所有非空子集
        //此处是枚举所有子集，既不是非空子集也不是真子集
        //当然，其实不必枚举所有子集
    }
    printf("%.15Lf",(dp[d][(1<<n)-1]*d-sum*sum)*1.0l/d/d);
    //此处需要注意精度问题，尽量在做除法时再转化成long double
    return 0;
}

由于枚举子集状态转移时间复杂度不低，我们有很大的把握认为，随机化卡不过不去，状压 dp 也不好受。

（优化第一维递推的状态压缩动态规划）还可以再优化一下上面的状态转移动态规划，我们发现袋子是一模一样的，而且不一定每次增加一个袋子作转移，观察上面给的 dp 转移方程：

\mathrm{dp}_{i,S}=\min\limits_{T\subseteq S}\mathrm{dp}_{i-1,S-T}+\mathrm{dp}_{1,T},\quad i\geqslant 2,

发现第一维不一定要从上一项转移而来，如果我们定义 \mathrm{dp}_{0,S}=\begin{cases}0,&S=\varnothing\\+\infty,&S\neq\varnothing\end{cases}，实际上有以下式子成立

\mathrm{dp}_{i,S}=\min\limits_{T\subseteq S}\mathrm{dp}_{i-j,S-T}+\mathrm{dp}_{j,T},\quad\forall 0\leqslant j\leqslant i,\forall S\subseteq U,

其中 U 是全集，忽略掉第二维，将转移方程写作：

\mathrm{dp}_i=\mathrm{dp}_{i-j}\oplus\mathrm{dp}_j,\forall 0\leqslant j\leqslant i,

此处 \oplus 不是异或运算，该运算就是第二维集合的转移运算：

\mathrm{dp}_{a+b,\varnothing}\\ \mathrm{dp}_{a+b,\{W_1\}}\\ \mathrm{dp}_{a+b,\{W_2\}}\\ \vdots\\ \mathrm{dp}_{a+b,U} \end{pmatrix}=\begin{pmatrix} \min\limits_{T\subseteq\varnothing}\mathrm{dp}_{a,\varnothing-T}+\mathrm{dp}_{b,T}\\ \min\limits_{T\subseteq\{W_1\}}\mathrm{dp}_{a,\{W_1\}-T}+\mathrm{dp}_{b,T}\\ \min\limits_{T\subseteq\{W_2\}}\mathrm{dp}_{a,\{W_2\}-T}+\mathrm{dp}_{b,T}\\ \vdots\\ \min\limits_{T\subseteq U}\mathrm{dp}_{a,U-T}+\mathrm{dp}_{b,T} \end{pmatrix}=\mathrm{dp}_a\oplus\mathrm{dp}_b.

容易发现这玩意满足结合律，于是我们还可以定义乘法：

\mathrm{dp}_a\otimes b=\underbrace{\mathrm{dp}_a\oplus\mathrm{dp}_a\oplus\cdots\oplus\mathrm{dp}_a}

括起来的部分共有 b 个 \mathrm{dp}_a 相加。

我们意识到这和倍增“快速幂”的思路不谋而合：

a^{19}=a^{10011_{(2)}}=a^{16}\times a^{2}\times a,

同样有

\mathrm{dp}_{19}=\mathrm{dp}_{10011_{(2)}}=\mathrm{dp}_{16}\oplus\mathrm{dp}_2\oplus\mathrm{dp}_1.

除此之外也可以将初始化 \mathrm{dp}_1 的代码写成对应时间复杂度 \Theta\left(2^N\right) 的算法，将上面的代数方法抽象成类/结构体后有了下面的代码，对应算法的时间复杂度就成了 \Theta\left(\log D\times 3^N\right)，即便 D 本身不大，但是整体效率依然提高了一些：

#include<cstdio>
#include<vector>
#include<algorithm>
using namespace std;
using ll=long long;
int n,d;vector<ll> w;
struct DP
{
    vector<ll> dp;
    void init_zero()//初始化为 dp_0，就是零元
    {
        dp=vector<ll>(1<<n,0x3f3f3f3f3f3f3f3f);
        dp[0]=0;
    }
    void init_one()//初始化为 dp_1，就是幺元
    {
        dp=vector<ll>(1<<n);
        for(int i=0;i<n;i++)dp[1<<i]=w[i];
        for(int s=0;s<1<<n;s++)dp[s]=dp[s^(s&-s)]+dp[s&-s];
        for(int s=0;s<1<<n;s++)dp[s]*=dp[s];
    }
    DP operator+(const DP&o)//新概念加法，也就是枚举子集转移方程
    {
        DP ans;ans.init_zero();
        for(int s=0;s<1<<n;s++)
        {
            int t=s;
            do{
                ans.dp[s]=min(ans.dp[s],dp[s^t]+o.dp[t]);
            }while((t=(t-1)&s)!=s);
        }
        return ans;
    }
};
DP doubling_mul(DP id,int b)//类似倍增“快速幂”/倍增“龟速乘”
{
    DP ans;ans.init_zero();
    for(;b;id=id+id,b>>=1)if(b&1)
    {
        ans=ans+id;
        if(b==1)break;//可以少算一次id+id
    }
    return ans;
}
int main()
{
    scanf("%d%d",&n,&d);ll sum=0;
    w=vector<ll>(n);for(auto&i:w)scanf("%d",&i),sum+=i;
    DP id;id.init_one();
    DP&&ans=doubling_mul(id,d);
    //id 是 dp_1，那么 dp 结果就包含在 id⊗d 中
    printf("%.15Lf",(ans.dp[(1<<n)-1]*d-sum*sum)*1.0l/d/d);
    return 0;
}