rand.c

Solaris环境下的数据挖掘算法：birch聚类算法。该算法适用于对大量数据的挖掘。

/****************************************************************************

     Copyright (c) 1993 1994 1995
     By Miron Livny, Madison, Wisconsin
     All Rights Reserved.

     UNDER NO CIRCUMSTANCES IS THIS PROGRAM TO BE COPIED OR DISTRIBUTED
     WITHOUT PERMISSION OF MIRON LIVNY

     obtained and modified on 3/15/96 by Tian Zhang
****************************************************************************/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <malloc.h>

#include "rand.h"

#define  DRAW	s1 = s1* 1103515245 + 12345
#define	 ABS(x)	(x & 0x7fffffff)
#define	 MASK(x) ABS(x)
#define  MAX_INT 2147483648

int	start = 1;
void	init();
int	binary_search(int, int, float*, float);

enum SeedType  { InitialSeed, LastSeed, NewSeed } ;
typedef int *IntArrayPtr;
typedef char *BooleanArrayPtr;

#define M1   2147483563
#define M2   2147483399
#define A1   40014
#define A2   40692
#define w    30 
#define A1W   1033780774
#define A2W   1494757890     
#define MAXG  1073741824    
#define UnSurM1   4.656613057e-10

IntArrayPtr Ig1, Ig2, Lg1, Lg2, Cg1, Cg2 ;
IntArrayPtr Antithetic ; 
int   A1VW, A2VW,  NumberOfGenerators;
float   ln10;

#define   H  32768
/*=================================================================*/
int MultMod( int a, int s, int m ) {
int   a0;
int   a1, q, qh, rh, k, p ;
   if ( a < H ) {
      a0 = a;  p = 0;
    }
   else {
      a1 = a / H;   a0 = a - H * a1;
      qh = m / H;   rh = m - H * qh;
      if ( a1 >= H ) { 
         a1 = a1 - H;  k = s / qh;
         p = H * (s - k * qh) - k * rh;
         while ( p < 0 ) {  p = p + m;  }
         }
	else {
         p = 0;
      };
      if ( a1 != 0 ) {
         q = m / a1;  k = s / q;
         p = p - k * (m - a1 * q);
         if ( p > 0 ) {  p = p - m ; };
         p = p + a1 * (s - k * q);
         while ( p < 0 ) {  p = p + m;  };
      };          
      k = p / qh;
      p = H * (p - k * qh) - k * rh;
      while ( p < 0 ) {  p = p + m;  };
   };            
   if ( a0 != 0 ) {
      q = m / a0;  k = s / q;
      p = p - k * (m - a0 * q);
      if ( p > 0 ) {  p = p - m;  };
      p = p + a0 * (s - k * q);
      while ( p < 0 ) {  p = p + m;  };
   };
   return p;
};

/*=================================================================*/
void GetState( int g , int  *S1, int *S2 ) {
   *S1 = Cg1[g];  *S2 = Cg2[g];
};

/*=================================================================*/
void AdvanceState( int g , int k ){
    int B1, B2, N ;
    int i ;
    B1 = A1;  B2 = A2;
    for ( i = 1; i <= k-1 ; i++) {
        B1 = MultMod( B1, A1, M1 );
        B2 = MultMod( B2, A2, M2 );
    };
   /* B1 = A1 ** k   and   B2 = A2 ** k.    */
   Cg1[g] = MultMod( B1, Cg1[g], M1 ) ;
   Cg2[g] = MultMod( B2, Cg2[g], M2 ) ;
   N = MultMod( A1W, Lg1[g], M1 ) ;
   while ( Cg1[g] >= N ) {
	  Lg1[g] = N ;
	  N = MultMod( A1W, Lg1[g], M1 ) ;
   };
   N = MultMod( A2W, Lg2[g], M2 ) ;
   while ( Cg2[g] >= N ) {
	  Lg2[g] = N ;
	  N = MultMod( A2W, Lg2[g], M2 );
       } ;
    };

/*=================================================================*/
void InitGenerator( int g , SeedType Where ){
    switch (Where) {
       case InitialSeed  :  
          Lg1[g] = Ig1[g];
          Lg2[g] = Ig2[g];
	  break;
       case NewSeed :
           Lg1[g] = MultMod( A1W, Lg1[g], M1 );
           Lg2[g] = MultMod( A2W, Lg2[g], M2 );
           break;
       case LastSeed :
           break;
       };
   Cg1[g] = Lg1[g];  
   Cg2[g] = Lg2[g];
   };

/*=================================================================*/
void SetSeed( int g , int S1, int S2 ) {
   Ig1[g] = S1;  Ig2[g] = S2;
   InitGenerator( g, InitialSeed );
};

/*=================================================================*/
void ResetGenerators( SeedType Where ){
    for (int g = 1;g <= NumberOfGenerators; g++) InitGenerator(g,Where);
};

/*=================================================================*/
void SetInitialSeed( int S1, int S2 ){
   Ig1[1] = S1;  Ig2[1] = S2;
   InitGenerator( 1, InitialSeed );
    for ( int g = 2 ; g <= NumberOfGenerators; g++ ) {
      Ig1[g] = MultMod( A1VW, Ig1[g-1], M1 );
      Ig2[g] = MultMod( A2VW, Ig2[g-1], M2 );
      InitGenerator( g, InitialSeed );
   }
}; 

/*=================================================================*/
void SetAntithetic ( int g, int B ) {
   Antithetic[g] = B;
};

/*=================================================================*/
int Random( int g ){
    int Z, k, s1, s2 ;
    s1 = Cg1[g];  s2 = Cg2[g];
    k  = s1 / 53668;
    s1 = 40014 * (s1 - k * 53668) - k * 12211;
    if ( s1 < 0 )  s1 = s1 + 2147483563 ;
    k  = s2 / 52774;
    s2 = 40692 * (s2 - k * 52774) - k * 3791;
    if ( s2 < 0 )  s2 = s2 + 2147483399 ;
    Cg1[g] = s1;  Cg2[g] = s2;
    Z = s1 - s2;
    if ( Z < 1 )  Z = Z + 2147483562 ;
    if ( Antithetic[g] )  Z = 2147483561 - Z;
    return Z ;
};

/*=================================================================*/
double Uniform( int g ) {
    double ReturnValue;
    ReturnValue = float (Random(g)) * UnSurM1 ;
    return ReturnValue ;

};

/*=================================================================*/
void CreateGenerators( int gnum){ 
    int v, i; 
    double longW; 
    float realV;
    delete(Ig1) ; delete(Ig2) ;
    delete(Lg1) ; delete(Lg2) ;
    delete(Cg1) ; delete(Cg2) ;
    delete(Antithetic); 
    /*  Memory Allocation   */ 
    NumberOfGenerators = gnum ; v = NumberOfGenerators * sizeof(int); 
    Ig1        = new int[gnum+1] ;
    Ig2        = new int[gnum+1] ; 
    Lg1        = new int[gnum+1] ; 
    Lg2        = new int[gnum+1] ;
    Cg1        = new int[gnum+1] ; 
    Cg2        = new int[gnum+1] ;
    Antithetic = new int[gnum+1] ;
    /*  v : Number of sub-sequences in each generator  */
    longW = w;
    realV = (longW -log10(NumberOfGenerators) ) / log10( 2.0 ) ;
    v = int (realV);
    /*  Calculate A1VW, A2VW  */
    A1VW = A1W ; A2VW = A2W ;
    for ( i = 1 ; i <=  v ; i++ ) {
         A1VW = MultMod( A1VW, A1VW, M1 ) ;
         A2VW = MultMod( A2VW, A2VW, M2 ) ;
        } ;
    for ( i = 1 ; i <= NumberOfGenerators ; i++) Antithetic[i] = 0;
   /*  Initialize the generators.
	   You can change the initial seeds afterwards by explicitly 
	   using the SetInitialSeed routine in your code  */
   SetInitialSeed( 1234567890, 123456789 );
};

/*=================================================================*/
int BddRandom( int g , int LB, int UB) {
int ReturnValue ;
float realUni;
   realUni = Uniform( g );
   ReturnValue = int(float( UB - LB ) * realUni )  + LB ;
   return ReturnValue;
};

/*=================================================================*/
void NewSeeds( int consume ) {
IntArrayPtr  shuffle;
int   i, temp, g, ndx ;
   shuffle = new int[NumberOfGenerators+1];
   g = consume % NumberOfGenerators ;
   AdvanceState( g, consume ) ;
    for ( i = 1 ; i <= NumberOfGenerators - 1; i++ ) {
	  shuffle[i] = BddRandom( g, i + 1, NumberOfGenerators + 1 ) ;
   } ;
    for ( i = 1 ; i <=  NumberOfGenerators - 1 ; i++ ) {
	  ndx = shuffle[i] ;
	  temp = Ig1[ndx] ;
	  Ig1[ndx] = Ig1[i] ;
	  Ig1[i] = temp ;
	  temp = Ig2[ndx] ;
	  Ig2[ndx] = Ig2[i] ;
	  Ig2[i] = temp ;
      InitGenerator( i, InitialSeed ); 
   } ;
   InitGenerator( NumberOfGenerators, InitialSeed ) ;
   delete(shuffle);
};

/*=================================================================*/
void init(){
int i; 

    ln10 = log(10.0);
/************

     By default set the number of generators to 1024. 
	 If the user needs more or less, he can use the routine 
	                      CreateGenerators( gnum ), 
	 where gnum is the number of the desired generators.

************/
   NumberOfGenerators = 1024 ;
   Ig1        = new int[NumberOfGenerators+1];
   Ig2        = new int[NumberOfGenerators+1];
   Lg1        = new int[NumberOfGenerators+1];
   Lg2        = new int[NumberOfGenerators+1];
   Cg1        = new int[NumberOfGenerators+1];
   Cg2        = new int[NumberOfGenerators+1];
   Antithetic = new int[NumberOfGenerators+1];

   A1VW = 2082007225 ;  A2VW = 784306273 ;
   /*  Initialize the Antithetic values  */
    for ( i = 1 ; i <=  NumberOfGenerators ; i++) {
      Antithetic[i] = 0;
   };

   /* Initialize the generators */
   SetInitialSeed(1234567890, 123456789);
};


// ====================================================================================

RandomStream::RandomStream(){
    count     = 0;
    statFlag  = 2;
    sum       = 0;
    sumOfSqu  = 0;
    if (start) { init(); start = 0;};
    };

IntRandomStream::IntRandomStream(){
    sampleMax = -MAX_INT;
    sampleMin = 1000000;
    };

FloatRandomStream::FloatRandomStream(){
    sampleMax = -1.0e29;
    sampleMin = 1.0e29;
    };

void IntRandomStream::Report() {
	float mean;
	float fCount = count ;
	mean = sum/count;
	printf(" %f (avg) ", mean);
	printf("%f (std) ", 
	    sqrt((fCount*sumOfSqu - sum*sum)/(fCount*(fCount - 1))));
	printf("%d (count)\n ", count);
    };

void FloatRandomStream::Report() {
	float mean;
	float fCount = count ;
	mean = sum/count;
	printf(" %f (avg) ", mean);
	printf("%f (std) ", 
	    sqrt((fCount*sumOfSqu - sum*sum)/(fCount*(fCount - 1))));
	printf("%d (count)\n ", count);
    };

void IntRandomStream::CollectStat(int value) {
    count++;