Schewchuk float

/*
 * Copyright 2013-2015 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 <stdio.h>
/* using printf */
#include <stdlib.h>
/* using malloc, free */
#include <math.h>
/* using fabsf */


/**
 * 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 { float* array; float 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.f;
}

/**
 * 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.f;
}

/**
 * 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 float _accumulator_update_partials(struct Accumulator* self,
                      float x,
                      size_t* tail)
{
  *tail = 0;
  float* 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) {
    float y = partials[partial_i];
    if(fabsf(x) < fabsf(y)) {
      float tmp = y;
      y = x;
      x = tmp;
    }
    float hi = x + y;
    float lo = y - (hi - x);
    if(lo != 0.f) {
      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
 */
static int _accumulator_reserve(struct Accumulator* self, size_t tail)
{
  if(tail < self->partials_cap)
    return 0;
  size_t reserved = (tail + 1) * 2;
  float* partials_ext;
  if(self->partials_cap == 1) {
    partials_ext = malloc(reserved * sizeof(float));
    if(! partials_ext)
      return 1;
    partials_ext[0] = self->partials.scalar;
  }
  else {
    partials_ext = realloc(self->partials.array, reserved * sizeof(float));
    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, float 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
 */
float accumulator_sum(const struct Accumulator* self)
{
  if(self->partials_cap == 1)
    return self->partials.scalar;
  float sum = 0.f;
  size_t i;
  for(i = 0; i < self->partials_n; ++i)
    sum += self->partials.array[i];
  return sum;
}

int main(void)
{
  float a = 1.0f;
  float b = 100000000.f;
  float c = -100000000.f;
  double d,e,f;
  struct Accumulator sum;

  d = (a + b) + c;
  e = a + (b + c);
  accumulator_init(&sum);
  accumulator_add(&sum, a);
  accumulator_add(&sum, b);
  accumulator_add(&sum, c);
  f = accumulator_sum(&sum);
  accumulator_destroy(&sum);
  printf("d=%20.20lf e=%20.20lf, f=%20.20lf\n",d, e, f);
  return 0;
}