math_add.c

一个很好用的可执行PSO程序

/*------------------------------------------------------------------- ALEA */
int alea(int min,int max) /* Random integer number between min and max */
{
// Note: "seed" is a global variable
int ir;
float r;

srand(seed); /* Initialize random sequence */

r=rand();
seed=r; /* Prepare reinitialization for the next time */

r=r/RAND_MAX; /* Normally, RAND_MAX = 32767 = 2^31-1 */
r=0.5+min+r*(max-min);
ir=r;
//printf("\n alea %i",ir);
return ir;
}


/*------------------------------------------------------------------- ALEA_DIFF */
int alea_diff(int min,int max, int num) /* Generate randomly a number # num */
{
int	temp1,temp2;

//printf("\n alea_diff, %i,%i,%i",min,max,num);

if (num==min)
	{
	temp1=alea(min+1,max);
	return temp1;
	}
	
if (num==max)
	{
	temp1=alea(min,max-1);
	return temp1;
	}

temp1=alea(min,num-1);
temp2=alea(num+1,max);

if (alea(0,100)<50)
	{
	return temp1;
	}
else
	{
	return temp2;
	}
}

/*------------------------------------------------------------------- ALEA_FLOAT */
float alea_float(float min, float max) /* Random integer number between min and max */
{
// Note: "seed" is a global variable

float r;

srand(seed); /* Initialize random sequence */

r=rand();
seed=r; /* Prepare reinitialization for the next time */

r=r/RAND_MAX; /* Normally, RAND_MAX = 32767 = 2^31-1 */
r=min+r*(max-min);

return r;
}

//======================================================== CHECK_SEQ
 int check_seq(struct seq s)
 // Check whether there is two times the same value
 {
  int check;
  int i,j;

  for (i=0;i<s.size-1;i++)
 {
     for (j=i+1;j<s.size;j++)
     {
        if (s.s[i]==s.s[j])
        {
           printf("\n    ERROR. Same value %i at ranks %i and %i",s.s[i],i,j);
           return 1;
       }
     }
   
 }

  return 0; 

 }


/*------------------------------------------------------------------- FACT */
int fact(int n) /* n! */
{
int i,f;

f=1;
if (n<2) return f;

for (i=2;i<=n;i++)
	f=f*i;

return f;	
}

/*------------------------------------------------------------------- MAX */
float MAX(float a,float b)
{
if (a>b) return a; return b;
}

/*------------------------------------------------------------------- MIN */
float MIN(float a,float b)
{
if (a<b) return a; return b;
}

/*------------------------------------------------------------------- MAX_EIG */
float max_eig(struct coeff coeff, float phi)
{
// Compute the max of the two eigenvalues (general case)

float 	a;
float 	alpha, beta, gamma, delta, eta;
float 	dis1,discrim;
//int		i;
float	maxneprim;
float	nemaxneprim;
float	neprim1,neprim2;

alpha=coeff.c[0];
beta=coeff.c[1];
gamma=coeff.c[2];
delta=coeff.c[3];
eta=coeff.c[4];

//printf("\n phi %f, coeffs:",phi);
//for (i=0;i<5;i++)
	//printf(" %f", coeff.c[i]);

dis1=(eta*phi)*(eta*phi)-4*beta*gamma*phi; //just to have a shorter line
discrim=(dis1+(alpha-delta)*(alpha-delta)+2*eta*phi*(alpha-delta))/4;
a=(alpha+delta-eta*phi)/2;

//printf("\n dis1 %f,  discrim %f, a %f",dis1,discrim,a);

if (discrim>0)
	{
	neprim1=fabs(a+sqrt(discrim));     //norm of the first eigenvalue
	neprim2=fabs(a-sqrt(discrim));     //norm of the second eigenvalue
	maxneprim=MAX(neprim1,neprim2);
	}
else    //In this case the eigenvalues are 'true' complex numbers
	{
	maxneprim=sqrt(a*a+fabs(discrim));
	nemaxneprim=maxneprim;      /*Compute some values, to evaluate convergence
                      			probability  */
	}


return maxneprim;
}

/*------------------------------------------------------------------- MAX_EIG_BASIC */
float max_eig_basic(float phi)
//---------------- Basic case (useful to compute constriction coefficients)
// Basic case (all coeffs equal to 1)
{
struct coeff	c;
int				i;
float 			maxne;

for (i=0;i<5;i++)
	c.c[i]=1;

//Compute the max of the norms of the two eigenvalues e1 and e2
maxne=max_eig(c, phi);

// Note. For phi<=4 ne1=ne2=1. For phi>4, ne1 tends to 0 and ne2 to infinity

return maxne;

}

/*------------------------------------------------------------------- ROTATE */
struct seq rotate(struct seq s,int k)
{
/* Rotate a sequence so that it begins on k
Except the last/extra element. Must be equal to the first one
*/

int i,j;
struct seq st;

if (s.s[0]==k && s.s[s.size]==k) return s;

if (s.s[0]==k && s.s[s.size]!=k)
{
   st=s;
   st.s[s.size]=k; 
   return st;
}

for (i=0;i<s.size;i++)
	{
	if (s.s[i]!=k) continue;
		for (j=i;j<s.size;j++)
			st.s[j-i]=s.s[j];
			
		for (j=0;j<i;j++)
			st.s[j+s.size-i]=s.s[j];
	}
st.s[s.size]=k; /* Needed when the sequence is a particle position */
st.size=s.size;
return st;
}



/*------------------------------------------------------------------- SEQ */
struct seq	seq(int k,int N)
{
/* Generates a sequence of N different integers, with second number equal to k */
int i,j;
struct seq seqt;

seqt.s[1]=k;

for (i=2;i<N;i++)
	{
	
	j=seqt.s[i-1]+1;
	if (j==N) j=1;
	seqt.s[i]=j;
	}
return seqt;
}