/* file sorts, various sort routines for different types of data */
 /* taken from Press etal */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "ana_structures.h"
 /*------------------------------------------------------------------------- */
#define ALN2I 1.442695022
#define TINY 1.0e-5
 /*------------------------------------------------------------------------- */
void shell_f(n,arr)
 float arr[];
 int n;
{
 int nn,m,j,i,lognb2;
 float t;

 lognb2=(log((double) n)*ALN2I+TINY);
 m=n;
 for (nn=1;nn<=lognb2;nn++) {
  m >>= 1;
  for (j=m;j<n;j++) {
   i=j-m;
   t=arr[j];
   while (i >= 0 && arr[i] > t) {
    arr[i+m]=arr[i];
    i -= m;
   }
   arr[i+m]=t;
  }
 }
}
 /*------------------------------------------------------------------------- */
void shell_b(n,arr)
 byte arr[];
 int n;
{
 int nn,m,j,i,lognb2;
 byte t;

 lognb2=(log((double) n)*ALN2I+TINY);
 m=n;
 for (nn=1;nn<=lognb2;nn++) {
  m >>= 1;
  for (j=m;j<n;j++) {
   i=j-m;
   t=arr[j];
   while (i >= 0 && arr[i] > t) {
    arr[i+m]=arr[i];
    i -= m;
   }
   arr[i+m]=t;
  }
 }
}
 /*------------------------------------------------------------------------- */
void shell_l(n,arr)
 int arr[];
 int n;
{
 int nn,m,j,i,lognb2;
 int t;

 lognb2=(log((double) n)*ALN2I+TINY);
 m=n;
 for (nn=1;nn<=lognb2;nn++) {
  m >>= 1;
  for (j=m;j<n;j++) {
   i=j-m;
   t=arr[j];
   while (i >= 0 && arr[i] > t) {
    arr[i+m]=arr[i];
    i -= m;
   }
   arr[i+m]=t;
  }
 }
}
 /*------------------------------------------------------------------------- */
void shell_w(n,arr)
 short arr[];
 int n;
{
 int nn,m,j,i,lognb2;
 short t;

 lognb2=(log((double) n)*ALN2I+TINY);
 m=n;
 for (nn=1;nn<=lognb2;nn++) {
  m >>= 1;
  for (j=m;j<n;j++) {
   i=j-m;
   t=arr[j];
   while (i >= 0 && arr[i] > t) {
    arr[i+m]=arr[i];
    i -= m;
   }
   arr[i+m]=t;
  }
 }
}
 /*------------------------------------------------------------------------- */
void shell_d(n,arr)
 double arr[];
 int n;
{
 int nn,m,j,i,lognb2;
 double t;

 lognb2=(log((double) n)*ALN2I+TINY);
 m=n;
 for (nn=1;nn<=lognb2;nn++) {
  m >>= 1;
  for (j=m;j<n;j++) {
   i=j-m;
   t=arr[j];
   while (i >= 0 && arr[i] > t) {
    arr[i+m]=arr[i];
    i -= m;
   }
   arr[i+m]=t;
  }
 }
}
 /*------------------------------------------------------------------------- */
void sort_f(n,ra)
 int n;
 float ra[];
{
 int l,j,ir,i;
 float rra;

 l = (n/2);
 ir = n-1;
 for (;;) {
  if (l > 0)
   rra = ra[--l];
  else {
   rra = ra[ir];
   ra[ir] = ra[0];
   if (--ir == 0) {
    ra[0] = rra;
    return;
   }
  }
  i = l;
  j = l + l + 1;
  while (j <= ir) {
   if (j < ir && ra[j] < ra[j+1]) j++;
   if (rra < ra[j]) {
    ra[i] = ra[j];
    j += (i=j) + 1;
   }
   else j = ir + 1;
  }
  ra[i] = rra;
 }
}
 /*------------------------------------------------------------------------- */
void sort_b(n,ra)
 int n;
 byte ra[];
{
 int l,j,ir,i;
 byte rra;

 l = (n/2);
 ir = n-1;
 for (;;) {
  if (l > 0)
   rra = ra[--l];
  else {
   rra = ra[ir];
   ra[ir] = ra[0];
   if (--ir == 0) {
    ra[0] = rra;
    return;
   }
  }
  i = l;
  j = l + l + 1;
  while (j <= ir) {
   if (j < ir && ra[j] < ra[j+1]) j++;
   if (rra < ra[j]) {
    ra[i] = ra[j];
    j += (i=j) + 1;
   }
   else j = ir + 1;
  }
  ra[i] = rra;
 }
}
 /*------------------------------------------------------------------------- */
void sort_w(n,ra)
 int n;
 short ra[];
{
 int l,j,ir,i;
 short rra;

 l = (n/2);
 ir = n-1;
 for (;;) {
  if (l > 0)
   rra = ra[--l];
  else {
   rra = ra[ir];
   ra[ir] = ra[0];
   if (--ir == 0) {
    ra[0] = rra;
    return;
   }
  }
  i = l;
  j = l + l + 1;
  while (j <= ir) {
   if (j < ir && ra[j] < ra[j+1]) j++;
   if (rra < ra[j]) {
    ra[i] = ra[j];
    j += (i=j) + 1;
   }
   else j = ir + 1;
  }
  ra[i] = rra;
 }
}
 /*------------------------------------------------------------------------- */
void sort_l(n,ra)
 int n;
 int ra[];
{
 int l,j,ir,i;
 int rra;

 l = (n/2);
 ir = n-1;
 for (;;) {
  if (l > 0)
   rra = ra[--l];
  else {
   rra = ra[ir];
   ra[ir] = ra[0];
   if (--ir == 0) {
    ra[0] = rra;
    return;
   }
  }
  i = l;
  j = l + l + 1;
  while (j <= ir) {
   if (j < ir && ra[j] < ra[j+1]) j++;
   if (rra < ra[j]) {
    ra[i] = ra[j];
    j += (i=j) + 1;
   }
   else j = ir + 1;
  }
  ra[i] = rra;
 }
}
 /*------------------------------------------------------------------------- */
void sort_d(n,ra)
 int n;
 double ra[];
{
 int l,j,ir,i;
 double rra;

 l = (n/2);
 ir = n-1;
 for (;;) {
  if (l > 0)
   rra = ra[--l];
  else {
   rra = ra[ir];
   ra[ir] = ra[0];
   if (--ir == 0) {
    ra[0] = rra;
    return;
   }
  }
  i = l;
  j = l + l + 1;
  while (j <= ir) {
   if (j < ir && ra[j] < ra[j+1]) j++;
   if (rra < ra[j]) {
    ra[i] = ra[j];
    j += (i=j) + 1;
   }
   else j = ir + 1;
  }
  ra[i] = rra;
 }
}
 /*------------------------------------------------------------------------- */
void sort_s(n,ra)
 /* sort a herd of strings */
 /* just moves the pointers around */
 int n;
 char * ra[];
{
 int l,j,ir,i;
 char * rra;

 l = (n/2);
 ir = n-1;
 for (;;) {
  if (l > 0)
   rra = ra[--l];
  else {
   rra = ra[ir];
   ra[ir] = ra[0];
   if (--ir == 0) {
    ra[0] = rra;
    return;
   }
  }
  i = l;
  j = l + l + 1;
  while (j <= ir) {
   if (j < ir && strcmp(ra[j],ra[j+1]) < 0 ) j++;
   if (strcmp(rra, ra[j]) < 0) {
    ra[i] = ra[j];
    j += (i=j) + 1;
   }
   else j = ir + 1;
  }
  ra[i] = rra;
 }
}
 /*------------------------------------------------------------------------- */
void indexx_d(n,ra,indx)
 int	n;
 double	ra[];
 int	indx[];
 {
  int l,j,ir,i,indxt;
  double q;
 
  for (i=0;i<n;i++) indx[i] = i;
  l = (n/2);
  ir = n-1;
  for (;;) {
   if (l > 0)
    q=ra[(indxt=indx[--l])];
   else {
    q=ra[(indxt=indx[ir])];
    indx[ir]=indx[0];
    if (--ir == 0) {indx[0]=indxt;return; }
   }
   i = l;
   j = l + l + 1;
   while (j <= ir) {
    if (j < ir && ra[indx[j]] < ra[indx[j+1]]) j++;
    if (q < ra[indx[j]]) {
     indx[i] = indx[j];
     j += (i=j) + 1;
    }
    else j = ir + 1;
   }
   indx[i] =indxt;
  }
 }
 /*------------------------------------------------------------------------- */
void indexx_b(n,ra,indx)
 int	n;
 byte	ra[];
 int	indx[];
 {
  int l,j,ir,i,indxt;
  byte q;
 
  for (i=0;i<n;i++) indx[i] = i;
  l = (n/2);
  ir = n-1;
  for (;;) {
   if (l > 0)
    q=ra[(indxt=indx[--l])];
   else {
    q=ra[(indxt=indx[ir])];
    indx[ir]=indx[0];
    if (--ir == 0) {indx[0]=indxt;return; }
   }
   i = l;
   j = l + l + 1;
   while (j <= ir) {
    if (j < ir && ra[indx[j]] < ra[indx[j+1]]) j++;
    if (q < ra[indx[j]]) {
     indx[i] = indx[j];
     j += (i=j) + 1;
    }
    else j = ir + 1;
   }
   indx[i] =indxt;
  }
 }
 /*------------------------------------------------------------------------- */
void indexx_w(n,ra,indx)
 int	n;
 short	ra[];
 int	indx[];
 {
  int l,j,ir,i,indxt;
  short q;
 
  for (i=0;i<n;i++) indx[i] = i;
  l = (n/2);
  ir = n-1;
  for (;;) {
   if (l > 0)
    q=ra[(indxt=indx[--l])];
   else {
    q=ra[(indxt=indx[ir])];
    indx[ir]=indx[0];
    if (--ir == 0) {indx[0]=indxt;return; }
   }
   i = l;
   j = l + l + 1;
   while (j <= ir) {
    if (j < ir && ra[indx[j]] < ra[indx[j+1]]) j++;
    if (q < ra[indx[j]]) {
     indx[i] = indx[j];
     j += (i=j) + 1;
    }
    else j = ir + 1;
   }
   indx[i] =indxt;
  }
 }
 /*------------------------------------------------------------------------- */
void indexx_l(n,ra,indx)
 int	n;
 int	ra[];
 int	indx[];
 {
  int l,j,ir,i,indxt;
  int q;
 
  for (i=0;i<n;i++) indx[i] = i;
  l = (n/2);
  ir = n-1;
  for (;;) {
   if (l > 0)
    q=ra[(indxt=indx[--l])];
   else {
    q=ra[(indxt=indx[ir])];
    indx[ir]=indx[0];
    if (--ir == 0) {indx[0]=indxt;return; }
   }
   i = l;
   j = l + l + 1;
   while (j <= ir) {
    if (j < ir && ra[indx[j]] < ra[indx[j+1]]) j++;
    if (q < ra[indx[j]]) {
     indx[i] = indx[j];
     j += (i=j) + 1;
    }
    else j = ir + 1;
   }
   indx[i] =indxt;
  }
 }
 /*------------------------------------------------------------------------- */
void indexx_f(n,ra,indx)
 int	n;
 float	ra[];
 int	indx[];
 {
  int l,j,ir,i,indxt;
  float q;
 
  for (i=0;i<n;i++) indx[i] = i;
  l = (n/2);
  ir = n-1;
  for (;;) {
   if (l > 0)
    q=ra[(indxt=indx[--l])];
   else {
    q=ra[(indxt=indx[ir])];
    indx[ir]=indx[0];
    if (--ir == 0) {indx[0]=indxt;return; }
   }
   i = l;
   j = l + l + 1;
   while (j <= ir) {
    if (j < ir && ra[indx[j]] < ra[indx[j+1]]) j++;
    if (q < ra[indx[j]]) {
     indx[i] = indx[j];
     j += (i=j) + 1;
    }
    else j = ir + 1;
   }
   indx[i] =indxt;
  }
 }
 /*------------------------------------------------------------------------- */
void indexx_s(n,ra,indx)
 /* sort a herd of strings */
 int	n;
 char ** ra;
 int	indx[];
 {
  int l,j,ir,i,indxt;
  char * q;
 
  for (i=0;i<n;i++) indx[i] = i;
  l = (n/2);
  ir = n-1;
  while (1) {
   if (l > 0)
    q=ra[(indxt=indx[--l])];
   else {
    q=ra[(indxt=indx[ir])];
    indx[ir]=indx[0];
    if (--ir == 0) {indx[0]=indxt;return; }
   }
   i = l;
   j = l + l + 1;
   while (j <= ir) {
    if (j < ir && strcmp(ra[indx[j]], ra[indx[j+1]]) < 0) j++;
    if (strcmp(q, ra[indx[j]]) < 0 ) {
     indx[i] = indx[j];
     j += (i=j) + 1;
    }
    else j = ir + 1;
   }
   indx[i] =indxt;
  }
 }
 /*------------------------------------------------------------------------- */
#undef ALN2I
#undef TINY