#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>

using namespace std ;


// examples of problematic sum lists
float list1[] = {1, 1e15, 1, -1e15, -1, 1e15, 1, 1e15, 1, -1e15,-1, -1e15 } ;
float list2[] = {1, 1e15, 1, -1e30,  1, 1e15, 1, 1e30, 1, -1e15, 1, -1e15, 1e-15, -1e-15 } ;

//---------------------------------------------------------------------------
// Typical naive summation routine using floating point data
//---------------------------------------------------------------------------
float naive_sum(float l[], int n) {
  float s = l[0] ;
  for(int i=1;i<n;++i)
    s += l[i] ;
  return s ;
}

//---------------------------------------------------------------------------
// Typical naive summation routine using a double precision summation register
//---------------------------------------------------------------------------
float naive_sumd(float l[], int n) {
  double s = l[0] ;
  for(int i=1;i<n;++i)
    s += l[i] ;
  return float(s) ;
}

//---------------------------------------------------------------------------
// Kahan compensated summation example
//---------------------------------------------------------------------------
float kahan_sum(float l[], int n) {
  float sum = 0.0 ; // running sum
  float c = 0.0 ; // Comensation factor (round off error)
  for(int i=0;i<n;++i) {
    float y = l[i]-c ; // add in compensation
    float t = sum + y ; // compute new total
    c = (t - sum) - y ; // estiamte error in sum for compensation
    sum = t ;
  }
  return sum ;
}

//---------------------------------------------------------------------------
// Neumaier version of a compensated sum
//---------------------------------------------------------------------------
float neumaier_sum(float l[], int n) {
  float sum = l[0] ;
  float c= 0 ;  // Running estimate of roundoff error 
  for(int i=1;i<n;++i) {
    float t = sum + l[i] ; // Compute sum
    if(fabs(sum) > fabs(l[i])) // Estimate roundoff error component
      c += (sum - t) + l[i] ;
    else
      c += (l[i] - t) + sum ;
    sum = t ;
  }
  return sum + c ; // Final sum is the combination of roundoff and total
}

//---------------------------------------------------------------------------
// Pairwise summation approach has lower error, break sum into a tree
// this makes the roundoff error growth only log(n)
//---------------------------------------------------------------------------
float pair_sum(float l[], int n) {
  vector<float> scratch(n) ;
  for(int i=0;i<n;++i)
    scratch[i] = l[i] ;
  while(n > 1) { // Repeat adding consecutive pairs until only 1 number remains
    //#define SHOW
#ifdef SHOW
    for(int i=0;i<n;++i)
      cout << scratch[i] << ' ' ;
    cout << endl ;
#endif
    int s = 0 ;
    // loop over pairs
    for(int i=0;i<n;i+=2) {
      if(i+1 < n)  // If pair sum
	scratch[s] = scratch[i]+scratch[i+1] ;
      else
	scratch[s] = scratch[i] ; // else copy
      s++ ;
    }
    n = s; // List now half the size
  }
  return scratch[0] ;
}

bool fabsop(float a, float b) {
  return fabs(a)<fabs(b) ;
}
	  
//---------------------------------------------------------------------------
// Pairwise summation error is even lower if input is sorted according to 
// number magnitudes.
//---------------------------------------------------------------------------
float pair_sum_sort(float l[], int n) {
  vector<float> scratch(n) ;
  for(int i=0;i<n;++i)
    scratch[i] = l[i] ;
  sort(scratch.begin(),scratch.end(),fabsop) ;
  return pair_sum(&scratch[0],n) ;
}

int main() {
  int l1sz = sizeof(list1)/sizeof(float) ;
  int l2sz = sizeof(list2)/sizeof(float) ;
  cout << "naive_sum1 = " << naive_sum(list1,l1sz) << endl ;
  cout << "naive_sum2 = " << naive_sum(list2,l2sz) << endl ;
  cout << "naive_sumd1 = " << naive_sumd(list1,l1sz) << endl ;
  cout << "naive_sumd2 = " << naive_sumd(list2,l2sz) << endl ;
  cout << "kahan_sum1 = " << kahan_sum(list1,l1sz) << endl ;
  cout << "kahan_sum2 = " << kahan_sum(list2,l2sz) << endl ;
  cout << "neumaier_sum1 = " << neumaier_sum(list1,l1sz) << endl ;
  cout << "neumaier_sum2 = " << neumaier_sum(list2,l2sz) << endl ;
  cout << "pair_sum1 = " << pair_sum(list1,l1sz) << endl ;
  cout << "pair_sum2 = " << pair_sum(list2,l2sz) << endl ;
  cout << "pair_sum_sort1 = " << pair_sum_sort(list1,l1sz) << endl ;
  cout << "pair_sum_sort2 = " << pair_sum_sort(list2,l2sz) << endl ;
}