mathstat.c

国外网站上的一些精典的C程序

/*
**  MATHSTAT.C - Statistical analysis in C and C++
**
**  Public domain by Bob Stout
**
**  Provided functions:
**  1.  Open/Initialize
**  2.  Close
**  3.  Flush/Reinitialize
**  4.  Add datum
**  5.  Delete datum
**  6.  "Olympic" filter - discards high and low data.
**  7.  "Olympiad" filter - an auto-windowing median filter, it iteratively discards
**      high and low data until a target coefficient of variation is achieved.
**  8.  Median Filter - discards data outside of a window around the median datum.
**  9.  Minumum datum
**  10. Maximum datum
**  11. Minimum error - percent deivation from arithmetic mean
**  12. Maximum error - percent deivation from arithmetic mean
**  13. Arithmetic mean
**  14. Geometric mean
**  15. Harmonic mean
**  16. Standard deviation (population)
**  17. Standard deviation (sample)
**  18. Variance (population)
**  19. Variance (sample)
**  20. Coefficient of variation (population)
**  21. Coefficient of variation (sample)
**
**  Uses SNIPPETS files:
**    MATHSTAT.H
**    SNIPTYPE.H
**    CAST.H
**    BINSRCH.H
**    BINSRCH.C
*/

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <float.h>
#include <assert.h>
#include "mathstat.h"
#include "cast.h"
#include "binsrch.h"

/*
**  We need a local compiler-independent NaN
*/

VAR64_ stat_nan_;
static Boolean_T nan_initialized = False_;

/*
**  For functions requiring min and max data, all data are stored in an array that
**  is sized upon object creation. Data are always added to the array in sorted order.
*/


/*
**  Initialize a stat structure
**
**  Parameters: 1 - Pointer to Stat_T structure type
**              2 - Data array depth (use 0 for no data array and no filtering)
**
**  Returns: Success_ or Error_
**
**  Side effects: Allocates data array
*/

#if !(__cplusplus)

Boolean_T stat_init(Stat_T *ptr, size_t qdepth)
{
      if (!nan_initialized)
      {
            stat_nan_.v32[0].l = -1L;
            stat_nan_.v32[1].l = -1L;
            nan_initialized = True_;
      }
      ptr->count   = 0;
      ptr->len     = qdepth;
      ptr->total   = 0.0;
      ptr->total2  = 0.0;
      ptr->recip   = 0.0;
      ptr->product = 1.0;
      if (qdepth)
      {
            if (NULL == (ptr->filt = calloc(qdepth, sizeof(double))))
                  return Error_;
      }
      else
      {
            ptr->filt = NULL;
            ptr->min = FLT_MAX;
            ptr->max = FLT_MIN;
      }
      return Success_;
}

#else /* C++ */

inline Boolean_T Stat::Stat(size_t qdepth)
{
      if (!nan_initialized)
      {
            stat_nan_.v32[0].l = -1L;
            stat_nan_.v32[1].l = -1L;
            nan_initialized = True_;
      }
      count   = 0;
      len     = qdepth;
      total   = 0.0;
      total2  = 0.0;
      recip   = 0.0;
      product = 1.0;
      if (qdepth)
      {
            if (NULL == (filt = calloc(qdepth, sizeof(double))))
                  return Error_;
      }
      else
      {
            filt = NULL;
            min = FLT_MAX;
            max = FLT_MIN;
      }
      return Success_;
}

#endif

/*
**  Flush/reinitialize a stat structure
**
**  Parameters: 1 - Pointer to Stat_T structure type
**
**  Returns: Nothing
*/

#if !(__cplusplus)

void stat_flush(Stat_T *ptr)
{
      ptr->count   = 0;
      ptr->total   = 0.0;
      ptr->total2  = 0.0;
      ptr->recip   = 0.0;
      ptr->product = 1.0;
      if (0 == ptr->len)
      {
            ptr->min = FLT_MAX;
            ptr->max = FLT_MIN;
      }
}

#else /* C++ */

inline void Stat::Flush()
{
      Boolean_T retval = Success_;

      count   = 0;
      total   = 0.0;
      total2  = 0.0;
      recip   = 0.0;
      product = 1.0;
      if (0 == ptr->len)
      {
            min = FLT_MAX;
            max = FLT_MIN;
      }
}

#endif

/*
**  Closee a stat structure
**
**  Parameters: 1 - Pointer to Stat_T structure type
**
**  Returns: Nothing
**
**  Side effects: Frees the data array
*/

#if !(__cplusplus)

void stat_done(Stat_T *ptr)
{
      if (ptr->filt)
            FREE(ptr->filt);
}

#else /* C++ */

inline Boolean_T ~Stat::Stat()
{
      if (ptr->filt)
            FREE(filt);
}

#endif

/*
**  Return the number of data points
**
**  Parameters: 1 - Pointer to Stat_T structure type
**
**  Returns: Number of data points
*/

#if !(__cplusplus)

size_t stat_count(Stat_T *ptr)
{
      return ptr->count;
}

#else /* C++ */

inline size_t Stat::Count()
{
      return count;
}

#endif

/*
**  Add a value for statistical analysis
**
**  Parameters: 1 - Datum to add
**              2 - Pointer to Stat_T structure type
**
**  Returns: Error_ if data array is full, else Success_
**
**  Side effects: Inserts the new datum into the data array in sorted order
*/

#if !(__cplusplus)

Boolean_T stat_add(double datum, Stat_T *ptr)
{
      int high, low, idx;

      if ((0 == ptr->len) || (ptr->count < ptr->len))
      {
            ptr->total   += datum;
            ptr->total2  += datum * datum;
            ptr->recip   += 1.0 / datum;
            ptr->product *= datum;

        /*
        **  If a data array is being used, maintain it
        */

            if (ptr->filt)
            {
                  /*
                  ** Place the data in sorted order
                  **
                  ** If the array is empty or the datum is a new max, simply append it
                  */
                  
                  if ((0 == ptr->count) || (ptr->filt[ptr->count - 1] <= datum))
                        ptr->filt[ptr->count++] = datum;
                  
                  /*
                  ** If the datum is a new min, shift the array up and insert it at the beginning.
                  */
                  
                  else if (ptr->filt[0] >= datum)
                  {
                        memmove(&ptr->filt[1], &ptr->filt[0], (ptr->count++) * sizeof(datum));
                        ptr->filt[0] = datum;
                  }
                  
                  /*
                  ** Use a binary search to insert the datum into the array
                  */
                  
                  else
                  {
                        for (low = 0, high = ptr->count - 1;  high - low > 1; )
                        {
                              idx = (high+low) / 2;
                              if (datum == ptr->filt[idx])
                                    break;
                              else if (datum < ptr->filt[idx])
                                    high = idx;
                              else  low  = idx;
                        }
                        if (datum != ptr->filt[idx])
                              idx = high;
                        memmove(&ptr->filt[idx+1], &ptr->filt[idx], (ptr->count++ - idx) * sizeof(datum));
                        ptr->filt[idx] = datum;
                  }
                  return Success_;
            }
            else
            {
                  ptr->count++;
                  if (datum > ptr->max)
                        ptr->max = datum;
                  if (datum < ptr->min)
                        ptr->min = datum;
            }
      }
      else  return Error_;
}

#else /* C++ */

Boolean_T Stat::Add(double datum)
{
      int high, low, idx;

        if ((0 == ptr->len) || (ptr->count < ptr->len))
      {
            total   += datum;
            total2  += datum * datum;
            recip   += 1.0 / datum;
            product *= datum;

            /*
            **  If a data array is being used, maintain it
            */

            if (filt)
            {
                  /*
                  ** Place the data in sorted order
                  **
                  ** If the array is empty or the datum is a new max, simply append it
                  */
                  
                  if ((0 == count) || (filt[count - 1] <= datum))
                        filt[count++] = datum;
                  
                  /*
                  ** If the datum is a new min, shift the array up and insert it at the beginning.
                  */
                  
                  else if (filt[0] >= datum)
                  {
                        memmove(&filt[1], &filt[0], (count++) * sizeof(datum));
                        filt[0] = datum;
                  }
                  
                  /*
                  ** Use a binary search to insert the datum into the array
                  */

                  else
                  {
                        for (low = 0, high = count - 1;  high - low > 1; )
                        {
                              idx = (high+low) / 2;
                              if (datum == filt[idx])
                                    break;
                              else if (datum < filt[idx])
                                    high = idx;
                              else  low  = idx;
                        }
                        if (datum != filt[idx])
                              idx = high;
                        memmove(&filt[idx+1], &filt[idx], (count++ - idx) * sizeof(datum));
                        filt[idx] = datum;
                  }
                  return Success_;
            }
            else
            {
                  count++;
                  if (datum > max)
                        max = datum;
                  if (datum < min)
                        min = datum;
            }
      }
      else  return Error_;
}

#endif

/*
**  Delete a value from a statistical analysis
**
**  Parameters: 1 - Datum to add
**              2 - Pointer to Stat_T structure type
**
**  Returns: Error_ if data array is empty or doesn't contain the datum, else Success_
**
**  Side effects: Removes the datum from the data array, leaving it in sorted order
*/

#if !(__cplusplus)

Boolean_T stat_delete(double datum, Stat_T *ptr)
{
      if (ptr->count)
      {
            int idx;

            ptr->total   -= datum;
            ptr->total2  -= datum * datum;
            ptr->recip   -= 1.0 / datum;
            ptr->product /= datum;

            /*
            **  If a data array is being used, maintain it
            */

            if (ptr->filt)
            {
                  /*
                  **  Remove low data by physically shifting the array. This is laborious,
                  **  but no more so than implementing a decent sort on a circular buffer.
                  */
                  
                  if (datum == ptr->filt[0])
                        memmove(&ptr->filt[0], &ptr->filt[1], (--ptr->count) * sizeof(double));
                  
                  /*
                  **  Remove high data. Actually, if the datum being removed was at the top
                  **  of the array, it may be effectively wiped out by decremented the count,
                  **  so the only reamining task to locate data within the array and delete it.
                  **  Again, this is done by physically shifting data in the array.
                  */
                  
                  else if (datum != ptr->filt[ptr->count-1])
                  {
                        if (-1 == (idx = BinSearchD(datum, &ptr->filt[0], ptr->count)))
                              return Error_;
                        memmove(&ptr->filt[idx], &ptr->filt[idx+1], (--ptr->count - idx) * sizeof(double));
                  }
                  else  --(ptr->count);
                  
                  return Success_;
            }
      }
      return Error_;
}