Shewchuk

/*
 * Copyright 2013 Florian Philipp
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0

 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <math.h>
/* using fabs */
#include <stdlib.h>
/* using malloc, realloc, free, NULL */
#include <stdio.h>
/* using printf, fopen, fwrite, fclose */


/**
 * Shewchuk summation algorithm.
 *
 * Derived from msum here:
 * http://code.activestate.com/recipes/393090/
 */
struct Accumulator
{
  /**
   * Number of currently used partial sums
   */
  size_t partials_n;
  /**
   * Capacity for partial sums

   * partials_cap == 1 indicates use of partials.scalar
   */
  size_t partials_cap;
  /**
   * Partial sums ordered by value
   *
   * To avoid dynamic allocation in the common case that only one sum is kept,
   * the scalar union member is used.
   * The additional if/else comparisons have negligible effect on the run time
   */
  union { double* array; double scalar; } partials;
};


/**
 * Initializes an Accumulator
 *
 * An Accumulator may only be initialized once to avoid memory leaks
 * unless it is destroyed between initializations.
 * See also: accumulator_reset
 */
void accumulator_init(struct Accumulator* self)
{
  self->partials_n = 0;
  self->partials_cap = 1;
  self->partials.scalar = 0.;
}

/**
 * Releases all resources associated with this Accumulator
 *
 * May only be called once to avoid segfaults.
 * A destroyed Accumulator may only be used after calling accumulator_init.
 * As a rule of thumb, every accumulator_init should be paired with an
 * accumulator_destroy
 */
void accumulator_destroy(struct Accumulator* self)
{
  if(self->partials_cap > 1)
    free(self->partials.array);
}

/**
 * Re-initializes an Accumulator in order to use it for another summation
 */
void accumulator_reset(struct Accumulator* self)
{
  self->partials_n = 0;
  if(self->partials_cap == 1)
    self->partials.scalar = 0.;
}

/**
 * Private method of Accumulator. Updates existing partial sums
 *
 * \param self an initialized Accumulator
 * \param x a value to be added
 * \param tail output argument. Returns the index of the last used partial sum.
 * Is less or equal to self->partials_n
 *
 * \return the last partial sum. Has to be stored in self->partials[tail]
 */
static double _accumulator_update_partials(struct Accumulator* self,
                       double x,
                       size_t* tail)
{
  *tail = 0;
  double* partials;
  if(self->partials_cap > 1)
    partials = self->partials.array;
  else
    partials = &self->partials.scalar;
  size_t partial_i;
  for(partial_i = 0; partial_i < self->partials_n; ++partial_i) {
    double y = partials[partial_i];
    if(fabs(x) < fabs(y)) {
      double tmp = y;
      y = x;
      x = tmp;
    }
    double hi = x + y;
    double lo = y - (hi - x);
    if(lo != 0.) {
      partials[*tail] = lo;
      *tail += 1;
    }
    x = hi;
  }
  return x;
}

/**
 * Private method of Accumulator. Extends array of partial sums if required
 *
 * Updates self->partials_cap and self->partials.
 * May move data from self->partials.scalar to self->partials.array
 *
 * \param self an initialized Accumulator
 * \param tail the last index that has to be used
 *
 * \return 0 on success. 1 if memory allocation failed
 */
int _accumulator_reserve(struct Accumulator* self, size_t tail)
{
  if(tail < self->partials_cap)
    return 0;
  size_t reserved = (tail + 1) * 2;
  double* partials_ext;
  if(self->partials_cap == 1) {
    partials_ext = malloc(reserved * sizeof(double));
    if(! partials_ext)
      return 1;
    partials_ext[0] = self->partials.scalar;
  }
  else {
    partials_ext = realloc(self->partials.array, reserved * sizeof(double));
    if(! partials_ext)
      return 1;
  }
  self->partials.array = partials_ext;
  self->partials_cap = reserved;
  return 0;
}

/**
 * Adds a value to the accumulated sum
 *
 * \return 0 on success. 1 if resource allocation failed
 */
int accumulator_add(struct Accumulator* self, double x)
{
  size_t tail;
  x = _accumulator_update_partials(self, x, &tail);
  if(_accumulator_reserve(self, tail))
    return 1;
  self->partials_n = tail + 1;
  if(self->partials_cap > 1)
    self->partials.array[tail] = x;
  else
    self->partials.scalar = x;
  return 0;
}

/**
 * Returns the sum of all added values
 *
 * Note that this call has some overhead. Results should be cached
 */
double accumulator_sum(const struct Accumulator* self)
{
  if(self->partials_cap == 1)
    return self->partials.scalar;
  double sum = 0.;
  size_t i;
  for(i = 0; i < self->partials_n; ++i)
    sum += self->partials.array[i];
  return sum;
}


/**
 * Kahan summation algorithm
 */
struct Kahan
{
  double sum, compensation;
};

void kahan_init(struct Kahan* self)
{
  self->sum = self->compensation = 0.;
}

void kahan_add(struct Kahan* self, double added)
{
  double y = added - self->compensation;
  double t = self->sum + y;
  self->compensation = (t - self->sum) - y;
  self->sum = t;
}

double kahan_sum(const struct Kahan* self)
{
  return self->sum;
}

void test_shewchuk(const double* items, size_t n, size_t repetitions)
{
  struct Accumulator accum;
  accumulator_init(&accum);
  size_t rep_i;
  for(rep_i = 0; rep_i < repetitions; ++rep_i) {
    size_t i;
    for(i = 0; i < n; ++i) {
      if(accumulator_add(&accum, items[i])) {
    printf("OOM\n");
    goto dtor;
      }
    }
  }
  printf("Shewchuk: %g\n", accumulator_sum(&accum));
 dtor:
  accumulator_destroy(&accum);
}

void test_kahan(const double* items, size_t n, size_t repetitions)
{
  struct Kahan kahan;
  kahan_init(&kahan);
  size_t rep_i;
  for(rep_i = 0; rep_i < repetitions; ++rep_i) {
    size_t i;
    for(i = 0; i < n; ++i)
      kahan_add(&kahan, items[i]);
  }
  printf("Kahan: %g\n", kahan_sum(&kahan));
}

void test_devnull(const double* items, size_t n)
{
  FILE* fd = fopen("/dev/null", "w");
  if(! fd)
    return;
  fwrite(items, sizeof(*items), n, fd);
  fclose(fd);
}

static void print_test(const double* items, size_t n_items, size_t repetitions)
{
  printf("Testing summation of {");
  if(n_items) {
    printf("%g", items[0]);
    size_t i;
    for(i = 1; i < n_items; ++i)
      printf(", %g", items[i]);
  }
  printf("} * %zd\n", repetitions);
}

#if ! (defined(WITH_SHEW) || defined(WITH_KAHAN) || defined(WITH_DEVNULL))
#  error No algorithm defined.
#  error Use -DWITH_SHEW and/or -DWITH_KAHAN, or -DWITH_DEVNULL
#endif

int main()
{
  static const double items[] = {1., 1e16, 1., -1e16};
  const size_t n_items = sizeof(items) / sizeof(items[0]);
#ifdef REPETITIONS
  const size_t repetitions = REPETITIONS;
#else
  const size_t repetitions = 100000000;
#endif
  print_test(items, n_items, repetitions);
#ifdef WITH_SHEW
  test_shewchuk(items, n_items, repetitions);
#endif
#ifdef WITH_KAHAN
  test_kahan(items, n_items, repetitions);
#endif
#ifdef WITH_DEVNULL
  test_devnull(items, n_items);
#endif
  return 0;
}