text1.cpp

该程序运用退火法解决了经典的圆排列问题

#include <stdio.h>
#include <math.h>
#include<iostream.h>
#include<stdlib.h>
#include <conio.h> 
#include <time.h> 
#define TFACTR  0.9
#define MAX 10000.0
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)
#define N 96
float weight[N+1][N+1];
int r[N+1]={0,
1,2,2,2,2,6,7,14,14,16,18,18,19,20,20,20,20,21,21,24,26,28,28,29,29,30,31,31,31,35,35,36,38,38,38,41,42,42,43,45,45,45,46,47,47,49,50,51,53,53,54,55,55,57,60,61,61,62,63,65,67,67,68,69,69,70,70,71,71,73,74,74,75,77,80,80,82,82,83,83,84,84,84,85,86,87,88,89,90,90,90,90,92,93,94,94
};
int a;
int b;
float l[N+1];
float x[N+1];
int p[N+1];
void anneal (int ncity)
{ 
float changdu(int *);
int randm();
int randj(int);
int irbit(unsigned long *iseed);
int metrop(float de,float t);
float ran(long *idum);
float reverse_len();
float transe_len();
int ans,nlimit;
int i,j,nsucc,idec;
unsigned long k;
unsigned long iseed;
float path,de,T;
nlimit=100*ncity;/*max succ try numbers */
r[1];
path=changdu(r);
iseed=111;
T=2.8;	
for(j=1;j<=15;j++)
{
     nsucc=0;
     for(k=1;k<=20000;k++)
	 {
	      idec=irbit(&iseed);
	      if(idec==0)
		  {
			 de=transe_len();
			 ans=metrop(de,T);
			 if (ans)
			 {
				 ++nsucc;
				 path+=de;				  
				 for(i=1;i<=N;i++)
					 r[i]=p[i];
				 for(i=1;i<N;i++)
					 l[i]=2*sqrt(r[i]*r[i+1]);
				 l[N]=r[1]+r[N];
			 }
		  }
		  else
		  {
			  de=reverse_len();
			  ans=metrop(de,T);
			  if (ans) {++nsucc;
			  path+=de;				  
			  for(i=1;i<=N;i++)
				  r[i]=p[i];
			  for(i=1;i<N;i++)
				  l[i]=2*sqrt(r[i]*r[i+1]);
			  l[N]=r[1]+r[N];
			  }
		  }
		  if(nsucc>=nlimit) break;
	 }
	 printf("\nT=%10.6f,path len=%10.6f,j=%d",T,path,j);
	 printf("successful moves: %6d\n",nsucc);
	 for(i=1;i<=N;i++) printf("%5d",r[i]);
	 r[1];
	 r[2];
	 T*=TFACTR;
	 //if(nsucc==0) return;
}
}

float changdu(int *r)
	{
		float cx,low,high;
		int i,j;
		for(i=0;i<=N;i++)x[i]=0;
		x[2]=2*sqrt(r[1]*r[2]);
		for(i=3;i<=N;i++)
			for(j=1;j<=i-1;j++)
			{
				cx=x[j]+2*sqrt(r[j]*r[i]);
				if(cx>x[i])x[i]=cx;
			}
		low=(-1)*r[1];high=x[N]+r[N];
		for(i=2;i<N;i++)
		{
			if(x[i]-r[i]<low)low=x[i]-r[i];
			if(x[i]+r[i]>high)high=x[i]+r[i];
		}
		cx=high-low;
		return cx;
	}
int randj(int k)
{//求随机数j（weight[][],rand()）
	int m,//概率的分子(随机数)
		n,//概率的分母
		j;
	float t=0.0,dmax=0;
	for(j=1;j<=N;j++)if(dmax<weight[k][j])dmax=weight[k][j];
	weight[k][k]=dmax;
	for(j=1;j<=N;j++)t+=(dmax-weight[k][j]);
	n=floor(t);
	srand((unsigned)time(NULL)*100); 
	m=rand()%n;
	t=0;
	for(j=1;j<=N;j++)
	{
		if(t<=m&&m<(t+dmax-weight[k][j]))return j;		
		t+=(dmax-weight[k][j]);	
	}
	return j;
}

int randi(void)
{
	int m,//概率的分子（随机数）
		n,//概率的分母
		i;
	float t=0.0;
	for(i=1;i<=N;i++)t+=l[i];
	n=(int)t;
	srand((unsigned)time(NULL)*100); 
	m=rand()%n;	
	for( i=1;i<=N;i++)
	{
		t=0.0;
		for(int j=1;j<i;j++)t+=l[j];
		if(floor(t)<=m&&m<floor(t+l[i]))
			return i;
	}
	return i;
}

float reverse_len(void)
{//将i和j间的数（包括i和j）翻转（randi（），randj（），r[],changdu()）
	int i,j,k;
    float  de;
    i=randi();
	j=randj(i);
	for(k=1;k<=N;k++)p[k]=r[k];
	if(i<j)
	{
		for(k=i;k<=j;k++)
			p[k]=r[i+j-k];
	}
	else
	{
		for(k=j;k<=i;k++)
			p[k]=r[i+j-k];
	}
	de=changdu(p)-changdu(r);
	return de;
}

float transe_len(void)
{//将j插入到i的后面（randi（），randj（），r[],changdu()）
	int i,j,k;
    float  de;
	i=randi();
	j=randj(i);
	for(k=1;k<=N;k++)p[k]=r[k];
	if(i+1<j)
	{
		p[i+1]=r[j];
		for(int k=i+1;k<j;k++)
	   	p[k+1]=r[k];		
	}
	else if(i>j)
	{
		for(int k=j;k<i;k++)
		p[k]=r[k+1];
		p[i]=r[j];
	}
	de=changdu(p)-changdu(r);
	return de;
}
   	    
 int metrop(float de,float t)
 {
      float ran(long *idum);
      static long gljdum=1;
      return de<0.0||ran(&gljdum)<exp(-de/t);
  }
float ran(long *idum)
{
 static int inext,inextp;
 static long ma[56];
 static int iff=0;
 long mj,mk;
 int i,ii,k;
 if(*idum<0||iff==0) {
		      iff=1;
		      mj=MSEED-(*idum<0?-*idum:*idum);
		      mj%=MBIG;
		      ma[55]=mj;
		      mk=1;
		      for(i=1;i<54;i++){
					ii=(21*i)%55;
					ma[ii]=mk;
					mk=mj-mk;
					if(mk<MZ)mk+=MBIG;
					mj=ma[ii];
			}
			for(k=1;k<=4;k++)
			     for(i=1;i<=55;i++){
				   ma[i]-=ma[1+(i+30)%55];
				   if(ma[i]<MZ) ma[i]+=MBIG;
				   }
			 inext=0;
			 inextp=31;
			 *idum=1;
		     }
    if(++inext==56)inext=1;
    if(++inextp==56) inextp=1;
    mj=ma[inext]-ma[inextp];
    if(mj<MZ) mj+=MBIG;
    ma[inext]=mj;
    return mj*FAC;
    }
 int irbit(unsigned long *iseed)
 {
     unsigned long newbit;
      newbit=(*iseed &131072)>>17
	    ^(*iseed&16)>>4
	      ^(*iseed &2)>>1
	      ^(*iseed &1);
	      *iseed=(*iseed<<1)|newbit;
	      return newbit;
}
void main(void)
{
  int ncity=N;
  int i,j;
for(i=1;i<=N;i++)
for(j=1;j<=N;j++)
if(i!=j)
weight[i][j]=2*sqrt(r[i]*r[j]);
else weight[i][j]=MAX;
for(i=1;i<N;i++)
	  		 l[i]=2*sqrt(r[i]*r[i+1]);
	         l[N]=r[1]+r[N];
anneal(ncity);
}