| ||||||||||
| Online Judge | Problem Set | Authors | Online Contests | User | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Web Board Home Page F.A.Qs Statistical Charts | Current Contest Past Contests Scheduled Contests Award Contest | |||||||||
我来可耻地贴一下代码吧:一个牛顿迭代的,一个模拟退火的……再次感叹一下模拟退火好简单!!In Reply To:其实我觉得这个算法还是叫“随机化变步长贪心”比较好,因为没有写出来“以一定的概率接受较差的解”=。= Posted by:Ruby931031 at 2012-10-02 14:03:14 //题目大意:已知n个点坐标,求其费马点,输出费马点到各个点距离之和,精确到个位,四舍五入
//算法:1.模拟退火;
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stdlib.h>
#define RANGE 10000
#define MAXN 110
#define eps 1e-3
const double pi=acos(-1.0);
struct Point{
double x,y;
}p[MAXN];
int n;
int SA(int cnt=3){ //cnt为执行模拟退火的次数,默认为3
double step,tmp,ans,theta;
double delta=pi/4,ret=1e40,rate=0.83; //降温速率
int i;
Point pt,cur;
srand(time(NULL));
while (cnt--) {
step = rate*RANGE; //初始温度
cur.x = (rand()%RANGE+1.0)/RANGE;
cur.y = (rand()%RANGE+1.0)/RANGE;
for (ans=i=0; i<n; ++i)
ans += sqrt( (cur.x-p[i].x)*(cur.x-p[i].x) + (cur.y-p[i].y)*(cur.y-p[i].y) );
while ( step>eps ){
for (theta=0; theta<2*pi+eps; theta+=delta) { //在该点附近产生8个点
pt.x = cur.x + step*cos(theta);
pt.y = cur.y + step*sin(theta);
for (tmp=i=0; i<n; ++i)
tmp += sqrt( (pt.x-p[i].x)*(pt.x-p[i].x) + (pt.y-p[i].y)*(pt.y-p[i].y) );
if (tmp<ans) {ans=tmp;cur=pt;}
}
step *= rate;
}
if (ans<ret) ret=ans;
}
return ret+0.5;
}
int main(){
scanf("%d", &n);
for (int i=0;i<n;++i) scanf("%lf %lf", &p[i].x , &p[i].y);
printf("%d\n", SA());
return 0;
}
//算法:2.牛顿迭代。
#include <stdio.h>
#include <math.h>
#define MAXN 110
#define eps 1e-3
struct Point{
double x,y;
}p[MAXN];
int n;
double dist[MAXN];
//令f(x,y)=Sigma( sqrt((x-xi)^2 + (y-yi)^2) ),本题要求f(x,y)的最小值。
//此函数是凸的。可验证极小值唯一。且当p0,p1..pn-1不在一条线上时,极小值点唯一。
//于是只需求得f'x=f'y=0的点即可。用牛顿迭代法解此二元方程组。
//简记disti=sqrt( (x-xi)^2 + (y-yi)^2 ),求导得:
//fx=Sigma( (x-xi)/disti) ); fy=Sigma( (y-yi)/disti );
//fxx=Sigma( (y-yi)^2/(dist^3) ); fyy=Sigma( (x-xi)^2/(dist^3) );
//fxy = -Sigma( (x-xi)*(y-yi)/(disti^3) );
//迭代公式为:
//x_k+1 = x_k + (fx*fyy-fy*fxy)/(fxy*fxy-fxx*fyy)|(x_k,y_k);
//y_k+1 = y_k + (fy*fxx-fx*fxy)/(fxy*fxy-fxx*fyy)|(x_k,y_k);
int Newton(){
double ret=1e40,ans=1e40,tmp=0;
double fx,fxx,fy,fyy,fxy; //各阶偏导数的值
double tx,ty,d3; //三个简化计算的临时变量
int i;
Point cur,pt;
cur.x = cur.y = 0;
for (i=0;i<n;++i) { cur.x += p[i].x; cur.y += p[i].y;}
cur.x = cur.x/n + 0.5; //取平均点作为迭代初值(+0.5是为了避免初值落在某个点上)
cur.y = cur.y/n + 0.5;
for (ans=i=0; i<n; ++i)
ans += dist[i] = sqrt( (cur.x-p[i].x)*(cur.x-p[i].x) + (cur.y-p[i].y)*(cur.y-p[i].y) );
while (1){
fx=fy=fxx=fyy=fxy=0; //计算各阶偏导数
for (i=0; i<n; ++i) {
tx = cur.x - p[i].x;
ty = cur.y - p[i].y;
d3 = dist[i]*dist[i]*dist[i];
fx += tx/dist[i];
fy += ty/dist[i];
fxx += ty*ty/d3;
fyy += tx*tx/d3;
fxy -= tx*ty/d3; //注意这里是减号!!(虽然数据太弱加号也能AC……=。=)
}
tmp = fxy*fxy-fxx*fyy; //迭代出新解
pt.x = cur.x + (fx*fyy-fy*fxy)/tmp;
pt.y = cur.y + (fy*fxx-fx*fxy)/tmp;
for (tmp=i=0; i<n; ++i)
tmp += dist[i] = sqrt( (pt.x-p[i].x)*(pt.x-p[i].x) + (pt.y-p[i].y)*(pt.y-p[i].y) );
if (fabs(tmp-ans)<eps) break;
else {cur=pt;ans=tmp;}
}
return ans+0.5;
}
int main(){
scanf("%d", &n);
for (int i=0;i<n;++i) scanf("%lf %lf", &p[i].x , &p[i].y);
printf("%d\n", Newton());
return 0;
}
Followed by: Post your reply here: |
All Rights Reserved 2003-2013 Ying Fuchen,Xu Pengcheng,Xie Di
Any problem, Please Contact Administrator