From:	MX%"BJ06@C53000.PETROBRAS.ANRJ.BR" 29-SEP-1994 16:43:48.79
To:	MX%"noel@erich.triumf.ca"
CC:	
Subj:	FD3 for DOS (FD3DOS.CPP) - Source code

Return-Path: 
Received: from fpsp.fapesp.br by Erich.Triumf.CA (MX V4.0-1 VAX) with SMTP;
          Thu, 29 Sep 1994 16:42:40 PST
Received: from DECNET-MAIL by fpsp.fapesp.br with PMDF#10108; Thu, 29 Sep 1994
          09:26 BSC (-0300 C)
Date: Thu, 29 Sep 1994 09:26 BSC (-0300 C)
From: "FAUSTO Arinos de Almeida Barbuto (Totxo)"
      
Subject: FD3 for DOS (FD3DOS.CPP) - Source code
To: noel@erich.triumf.ca
Message-ID: <04B80DB3A0009621@fpsp.fapesp.br>
X-Envelope-to: noel@erich.triumf.ca
X-VMS-To: @NOEL
References: ANSP network   HEPnet SPAN Bitnet Internet gateway
Comments: @FPSP.FAPESP.BR - @FPSP.HEPNET - @BRFAPESP.BITNET - .BR gateway

/* BEGIN NOTICE

Copyright (c) 1992 by John Sarraille and Peter DiFalco
(john@ishi.csustan.edu)

Permission to use, copy, modify, and distribute this software
and its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that copyright notice and this permission
notice appear in supporting documentation.

The algorithm used in this program was inspired by the paper
entitled "A Fast Algorithm To Determine Fractal Dimensions By
Box Counting", which was written by Liebovitch and Toth, and
which appeared in the journal "Physics Letters A", volume 141,
pp 386-390, (1989).

This program is not warranteed: use at your own risk.

-------------------------------------------------------
DOS Version created by Fausto Arinos de Almeida Barbuto
(BJ06@C53000.PETROBRAS.ANRJ.BR), September 23, 1994
(Rio de Janeiro, Federal Republic of BRAZIL)
-------------------------------------------------------

END NOTICE */

#include 
#include 
#include "fd.h"
#include 

int  empty_Av(long int);
void push_Av(long int *, long int);
void pop_Av(long int *, long int*);
void create_Av(long int *);
void print_Av(long int);

int get_e_dim(char *);
void max_min(char*, double *, double *);
void rad_pass(int, int);
void radixsort(void);
void sweep(ulong[], double[], double[], double[]);
float fracdim();

int empty_Q(QHeader);
void en_Q(QHeader *, long int);
void de_Q(QHeader *, long int *);
void create_Q(QHeader *);
void transf_Q(QHeader *, QHeader *);
void print_Q(QHeader *);
void GetDims(double[], double[], double[], ulong, double*, double*, double*);
void findMark(ulong *, ulong[]);
void fitLSqrLine (long, long, double[], double[], double *, double *);

int empty_Av(long int avail)
//       long int avail;
{
  return (avail == -1);
}

/*
    Puts the info pointed to by "pointer" into the avail list.  Use
    this after dequeuing an element if you want to recycle it.  Note the
    need to access the global array "next_after", which holds all the
    pointers.
*/
void push_Av(long int *p_avail, long int pointer)
//	long int *p_avail, pointer ;
{
  extern long int  *next_after;   /*  array for pointers to data. */
  next_after[pointer] = *p_avail;
  *p_avail = pointer;
}

/*
   Pops the info pointed to by "p_pointer" from the avail list.  Use
   this if you need an element to place on a queue.
*/
void pop_Av(long int *p_avail, long int *p_pointer)
//	long int *p_avail, *p_pointer;
{
  extern long int  *next_after;   /*  array for pointers to data. */
  *p_pointer = *p_avail;
  *p_avail = next_after[*p_avail];
}

/*
    Creates an avail list with room for all the data, plus two more
    elements to be used as markers.
*/
void create_Av(long int *p_avail)
//	long int *p_avail;
{
  extern long int  *next_after;   /*  array for pointers to data. */
  extern ulong dataLines;        /* number of data points in input file */
  unsigned long int i;
  for (i=0;iQfront = pointer;
  else
     next_after[p_QHeader->Qrear] = pointer;
  p_QHeader->Qrear = pointer;
}

/*
    This Dequeuing function does NOT check to see if the queue is
    empty first, nor does it return the dequeued info to the avail
    list.  The caller MUST take care of these things!  Note also
    that the function needs to access the array called "next_after"
    globally.
*/
void de_Q(QHeader *p_QHeader,long int *p_pointer)
// 	QHeader *p_QHeader;		/* Queue to be dequeued. */
//	long int *p_pointer;	/* Pointer to the data to be returned */
				/* to the calling function. */	
{
  extern long int  *next_after;   /*  array for pointers to data. */
  *p_pointer = p_QHeader->Qfront;
  if ((p_QHeader->Qfront) == (p_QHeader->Qrear))
     p_QHeader->Qrear = -1 ;
  p_QHeader->Qfront = next_after[p_QHeader->Qfront];
}

void create_Q (QHeader *p_QHeader)
//	QHeader *p_QHeader;	/* Queue to be created. */
{
  p_QHeader->Qfront = -1;
  p_QHeader->Qrear = -1;
}

/*
    This is a function to rapidly move an "element" from the source
    to the target queue.  It bypasses the actions of placing the info
    into the avail list, and taking it back out again.  It also
    avoids testing to see if the source queue is being emptied, and
    testing to see if the target queue is empty before the enqueue
    is done.

    Such a function is useful because there will be a great many of
    these operations done in the radix sort of a large set.  For
    further time-savings, one should think about putting this code
    in-line, so that the overhead of function-calling can be saved.

    USE with EXTREME CAUTION.  The effect of the function will be
    INCORRECT if the source has one element left, or the target is
    empty!

    Note the reliance on the global array "next_after".
*/
void transf_Q(QHeader *p_Q1, QHeader *p_Q2)
//        QHeader *p_Q1, *p_Q2;  /* source and target queues */
{
  extern long int  *next_after;   /*  array for pointers to data. */
  long int temp;
  temp = p_Q1->Qfront;                     /* save front of Q1 */
  p_Q1->Qfront = next_after[p_Q1->Qfront];  /* fast dequeue */
  next_after[temp] = -1;               /* fast enqueue */
  next_after[p_Q2->Qrear] = temp;        /* fast enqueue */
  p_Q2->Qrear = temp;                    /* fast enqueue */
}

/*
    This may come in handy when de-bugging.  Note the reliance on
    the global arrays "next_after" and "data".
*/
void print_Q(QHeader *p_Q)
//     QHeader *p_Q;
{
  extern long int  *next_after;   /*  array for pointers to data. */
  extern int  *marker;    /*  To mark queues during radix sort */
  extern int   embed_dim;         /* How many coordinates points have. */
  extern ulong **data;	/* array of data points */  
  long int temp;
  int coord;
  printf("Starting to print the queue:\n\n");
  printf("Index\tNextField\tMarker?\t\tCoords\n\n");
  temp = p_Q->Qfront;
  while (temp != -1)
    {
        printf("%ld\t%ld\t\t%d\t", temp, next_after[temp], marker[temp]);
	for (coord=0;coord *pmax)
                        *pmax = temp;
                else
                        if(temp < *pmin)
                                *pmin = temp;
 if (debugging) printf ("Pmin is now: %lf, Pmax is now: %lf\n", *pmin, *pmax);
        }
        /* check to see if the maximum equals the minimum -- this is a
           degenerate case -- or it might mean that the input file is
}          faulty. */
        if(*pmax == *pmin)
        {
          printf ("The input file, %s, is confusing!\n\n", userinputfile);
          printf ("Either all the points in %s have the same ", userinputfile);
          printf ("coordinates,\n");
          printf ("(THE FRACTAL DIMENSION IS ZERO IN THIS CASE.)\n\n");
          printf ("or %s is simply of the wrong form for an\n",userinputfile);
          printf ("input file to this program -- please check.\n");
          exit(-1);
        }
        /*  close the input file */
        fclose(infile);
}
/* ################################################################## */
/* ################################################################## */
/*
   This procedure performs one pass of the radix sort.  It assumes that the
   queues each have a marker in front at the time of the call, and this is
   the condition the procedure LEAVES the queues in when it terminates.
*/
void rad_pass(int coord, int bitpos)
//        int coord,   /*  The coordinate we are doing this pass on */
//            bitpos;  /*  The bit position we are doing this pass on */
{  extern int  *marker;           /*  To mark queues during radix sort */
   extern QHeader Q[2];       /* array of two queues for the radix sort */
   extern ulong **data;   /* initial array of data points */
   int queue_num;
   long int index ;
    /*  Move the markers to the rear of the queues */
 if (debugging)
    printf ("Starting pass on bit %d of coord %d.\n",bitpos,coord);
 if (debugging) printf ("Moving markers to the rear of queues.\n");
 for (queue_num=0;queue_num<2;queue_num++) {
     de_Q(&(Q[queue_num]),&index);
     en_Q(&(Q[queue_num]),index);
 }
 if (debugging) print_Q(&Q[0]);
 if (debugging) print_Q(&Q[1]);

    /*  Move all non-marker elements to the appropriate queue. */
 if (debugging) printf ("Starting main part of pass.\n");
 for (queue_num=0;queue_num<2;queue_num++)
   {    /* Peek first to see if Qfront is a marker -- this violates
           information hiding, and would be "cleaner" if done with a
           procedure in the queue package.  */
     while ( marker[Q[queue_num].Qfront] == 0 )
       {     /*  Peeking again! */
         if (IS_A_ONE(data[coord][Q[queue_num].Qfront],bitpos))
               /* Directly transfer from the source to target queue */
          transf_Q( &(Q[queue_num]),&(Q[1]) );
         else transf_Q( &(Q[queue_num]),&(Q[0]) ) ;
        if (debugging) printf("A queue transfer is done.  \n");
        if (debugging) print_Q(&Q[0]);
        if (debugging) print_Q(&Q[1]);
       }  }  }

/* ################################################################## */
/* ################################################################## */
/*
THIS SORT IS TO BE USED DIRECTLY AFTER FDDRIVER DOES THE INITIAL LOADING OF
THE DATA INTO THE QUEUES.  IT LEAVES THE DATA ESSENTIALLY SORTED, WITH ALL
THE DATA WHOSE X-COORD STARTS WITH 0 IN Q[0], IN ORDER, AND ALL THE DATA
WHOSE X-COORD STARTS WITH 1 IN Q[1], IN ORDER.  THUS THE SWEEP THAT COMES
NEXT MUST TRAVERSE Q[0], AND THEN Q	[1].
*/
void radixsort()
{
  extern QHeader Q[2];       /* ARRAY OF TWO QUEUES FOR THE RADIX SORT */
  extern int   embed_dim;      /* HOW MANY COORDINATES POINTS HAVE. */
  extern long int  avail;             /*  AVAIL LIST POINTER */
  int  bitpos, coord, queue_num;
  long int index;

    /* FINISH UP ON THE ZERO-BIT -- LOADING DATA TOOK CARE OF FIRST PASS. */
  for (coord=embed_dim-2;coord>=0;coord--) rad_pass(coord,0);

      /* NOW SORT ON THE REST OF THE BITS. */
  for (bitpos=1;bitpos=0;coord--) rad_pass(coord,bitpos);

      /*  LASTLY, GET THE MARKERS OUT OF THE QUEUES. */
  for (queue_num=0;queue_num<2;queue_num++)
    {
      de_Q(&(Q[queue_num]),&index);
      push_Av(&avail,index);
    }}
/* ################################################################## */
/* ################################################################## */
/*
THIS PROCEDURE TRAVERSES THE SORTED DATA POINT LIST, AND EXTRACTS THE
INFORMATION NEEDED TO COMPUTE THE CAPACITY, INFORMATION, AND CORRELATION
DIMENSIONS OF THE DATA.  DATA IS CORRECTED ON THE FLY DURING THE TRAVERSAL,
AND THEN "MASSAGED".  AFTER "SWEEP" RUNS, WE NEED ONLY FIT A SLOPE TO THE
DATA TO OBTAIN THE FRACTAL DIMENSIONS.
*/
void sweep(ulong boxCountS[], double negLogBoxCountS[], double logSumSqrFreqS[],
	   double informationS[])

//   ulong   boxCountS[numbits+1];  /* COUNT BOXES OF EACH SIZE */

//   double
	       /* FOR EACH BOX SIZE #s, negLogBoxCountS[s] WILL EVENTUALLY BE
		  SET TO THE NEGATIVE OF THE LOG (BASE TWO) OF THE NUMBER OF
		  BOXES OF SIZE #s THAT ARE OCCUPIED BY ONE OR MORE DATA
		  POINTS.
		*/

//	   negLogBoxCountS[numbits+1],

	       /* FOR EACH BOX SIZE #s, logSumSqrFreqS[s] WILL EVENTUALLY BE
		  SET TO THE LOG (BASE TWO) OF THE PROBABILITY THAT TWO
		  RANDOMLY CHOSEN DATA POINTS ARE IN THE SAME BOX OF SIZE
		  #s.
	       */

//	   logSumSqrFreqS[numbits+1],

	      /* FOR EACH BOX SIZE #s, informationS[s] WILL EVENTUALLY BE
		 SET TO THE INFORMATION (BASE TWO) IN THE DISTRIBUTION OF
		 DATA POINTS IN THE BOXES OF SIZE #s.
	      */
//	   informationS[numbits+1] ;

{ extern int  embed_dim;         /* HOW MANY COORDINATES POINTS HAVE. */
  extern ulong **data,           /* INITIAL ARRAY OF DATA POINTS */
	       dataLines ;
  extern QHeader Q[2];           /* ARRAY OF TWO QUEUES FOR THE RADIX SORT */
  extern ulong  *diff_test;      /* XOR'S OF PAIRS OF COORDINATES */
  extern long int  *next_after;  /*  ARRAY FOR POINTERS TO DATA. */

  int       bitpos, countpos, coord, queue_num, found;
  ulong     pointCount[numbits+1] ;
  long int  current, previous;
  double    sumSqrFreq[numbits+1], freq, log2=log(2.0) ;

    /* GET A POINTER TO THE FIRST DATA VALUE */
  if    (Q[0].Qfront != -1) previous = Q[0].Qfront;
  else  previous = Q[1].Qfront;

     /* INIT boxCountS, pointCountS, sumSqrFreq, AND informationS.*/
  for(countpos=0;countpos<=numbits;countpos++) {
       boxCountS[countpos]=1;
       sumSqrFreq[countpos]=0.0;
       pointCount[countpos]=1;
       informationS[countpos]=0.0;
  }

  for (queue_num=0;queue_num<2;queue_num++)
  { current = Q[queue_num].Qfront;
    while (current != -1)
    { found = 0 ;
        /* START BY LOOKING AT THE BIGGEST BOX SIZE */
      bitpos=numbits-1;
      for (coord=embed_dim-1;coord>=0;coord--)
        diff_test[coord] = data[coord][previous]^ data[coord][current];
      do {coord = embed_dim - 1;
          do {/* IF THE CURRENT POINT AND PREVIOUS POINTS ARE IN DIFFERENT
	         BOXES OF THIS SIZE, */
	       if ( IS_A_ONE(diff_test[coord],bitpos) )
               {/* THEN THE CURRENT POINT IS IN NEW BOXES OF ALL SMALLER
		   SIZES TOO, AND STILL IN THE SAME BOXES OF LARGER SIZES,
		   SO ... */
	        for (countpos=bitpos;countpos>=0;countpos--)
		  {/* CALCULATE FREQUENCY OF POINTS IN THE BOX, ASSUMING FOR
		      NOW THAT THE NUMBER OF DATA LINES IN THE INPUT FILE IS
		      THE NUMBER OF DISTINCT POINTS IN THE DATA SET.  WE
		      ADJUST THIS AT THE END OF THIS FUNCTION. */
                   if (debugging)
		    printf("pointCount[%d] is %d...\n", countpos, pointCount[countpos]);
		   freq = pointCount[countpos]/(double)dataLines ;
		     /* WE WILL ENCOUNTER NO MORE OF THE POINTS IN THE BOX
			WE JUST LEFT (THE SPECIAL ORDERING OF THE SORT WE
			USED ABOVE GUARANTEES THIS!), SO WE COMPUTE WHAT
			THIS BOX CONTRIBUTES TO THE RUNNING SUMS. */
                   sumSqrFreq[countpos] += (freq * freq) ;
                   informationS[countpos] += ( freq*log(freq)/log2 ) ;
 		     /* WE HAVE GOTTEN INTO A NEW BOX AT THIS LEVEL, SO WE
			REFLECT THE NEW BOX IN THE COUNT */
                   boxCountS[countpos]++ ;
                     /* SINCE WE HAVE A NEW BOX AT THIS LEVEL, THERE IS ONLY
		        ONE KNOWN POINT IN IT SO FAR -- THE CURRENT POINT */
		   pointCount[countpos]=1;
                  }
                for (countpos=bitpos+1;countpos<=numbits;countpos++)
		    /* THE CURRENT POINT IS IN THE BOXES AT THESE LEVELS, SO
		       JUST INCREMENT THE POINT COUNTER.  */
                  pointCount[countpos]++ ;
                found = 1;
               }
               else coord-- ;
             }
          while ( (found == 0) && (coord > -1) );
          bitpos-- ;
         }
      while ( (found == 0) && (bitpos > -1) );
      previous = current; current = next_after[current];
    } }
    
    /* NOW ADD IN THE CONTRIBUTION DUE TO THE COUNTS REMAINING AFTER THE
       LAST POINT HAS BEEN FOUND, RENORMALIZE WITH BOXCOUNT[0], AND MASSAGE
       THE RAW DATA FROM THE TRAVERSAL SO THAT IS IS READY FOR THE LEAST
       SQUARES FITTING. */
       
  for (countpos=numbits;countpos>=0;countpos--)
  {
    negLogBoxCountS[countpos] = -log((double)boxCountS[countpos])/log(2.0);

    if (debugging)
      printf("pointCount[%d] is %d...\n", countpos, pointCount[countpos]);
    freq = pointCount[countpos]/(double)dataLines ;

    sumSqrFreq[countpos] += (freq * freq) ;
    sumSqrFreq[countpos] *= (dataLines/(double) boxCountS[0]) ;
    sumSqrFreq[countpos] *= (dataLines/(double) boxCountS[0]) ;

       /* sumSqrFreq[countpos] NOW CONTAINS THE SUM OF THE SQUARES OF THE
	  FREQUENCIES OF POINTS IN ALL OCCUPIED BOXES OF THE SIZE
	  CORRESPONDING TO countpos. */
	  
    logSumSqrFreqS[countpos] = log(sumSqrFreq[countpos])/log(2.0) ;
    informationS[countpos] += ( freq*log(freq)/log(2.0) ) ;
    informationS[countpos] *= (dataLines/(double)boxCountS[0]) ;
    informationS[countpos] +=
       ( log((double)dataLines)-log((double)boxCountS[0]) ) / log(2.0) ;

       /* information[countpos] NOW CONTAINS THE INFORMATION SUM FOR ALL THE
          OCCUPIED BOXES OF THIS SIZE. */

   }
}
/* ################################################################## */
/* ################################################################## */
  /*  FIT LEAST SQUARE LINE TO DATA IN X,Y.  NO PROTECTION AGAINST OVERFLOW
      HERE.  IT IS ASSUMED THAT LAST > FIRST AND THAT THE X'S ARE NOT ALL THE
      SAME -- ELSE DIVISION BY ZERO WILL OCCUR.  */
void fitLSqrLine (long first, long last, double X[], double Y[],
		  double *slopePtr, double *interceptPtr)
//     long   first, last ;
//     double X[], Y[], *slopePtr, *interceptPtr ;
{
  int    index , pointCount ;
  double Xsum=0, Ysum=0, XYsum=0, XXsum=0, Xmean=0, Ymean=0,
         Xtemp, Ytemp;
  for (index=first; index<=last; index++)
    { Xtemp = X[index]; Ytemp = Y[index];
      Xsum += Xtemp; Ysum += Ytemp;
      XYsum += (Xtemp * Ytemp); XXsum += (Xtemp * Xtemp);  }
  pointCount = last - first + 1 ;
  Xmean = Xsum/pointCount;  Ymean = Ysum/pointCount;
  *slopePtr = (XYsum - Xsum * Ymean)/(XXsum - Xsum * Xmean) ;
  *interceptPtr = Ymean - *slopePtr * Xmean ;
}
/* ################################################################## */
/* ################################################################## */
/*
MARK GREATEST INDEX WHERE COUNT > boxCountF[0]/cutOff_factor.

COUNTS AT LESSER INDEXES WILL NOT BE USED IN THE ESTIMATE OF FRACTAL
DIMENSION -- DISTORTION DUE TO SATURATION IS THE CONCERN.

NOTE THAT boxCountF[0] IS THE NUMBER OF BOXES OF SIZE 1 (THE SMALLEST SIZE)
THAT CONTAIN A POINT OF THE SET.  FOR ALL PRACTICAL PURPOSES, boxCountF[0]
WILL EQUAL THE NUMBER OF DISTINCT POINTS IN THE INPUT FILE, BECAUSE THESE
BOXES ARE REALLY SMALL COMPARED TO THE SIZE OF THE BIGGEST BOX (ABOUT 4
BILLION IF AN UNSIGNED LONG INT IS 32 BITS TO THE PLATFORM COMPILER.  THE
POINTS ARE SCALED BY THE PROGRAM SO THAT THE SET IS TOO "LARGE" TO FIT IN
THE NEXT SMALLEST BOX SIZE, SO THAT "1" IS THE SMALLEST DIFFERENCE IN
VALUE THAT CAN BE RESOLVED.) ONE BOX, IN EFFECT, COVERS ONLY A SINGLE POINT
OF THE INPUT SET BECAUSE THE PROGRAM CAN'T RESOLVE POINTS WITH A SMALLER
DIFFERENCE.

WE THINK IT WOULD BE A BAD IDEA TO USE dataLines/cutOff_factor AS THE LIMIT
BECAUSE IN CASES WHERE THERE WERE MANY DUPLICATE POINTS, WE WOULD SERIOUSLY
OVER-ESTIMATE THE NUMBER OF DISTINCT POINTS, AND THUS USE SATURATED DATA TO
BASE THE ESTIMATE OF FRACTAL DIMENSION UPON.  WHEN TESTING THE PROGRAM WITH
RANDOM DATA SAMPLED WITH REPLACEMENT, THIS COULD THROW THE RESULTS WAY OFF.
(THIS HAPPENED TO US, AND IT TOOK US A WHILE TO FIGURE OUT WHY.  AFTERWARDS,
WE STOPPED USING dataLines/cutOff_factor, AND CHANGED TO
boxCountF[0]/cutOff_factor.)
*/
void findMark(ulong *markPtr, ulong boxCountM[])
//   ulong *markPtr, boxCountM[numbits+1] ;
{
    int     i,  cutOff_factor=1;
    
   /* Calculate cutOff_factor = 2^(embed_dim) + 1 */
 for (i=1;i<=embed_dim;i++) cutOff_factor = cutOff_factor * 2;
 cutOff_factor++;

 *markPtr=0;
 for(i=0;i boxCountM[0]/cutOff_factor) *markPtr=i; }
}

/* ################################################################## */
/* ################################################################## */
void GetDims(double negLogBoxCountF[], double logSumSqrFreqF[],
	      double informationF[], ulong markF,
	      double *capDimPtr, double *infDimPtr, double *corrDimPtr)
//      ulong   markF ;
//      double  negLogBoxCountF[numbits+1], informationF[numbits+1],
//              logSumSqrFreqF[numbits+1],
//              *capDimPtr, *infDimPtr, *corrDimPtr;
{
    int     i;
    double  logEps[numbits+1], slope, intercept;
    
    /* GET LOG (BASE 2) OF THE DIAMETER OF THE I'TH SIZE OF BOX. */
  for(i=numbits; i>=0; i--)  logEps[i] = i;

/* fitLSqrLine (markF, numbits-4, logEps, negLogBoxCountF, &slope,&intercept);*/
  fitLSqrLine (markF, numbits-2, logEps, negLogBoxCountF, &slope, &intercept);
  *capDimPtr = slope ;
/*fitLSqrLine(markF, numbits-4, logEps, informationF, &slope, &intercept);*/
  fitLSqrLine(markF, numbits-2, logEps, informationF, &slope, &intercept);
  *infDimPtr = slope ;
/*fitLSqrLine(markF,numbits-4, logEps, logSumSqrFreqF, &slope, &intercept);*/
  fitLSqrLine(markF,numbits-2, logEps, logSumSqrFreqF, &slope, &intercept);
  *corrDimPtr = slope ;
}
/* ################################################################## */
/* ################################################################## */

void main(int argc, char *argv[])
//        int     argc;
//        char    *argv[];
{
  extern FILE  *infile;
  extern int   embed_dim,         /* How many coordinates points have. */
               *marker;           /*  To mark queues during radix sort */
  extern ulong **data,    /* initial array of data points */
               dataLines,        /* number of data points in input file */
              *diff_test;        /* XOR's of pairs of coordinates */
  extern long int  *next_after,    /*  array for pointers to data. */
                   avail;             /*  avail list pointer */
  extern QHeader Q[2];       /* array of two queues for the radix sort */


  int    i,j,m;
  ulong  maxdiam=0,              /* 2^numbits - 1: determines the scaling. */
         boxCount[numbits+1],  /* Array to keep track of the box
                                         counts of each size.  Each index
                                         is the base-2 log of the
                                         corresponding box size. */
         pointCount[numbits+1],  /*  Array to keep track of how many points
                                     of the data set are contained in each
                                     box */
	 mark;                  /* Marks smallest usable box count */
  double    max, min,   /* Maximum and minimum values in input file-- we
                           need to know these in order to scale the
                           input points. */
            buf,        /* A buffer for inputs to stay in until they are
                           scaled and converted to integers.  */
	    negLogBoxCount[numbits+1],
            logSumSqrFreq[numbits+1],
            information[numbits+1],
	    capDim, infDim, corrDim;
  long int  pointer;    /* holds indexes used as pointers to records */

printf("\n\n") ;
printf("******************************************************************\n");
printf("  FRACTAL DIMENSION REPORT -- by fd software (DiFalco/Sarraille)\n");
printf("******************************************************************\n");
printf("\n") ;

     /*  IDENTIFY THE INPUT FILE */
  printf ("\nReporting on file named: %s.\n\n", argv[1]) ;

     /* FIND OUT HOW MANY COORDINATES POINTS HAVE -- 1?, 2?, 3?, MORE? */
  printf("Getting the embedding dimension ...\n");
  embed_dim = get_e_dim(argv[1]);
  printf("Embedding Dimension taken to be: %d\n\n",embed_dim);

/*
   find max and min in input file.  This function should open the input file
   for reading, scan it to get the maximum and minimum data values it
   contains, and then seek to the beginning of the input file so that the
   code below can proceed to read the data sequentially from the beginning.
   We are having trouble getting that to work, so we have temporarily taken
   the expedient of closing the file in the function max_min, and opening it
   again below.
*/
  printf("Finding max and min values in data ...");
  fflush(stdout);
  max_min (argv[1], &max, &min) ;
  printf ("\n\nMinimum value in input file is: %lf ...\n", min) ;
  printf ("Maximum value in input file is: %lf ...\n\n", max) ;
     /*
        set maxdiam=2^numbits-1 -- this is the largest value that an
        unsigned integer can have with the host compiler.
        (I checked the value obtained in this manner, and it
         is correct -- j.s.)
     */
  if (debugging || checking) printf("numbits is: %d\n",numbits);
  printf("%d different cell sizes will be used ...\n",numbits);  
  for (m=0;m= numbits-2)
  {
    printf("\nInsufficient Data.  More points are needed.\n") ;
    printf("Cannot assign a Fractal Dimension to this set.\n");
    printf("Examine the data above to make your own conjecture.\n");
    exit(-1) ;
  }

  printf ("\n\n\n%d is the smallest cell size used in ", mark);
  printf ("the overall dimension estimates\n");
  printf ("below.  The largest cell size is ");
  printf ("%d.  Data above corresponding to\n", numbits-2);
  printf ("this range is between rows of asterisks.\n\n");

  GetDims(negLogBoxCount, logSumSqrFreq, information, mark,
	  &capDim, &infDim, &corrDim) ;

  printf("\n\nLeast-Square Estimates based on Indicated Cell Range:\n\n");
  printf("Fractal Dimension  (Capacity)   =  %.5f\n", capDim );
  printf("Fractal Dimension (Information) =  %.5f\n", infDim );
  printf("Fractal Dimension (Correlation) =  %.5f\n", corrDim);
  printf("\n\n****************************************\n");    
}