/**********************************************************/
/*                                                        */
/* Program: qsortd.c                                      */
/*                                                        */
/* Test program for non recursive quick sort.             */
/*                                                        */
/* Array type: double precision                           */
/*                                                        */
/* This is the fastest sorting algorithm known to the     */
/* author of this test program.                           */
/* The person that developed the algorithm is unknown.    */
/*                                                        */
/* Author: Stefan Spaennare, Lund, Sweden                 */
/*                                                        */
/* E-mail: stefan@spaennare.se                            */
/*                                                        */
/* Latest update: 2003-02-19                              */
/*                                                        */
/*                                                        *************/
/* Benchmark                                                          */
/* =========                                                          */
/*                                                                    */
/* The test was run on a Intel Celeron (Tualatin core) computer       */
/* with 256 kbyte cache memory and 512 Mbyte primary memory.          */
/* The operating system was Red Hat 7.2 Linux. Compiler tested        */
/* was GNU GCC 2.96 (gcc) for Linux. Under Windows XP the compiler    */
/* GNU DJGPP 3.0.4 (gcc, djgpp) for DOS was used.                     */
/*                                                                    */
/* Numbers tested:                                                    */
/* ---------------                                                    */
/*                                                                    */
/* seed = 3                                                           */
/*                                                                    */
/* na =    1000000 =     10^6                                         */
/*                                                                    */
/* nb =   10000000 =  10*10^6                                         */
/*                                                                    */
/* Results:                                                           */
/* --------                                                           */
/*                                                                    */
/* n:  Computer:         CPU-time (s):  Compiler options:             */
/* ================================================================   */
/*                                                                    */
/* na  Celeron 1200 MHz       0.58      gcc -O3                       */
/* na  Celeron 1200 MHz       0.55      gcc -O2 (djgpp)               */
/*                                                                    */
/* nb  Celeron 1200 MHz       7.61      gcc -O3                       */
/* nb  Celeron 1200 MHz       7.20      gcc -O2 (djgpp)               */
/*                                                                    */
/*                                                                    */
/*                                                        *************/
/* Old benchmark                                          */
/* =============                                          */
/*                                                        */
/* seed=3; n=1000000                                      */
/*                                                        */
/* Computer:           CPU-time (s):      Compilation:    */
/*                                                        */
/* Pentium 180 MHz         5.02             gcc -O3       */
/*                                                        */
/**********************************************************/

/********************************************************************************/
/*                                                                              */
/* Notice                                                                       */
/* ======                                                                       */
/*                                                                              */
/* I make no warranties that this program is (1) free of errors, (2) consistent */
/* with any standard merchantability, or (3) meeting the requirements of a      */
/* particular application. This software shall not, partly or as a whole,       */
/* participate in a process, whose outcome can result in injury to a person or  */
/* loss of property. It is solely designed for analytical work. Permission to   */
/* use, copy and distribute is hereby granted without fee, providing that the   */
/* header above including this notice appears in all copies.                    */
/*                                                                              */
/*                                                            Stefan Spaennare  */
/*                                                                              */
/********************************************************************************/


#include <math.h>
#include <stdio.h>
#include <unistd.h>
#include <time.h>
#include <stdlib.h>

double timetic=(double)(CLOCKS_PER_SEC);

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran2(idum)
long *idum;
{
   int j;
   long k;
   static long idum2=123456789;
   static long iy=0; 
   static long iv[NTAB];
   double temp;

   if (*idum <= 0) {
      if (-(*idum) < 1) {
         *idum=1;
      }
      else {
         *idum=-(*idum);
      } /* if */
      idum2=(*idum);
      for (j=NTAB+7; j>=0; j--) {
         k=(*idum)/IQ1;
	 *idum=IA1*(*idum-k*IQ1)-k*IR1;
	 if (*idum < 0) {
	    *idum += IM1;
	 } /* if */
	 if (j < NTAB) {
	    iv[j]=*idum;
	 } /* if */
      } /* for j */
      iy=iv[0];
   } /* if */
   k=(*idum)/IQ1;
   *idum=IA1*(*idum-k*IQ1)-k*IR1;
   if (*idum < 0) {
      *idum += IM1;
   } /* if */
   k=idum2/IQ2;
   idum2=IA2*(idum2-k*IQ2)-k*IR2;
   if (idum2 < 0) {
      idum2 += IM2;
   } /* if */
   j=iy/NDIV;
   iy=iv[j]-idum2;
   iv[j] = *idum;
   if (iy < 1) {
      iy += IMM1;
   } /* if */
   temp=AM*(double)(iy);
   if (temp > RNMX) {
      return(RNMX);
   }
   else {
      return(temp);
   } /* if */

} /* ran2 */

void qqsort(f,sr,sl,n)
double *f;
int *sr;
int *sl;
int n;
{
  int l,r,s,i,j;

  double m,temp;

  s=1;
  sl[s]=1;
  sr[s]=n;

  do {
     l=sl[s];
     r=sr[s];
     s--;
    do { 
     i=l;
     j=r;
     m=f[(l+r)/2];
     do {
        while (f[i] < m) {
           i++;
        } /* while */
        while (f[j] > m) {
           j--;
        } /* while */
        if (i <= j) {
           temp=f[i];
           f[i]=f[j];
           f[j]=temp;
           i++;
           j--;
        } /* if */
     } while (i <= j);
     if (i < r) {
        s++;
        sl[s]=i;
        sr[s]=r;
     } /* if */
     r=j;
    } while (l < r);
  } while (s > 0);
   
} /* qqsort */



int main(argc,argv) 
int argc;
char *argv[];

{
  
   int i,n,seed;
   long idum;

   double *v;

   int *sr,*sl;

   clock_t time1,timediff;


   if (argc != 3) {
      printf("Usage: %s seedi ni\n",argv[0]);
      exit(0);
   } /* if */


seed=atoi(argv[1]);
n=atoi(argv[2]);

v = (double *)calloc(n+1,sizeof(double));
sr = (int *)calloc((int)(log((double)(n))/log(2.0))+10,sizeof(int));
sl = (int *)calloc((int)(log((double)(n))/log(2.0))+10,sizeof(int));


idum=-abs(seed)-1;

for (i=1; i<=n; i++) {
   v[i]=200.0*(ran2(&idum)-0.5);
} /* for i */


if (n <=20) {
   for (i=1; i<=n; i++) {
      printf("%14.5f\n",v[i]);
   } /* for i */
} /* if */


time1=clock();

qqsort(v,sr,sl,n);


timediff=clock()-time1;
printf("\n%10.2f\n\n",(double)(timediff)/timetic);

if (n <=20) {
   for (i=1; i<=n; i++) {
      printf("%14.5f\n",v[i]);
   } /* for i */
} /* if */

} /* End  */