01分数规划

问题：给定一个正整数m和两个正数数组a[n]和b[n]。对于数组x[n]，x[n]中有m个i使得x[i] = 0，其余x[i] = 1，并且令R = segma(a[i]*x[i]) / segma(b[i]*x[i])。求一个满足条件的x[n]，使得R的值最大，输出R的最大值。

分析：构造函数F(L) = segma(a[i]*x[i]) - L * segma(b[i]*x[i]) = segma((a[i] - L * b[i]) * x[i])。令c[i] = a[i] - L * b[i]。

　　　观察F(L)，发现以下性质：1、若L固定，则可以通过一次对c[i]的排序，找到使得F(L)取最值的数组x[n]；2、若x[n]固定，则F(L)随L的增大而减小；3、在2的基础上发现，即使x[n]不固定，F(L)也随着L的增大而减小；4、对于L0，若有F(L0) > 0，则由代数变形可知segma(a[i]*x[i]) / segma(b[i]*x[i]) > L0，所以R > L0。同理，若F(L0) < 0，则R < L0，若F(L0) = 0，则R = L0。

　　　则由性质3，4可知，该问题能用二分法得到解决。

　　　-------------------------------------------------------------二分法-------------------------------------------------------

　　　设置下界l和上届r，mid = (l+r)/2。另c[i] = a[i] - mid*b[i]，选出最大的(n-m)个c[i]之和，若和大于等于0，则l = mid，否则r = mid，直到满足精度为止。

　　　-----------------------------------------------------Dinkelbach算法分析------------------------------------------------

　　　但是，二分法解决此问题速度比较慢，还有一种难懂但是速度更快的算法Dinkelbach。

　　　对F(L)，设F(R) = 0。首先，随机选择一个数k1有两种情况：1、若F(k1) > 0，则k1 < R，记使得F(k1)取最大值的数组x[n]为x1[n]。令k2 = segma(a[i]*x1[i]) / segma(b[i]*x1[i])，由题设易知，k2 <= R，而F(k1) > 0 = F(k2)，故k2 > k1，所以有k1 < k2 <= R，故令k1 = k2可更接近R，且下一次仍然有F(k1) > 0；2、F(k1) < 0，则k1 > R，记使得F(k1)取得最大值的数组x[n]为x2[n]。令k2 = segma(a[i]*x2[i]) / segma(b[i]*x2[i])，由题设知k2 <= R，则令k1 = k2，从下依次开始，一直有F(k1) > 0，由第一种情况组成循环。故可以不断通过k1 = k2的方式使得k1 更接近R。

　　　-------------------------------------------------------Dinkelbach算法--------------------------------------------------

　　　随机取一个数k，令c[i] = a[i] - k * b[i]并对c[i]进行排序，选出最大的m个，令k' = segma(a'[]) / segma(b'[])，(其中a'[]和b'[]事排序后选出来的最大的m个)，若k'与k的差距在精度局限内则输出，否则k = k'，循环上诉步骤。

参考资料：http://www.cnblogs.com/perseawe/archive/2012/05/03/01fsgh.html

　　　　　http://blog.sina.com.cn/s/blog_6383bcba0100xf4z.html

　　　　　http://hi.baidu.com/misshanli/item/ad81b23110f7ccb3134b14c3

例题：

1、POJ 2976 Dropping tests

题意：给定一个正整数m和两个正数数组a[n]和b[n]。对于数组x[n]，x[n]中有m个i使得x[i] = 1，其余x[i] = 0，并且令R = segma(a[i]*x[i]) / segma(b[i]*x[i])。求一个满足条件的x[n]，使得R的值最大，输出R的最大值。(注意看，与上面的题目有所不同)

　　　1 <= n <= 1000，0 <= m < n，0 <= ai <= bi <= 10^9。

解法：不多说，完全就是上面那样。直接看代码。

tag:math, 01分数规划

二分法：

/*
 * Author:  Plumrain
 * Created Time:  2013-10-13 10:18
 * File Name: math-POJ-2976.cpp
 * 二分法
 */
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>

using namespace std;

#define CLR(x) memset(x, 0.0, sizeof(x))

const double eps = 1e-4;
const int N = 1000;

int n, m;
int a[N+5], b[N+5];
double c[N+5];

bool cmp(double x, double y)
{
    return x > y;
}

double gao()
{
    double l = 0.0, r = 1.0;
    while (l + eps < r){
        double mid = (l+r) / 2.0;
        CLR (c);
        for (int i = 0; i < n; ++ i)
            c[i] = a[i] - mid * b[i];
        sort(c, c+n, cmp);
        double tmp = 0.0;
        for (int i = 0; i < n-m; ++ i)
            tmp += c[i];
        if (tmp >= eps) l = mid;
        else r = mid;
    }
    return l;
}

int main()
{
    while (scanf ("%d%d", &n, &m) != EOF && n){    
        for (int i = 0; i < n; ++ i)
            scanf ("%d", &a[i]);
        for (int i = 0; i < n; ++ i)
            scanf ("%d", &b[i]);
        
        printf ("%d
", (int)(gao()*100.0 + 0.5));
    }
    return 0;
}

Dinkelbach算法：

/*
 * Author:  Plumrain
 * Created Time:  2013-10-13 17:29
 * File Name: math-POJ-2976-2.cpp
 * Dinkelbach算法
 */
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<utility>

using namespace std;

#define CLR(x) memset(x, 0.0, sizeof(x))
#define PB push_back

const double eps = 1e-4;
const int N = 1000;

int n, m;
int a[N+5], b[N+5];
vector<pair<double, int> > cur;

bool cmp(pair<double, int> x, pair<double, int> y)
{
    return x.first > y.first;
}

double gao()
{
    double p = 0.5, q = 0.0;
    while (1){
        cur.clear();
        for (int i = 0; i < n; ++ i){
            double c = a[i] - p * b[i];
            cur.PB (make_pair(c, i));
        }
        sort (cur.begin(), cur.end(), cmp);
        double tmp1 = 0.0, tmp2 = 0.0;
        for (int i = 0; i < n-m; ++ i){
            tmp1 += a[cur[i].second];
            tmp2 += b[cur[i].second];
        }
        q = tmp1 / tmp2;
        if (fabs(p-q) < eps) break;
        p = q;
    }
    return q;
}

int main()
{
    while (scanf ("%d%d", &n, &m) != EOF && n){
        for (int i = 0; i < n; ++ i)
            scanf ("%d", &a[i]);
        for (int i = 0; i < n; ++ i)
            scanf ("%d", &b[i]);

        printf ("%d
", (int)(gao()*100.0+0.5));
    }
    return 0;
}