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