XEphem/libip/median.c

117 lines
2.3 KiB
C

/* Find median. Algorithm from N. Wirth's book, implementation by N. Devillard.
* adapted for XEphem by Elwood Downey
* #define TMAIN for test program, see end.
*/
#include "ip.h"
/*---------------------------------------------------------------------------
Function : kth_smallest()
In : array of elements, # of elements in the array, rank k
Out : one element
Job : find the kth smallest element in the array
Reference:
Author: Wirth, Niklaus
Title: Algorithms + data structures = programs
Publisher: Englewood Cliffs: Prentice-Hall, 1976
Physical description: 366 p.
Series: Prentice-Hall Series in Automatic Computation
---------------------------------------------------------------------------*/
static double
dkth_smallest(double a[], int n, int k)
{
int i,j,l,m ;
double t, x ;
l=0 ; m=n-1 ;
while (l<m) {
x=a[k] ;
i=l ;
j=m ;
do {
while (a[i]<x) i++ ;
while (x<a[j]) j-- ;
if (i<=j) {
t=a[i]; a[i]=a[j]; a[j]=t;
i++ ; j-- ;
}
} while (i<=j) ;
if (j<k) l=i ;
if (k<i) m=j ;
}
return a[k] ;
}
/* find median of n double array.
* N.B. array is rearranged IN PLACE.
*/
double
dmedian(double a[], int n)
{
return (dkth_smallest (a, n, n/2));
}
static CamPix
ckth_smallest(CamPix a[], int n, int k)
{
int i,j,l,m ;
CamPix t, x ;
l=0 ; m=n-1 ;
while (l<m) {
x=a[k] ;
i=l ;
j=m ;
do {
while (a[i]<x) i++ ;
while (x<a[j]) j-- ;
if (i<=j) {
t=a[i]; a[i]=a[j]; a[j]=t;
i++ ; j-- ;
}
} while (i<=j) ;
if (j<k) l=i ;
if (k<i) m=j ;
}
return a[k] ;
}
/* find median of n CamPix array.
* N.B. array is rearranged IN PLACE.
*/
CamPix
cmedian(CamPix a[], int n)
{
return (ckth_smallest (a, n, n/2));
}
#ifdef TMAIN
#include <stdio.h>
#include <stdlib.h>
int
main (int ac, char *av[])
{
double *a = malloc (1);
char buf[128];
int n;
for (n = 0; fgets (buf, sizeof(buf), stdin); n++)
(a = realloc (a, (n+1)*sizeof(double)))[n] = atof(buf);
printf ("%g\n", dmedian(a, n));
return(0);
}
#endif