2012-A3-0004

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

题目要求四个点为圆心，确定每个圆的半径，使每个圆与其他圆无公共部分面积，且半径之和最大.

这题线性规划的模型非常明显.

每个圆的半径设为ri,

约束是

对任意i!=j, ri + rj <= dist(Pi,Pj) 

且有 ri >=0

目标是 max(sum(ri; i = 1..4))

可以用单纯形模版做此题

由于约束形式和目标函数简单，

也可以直接枚举单纯形的顶点求得最值。


{{{
#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
const int MAXS = 10;//sum of n and m
const double eps = 1e-7;
const double oo = 1e50;//infinite
#define clr(a) memset(a,0,sizeof(a))
inline bool equ(double a, double b) {
return (a < b) ? (a + eps >= b) : (a <= b + eps);
}
class LinearProgramming {
public:
double tA[MAXS+1][MAXS+1], tb[MAXS+1], tc[MAXS+1];

double A[MAXS+1][MAXS+1], b[MAXS+1], c[MAXS+1], v;
int N[MAXS+2], B[MAXS+2], n, m, cnt, failtype;
void init(int _n, int _m) {
n = _n; m = _m; cnt = 1; failtype = 0;
memset(c, 0, sizeof(double) * (n + m + 1));
clr(tA);clr(tb);clr(tc);clr(A),clr(b);
clr(N);clr(B);
}
void constraint(double _A[], double _b) {//_A * x <= b
memcpy(A[n + cnt] + 1, _A, sizeof(double) * n);
b[n + cnt++] = _b;
}
void coefficient(double _c[]) { memcpy(c + 1, _c, sizeof(double) * n); }
void pivot(int l, int e) {
int i, j;
tb[e] = b[l] / A[l][e];
tA[e][l] = 1 / A[l][e];
for(i = 1; i <= N[0]; i++)
if(N[i] != e)
tA[e][N[i]] = A[l][N[i]] / A[l][e];
for(i = 1; i <= B[0]; i++) {
tb[B[i]] = b[B[i]] - A[B[i]][e] * tb[e];
tA[B[i]][l] = -A[B[i]][e] * tA[e][l];
for(j = 1; j <= N[0]; j++)
if(N[j] != e)
tA[B[i]][N[j]] = A[B[i]][N[j]] - tA[e][N[j]] * A[B[i]][e];
}//for(i)
v += tb[e] * c[e];
tc[l] = -tA[e][l] * c[e];
for(i = 1; i <= N[0]; i++)
if(N[i] != e)
tc[N[i]] = c[N[i]] - tA[e][N[i]] * c[e];
for(i = 1; i <= N[0]; i++) if(N[i] == e) N[i] = l;
for(i = 1; i <= B[0]; i++) if(B[i] == l) B[i] = e;
for(i = 1; i <= B[0]; i++) {
for(j = 1; j <= N[0]; j++) A[B[i]][N[j]] = tA[B[i]][N[j]];
b[B[i]] = tb[B[i]];
}
for(i = 1; i <= N[0]; i++) c[N[i]] = tc[N[i]];
}
bool opt(){//false stands for unbounded
while (true) {
int l, e, tl, te;
double maxUp = -1;
for(int ie = 1; ie <= N[0]; ie++) {
if(c[te = N[ie]] <= eps) continue;
double d = oo;
tl = MAXS + 1;
for(int i = 1; i <= B[0]; i++)
if(A[B[i]][te] > eps) {
double temp = b[B[i]] / A[B[i]][te];
if(equ(d, oo) || temp < d || equ(temp, d) && B[i] < tl) {
d = temp; tl = B[i];
}
}
if(tl == MAXS + 1) return false;
if(d * c[te] > maxUp) {
maxUp = d * c[te];
l = tl;
e = te;
}
}//for(ie)
if(equ(maxUp, -1)) break;
pivot(l, e);
}
return true;
}
void delete0() {
int p, i;
for(p = 1; p <= B[0]; p++) if(B[p] == 0) break;
if(p <= B[0]) {
for(i = 1; i <= N[0]; i++)
if(A[0][N[i]] > eps || A[0][N[i]] < -eps) break;
pivot(0, N[i]);
}
for(p = 1; p <= N[0]; p++) if(N[p] == 0) break;
for(i = p; i < N[0]; i++) N[i] = N[i+1];
N[0]--;
}
bool initialize() {//false stands for infeasible
int i;
N[0] = n; B[0] = m; v = 0;
for(i = 1; i <= n; i++) N[i] = i;
for(i = 1; i <= m; i++)
B[i] = n + i;
int l = B[1];
for(i = 2; i <= m; i++) if(b[B[i]] < b[l]) l = B[i];
if(b[l] >= 0) return true;
double origC[MAXS+1];
memcpy(origC, c, sizeof(double) * (n + m + 1));
N[++N[0]] = 0;
for(i = 1; i <= B[0]; i++) A[B[i]][0] = -1;
memset(c, 0, sizeof(double) * (n + m + 1));
c[0] = -1; pivot(l, 0); opt();
if(v < -eps) return false;
delete0();
memcpy(c, origC, sizeof(double) * (n + m + 1));
bool inB[MAXS+1];
memset(inB, false, sizeof(bool) * (n + m + 1));
for(i = 1; i <= B[0]; i++) inB[B[i]] = true;
for(i = 1; i <= n + m; i++)
if(inB[i] && c[i] != 0) {
v += c[i] * b[i];
for(int j = 1; j <= N[0]; j++)
c[N[j]] -= A[i][N[j]] * c[i];
c[i] = 0;
}
return true;
}
bool simplex(bool calcEachAns = false) {
if(!initialize()) { failtype = 1; return false; }
if(!opt()) { failtype = 2; return false; }
if(!calcEachAns) return true;
bool inN[MAXS+1];
memset(inN, false, sizeof(bool) * (n + m + 1));
for(int i = 1; i <= N[0]; i++) inN[N[i]] = true;
for(int i = 1; i <= n; i++) if(inN[i]) b[i] = 0;
return true;
}
}mt;

int x[4],y[4];
int main(){
    while(scanf("%d%d",x,y)==2){
        mt.init(4,6);
        for(int i=1;i<4;++i)scanf("%d%d",x+i,y+i);
        int k;k=1;
        for(int i=0;i<4;++i)
            for(int j=i+1;j<4;++j){
                mt.A[4+k][i+1]=mt.A[4+k][j+1]=1;
                mt.b[4+k]=hypot((double)(x[i]-x[j]),(double)(y[i]-y[j]));
                k++;
            }
        mt.c[4]=mt.c[1]=mt.c[2]=mt.c[3]=1;
        mt.simplex(true);
        printf("%.8lf\n",mt.v);
    }
}
}}}

题目要求四个点为圆心，确定每个圆的半径，使每个圆与其他圆无公共部分面积，且半径之和最大.

这题线性规划的模型非常明显.

每个圆的半径设为ri,

约束是

对任意i!=j, ri + rj <= dist(Pi,Pj)

且有 ri >=0

目标是 max(sum(ri; i = 1..4))

可以用单纯形模版做此题

由于约束形式和目标函数简单，

也可以直接枚举单纯形的顶点求得最值。

#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
const int MAXS = 10;//sum of n and m
const double eps = 1e-7;
const double oo = 1e50;//infinite
#define clr(a) memset(a,0,sizeof(a))
inline bool equ(double a, double b) {
return (a < b) ? (a + eps >= b) : (a <= b + eps);
}
class LinearProgramming {
public:
double tA[MAXS+1][MAXS+1], tb[MAXS+1], tc[MAXS+1];
double A[MAXS+1][MAXS+1], b[MAXS+1], c[MAXS+1], v;
int N[MAXS+2], B[MAXS+2], n, m, cnt, failtype;
void init(int _n, int _m) {
n = _n; m = _m; cnt = 1; failtype = 0;
memset(c, 0, sizeof(double) * (n + m + 1));
clr(tA);clr(tb);clr(tc);clr(A),clr(b);
clr(N);clr(B);
}
void constraint(double _A[], double _b) {//_A * x <= b
memcpy(A[n + cnt] + 1, _A, sizeof(double) * n);
b[n + cnt++] = _b;
}
void coefficient(double _c[]) { memcpy(c + 1, _c, sizeof(double) * n); }
void pivot(int l, int e) {
int i, j;
tb[e] = b[l] / A[l][e];
tA[e][l] = 1 / A[l][e];
for(i = 1; i <= N[0]; i++)
if(N[i] != e)
tA[e][N[i]] = A[l][N[i]] / A[l][e];
for(i = 1; i <= B[0]; i++) {
tb[B[i]] = b[B[i]] - A[B[i]][e] * tb[e];
tA[B[i]][l] = -A[B[i]][e] * tA[e][l];
for(j = 1; j <= N[0]; j++)
if(N[j] != e)
tA[B[i]][N[j]] = A[B[i]][N[j]] - tA[e][N[j]] * A[B[i]][e];
}//for(i)
v += tb[e] * c[e];
tc[l] = -tA[e][l] * c[e];
for(i = 1; i <= N[0]; i++)
if(N[i] != e)
tc[N[i]] = c[N[i]] - tA[e][N[i]] * c[e];
for(i = 1; i <= N[0]; i++) if(N[i] == e) N[i] = l;
for(i = 1; i <= B[0]; i++) if(B[i] == l) B[i] = e;
for(i = 1; i <= B[0]; i++) {
for(j = 1; j <= N[0]; j++) A[B[i]][N[j]] = tA[B[i]][N[j]];
b[B[i]] = tb[B[i]];
}
for(i = 1; i <= N[0]; i++) c[N[i]] = tc[N[i]];
}
bool opt(){//false stands for unbounded
while (true) {
int l, e, tl, te;
double maxUp = -1;
for(int ie = 1; ie <= N[0]; ie++) {
if(c[te = N[ie]] <= eps) continue;
double d = oo;
tl = MAXS + 1;
for(int i = 1; i <= B[0]; i++)
if(A[B[i]][te] > eps) {
double temp = b[B[i]] / A[B[i]][te];
if(equ(d, oo) || temp < d || equ(temp, d) && B[i] < tl) {
d = temp; tl = B[i];
}
}
if(tl == MAXS + 1) return false;
if(d * c[te] > maxUp) {
maxUp = d * c[te];
l = tl;
e = te;
}
}//for(ie)
if(equ(maxUp, -1)) break;
pivot(l, e);
}
return true;
}
void delete0() {
int p, i;
for(p = 1; p <= B[0]; p++) if(B[p] == 0) break;
if(p <= B[0]) {
for(i = 1; i <= N[0]; i++)
if(A[0][N[i]] > eps || A[0][N[i]] < -eps) break;
pivot(0, N[i]);
}
for(p = 1; p <= N[0]; p++) if(N[p] == 0) break;
for(i = p; i < N[0]; i++) N[i] = N[i+1];
N[0]--;
}
bool initialize() {//false stands for infeasible
int i;
N[0] = n; B[0] = m; v = 0;
for(i = 1; i <= n; i++) N[i] = i;
for(i = 1; i <= m; i++)
B[i] = n + i;
int l = B[1];
for(i = 2; i <= m; i++) if(b[B[i]] < b[l]) l = B[i];
if(b[l] >= 0) return true;
double origC[MAXS+1];
memcpy(origC, c, sizeof(double) * (n + m + 1));
N[++N[0]] = 0;
for(i = 1; i <= B[0]; i++) A[B[i]][0] = -1;
memset(c, 0, sizeof(double) * (n + m + 1));
c[0] = -1; pivot(l, 0); opt();
if(v < -eps) return false;
delete0();
memcpy(c, origC, sizeof(double) * (n + m + 1));
bool inB[MAXS+1];
memset(inB, false, sizeof(bool) * (n + m + 1));
for(i = 1; i <= B[0]; i++) inB[B[i]] = true;
for(i = 1; i <= n + m; i++)
if(inB[i] && c[i] != 0) {
v += c[i] * b[i];
for(int j = 1; j <= N[0]; j++)
c[N[j]] -= A[i][N[j]] * c[i];
c[i] = 0;
}
return true;
}
bool simplex(bool calcEachAns = false) {
if(!initialize()) { failtype = 1; return false; }
if(!opt()) { failtype = 2; return false; }
if(!calcEachAns) return true;
bool inN[MAXS+1];
memset(inN, false, sizeof(bool) * (n + m + 1));
for(int i = 1; i <= N[0]; i++) inN[N[i]] = true;
for(int i = 1; i <= n; i++) if(inN[i]) b[i] = 0;
return true;
}
}mt;
int x[4],y[4];
int main(){
    while(scanf("%d%d",x,y)==2){
        mt.init(4,6);
        for(int i=1;i<4;++i)scanf("%d%d",x+i,y+i);
        int k;k=1;
        for(int i=0;i<4;++i)
            for(int j=i+1;j<4;++j){
                mt.A[4+k][i+1]=mt.A[4+k][j+1]=1;
                mt.b[4+k]=hypot((double)(x[i]-x[j]),(double)(y[i]-y[j]));
                k++;
            }
        mt.c[4]=mt.c[1]=mt.c[2]=mt.c[3]=1;
        mt.simplex(true);
        printf("%.8lf\n",mt.v);
    }
}