/*
source: http://www.astro.umd.edu/~dcr/Courses/ASTR615/ps2sol/node1.html
(c) Derek C. Richardson 2009-10-07 
 
** kahan.c
** =======
** Computes sums using the Kahan correction algorithm.
**
** Compile as:
**    gcc -Wall -o kahan kahan.c -lm
** with optional flag:
**    -DDOUBLE to use doubles instead of floats
*/

#include <stdio.h>
#include <stdlib.h> /* for rand() */
#include <unistd.h> /* for getpid() */
#include <time.h> /* for time() */
#include <math.h>
#include <assert.h>

/* construct REAL "type," depending on desired precision */

#ifdef DOUBLE
#define REAL double
#else
#define REAL float
#endif

/* prototype and macro definition for absolute value function */

REAL ABS(REAL);
#ifdef DOUBLE
#define ABS(x) fabs(x)
#else
#define ABS(x) fabsf(x) /* single-precision version */
#endif

REAL sum_kahan(const REAL f[],int N)
{
    /*
    ** Implementation of the Kahan summation algorithm, taken from the
    ** following Wikipedia entry:
    ** http://en.wikipedia.org/wiki/Kahan_summation_algorithm
    */

    REAL sum = f[0];
    REAL c = 0.0,y,t;
    int i;

    for (i=1;i<N;i++) {
        y = f[i] - c;
        t = sum + y;
        c = (t - sum) - y;
        sum = t;
        }

    return sum;
    }

REAL sum_brute(const REAL f[],int N)
{
    /* brute force summation */

    REAL sum = 0.0;
    int i;

    for (i=0;i<N;i++)
        sum += f[i];

    return sum;
    }

void fill_array(REAL f[],int N)
{
    /* fill array with random values between 0 and 1 */

    int i;

    for (i=0;i<N;i++)
        f[i] = (REAL) rand()/RAND_MAX;
    }

void show_sums(const REAL f[],int N)
{
    REAL rSumBrute = sum_brute(f,N);
    REAL rSumKahan = sum_kahan(f,N);

    assert(rSumKahan > 0.0); /* to ensure no division by zero */

#ifdef DOUBLE
    printf("Brute force sum = %.16e, Kahan = %.16e.\n"
           "Fractional difference = %.16e\n",
           rSumBrute,rSumKahan,ABS(rSumKahan - rSumBrute)/rSumKahan);
#else
    printf("Brute force sum = %.7e, Kahan = %.7e.\n"
           "Fractional difference = %.7e.\n",
           rSumBrute,rSumKahan,ABS(rSumKahan - rSumBrute)/rSumKahan);
#endif

    }

int compar_asc(const void *v1,const void *v2)
{
    /* comparison function for qsort(); cf. sort_array() */

    REAL r1 = *((REAL *) v1),r2 = *((REAL *) v2);

    return r1 > r2 ? 1 : r1 < r2 ? -1 : 0; /* ascending order */
    }

int compar_des(const void *v1,const void *v2)
{
    /* comparison function for qsort(); cf. sort_array() */

    REAL r1 = *((REAL *) v1),r2 = *((REAL *) v2);

    return r1 > r2 ? -1 : r1 < r2 ? 1 : 0; /* descending order */
    }

void sort_array(REAL f[],int N,int iDir)
{
    /* sort array in ascending (iDir == 1) or descending (-1) order */

    qsort((void *) f,N,sizeof(REAL),iDir == 1 ? compar_asc : compar_des);
    }

int main(int argc,char *argv[])
{
    printf("Kahan summation algorithm.\n");

    /* check arguments */

    if (argc != 2) {
        fprintf(stderr,"Usage: %s array-size.\n",argv[0]);
        return 1;
        }

    int N = atoi(argv[1]);

    if (N < 1) {
        fprintf(stderr,"Array size must be positive.\n");
        return 1;
        }

    /* show relevant compile-time options, if any */

    printf("Compile-time options: ");
#ifdef DOUBLE
    printf("DOUBLE");
#else
    printf("none");
#endif
    printf(".\n");

    /* seed the random number generator based on the current time and process IDs */

    int iSeed =  time(NULL)%getpid() + getppid();
    srand(iSeed);
    printf("Random number seed = %i.\n",iSeed);

    /* assign space to array */

    REAL f[N];
    printf("Input: array size = %i.\n",N);

    /* generate then sum array in different ways */

    fill_array(f,N);
    printf("Unsorted: ");
    show_sums(f,N);
    sort_array(f,N,1);
    printf("Sorted ascending: ");
    show_sums(f,N);
    sort_array(f,N,-1);
    printf("Sorted descending: ");
    show_sums(f,N);

    return 0;
    }