Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /*
- * 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;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement