From ad2815434ac1b5391908d6e5464dfda1ab07c66f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Tue, 13 Mar 2018 10:54:32 +0200
Subject: [PATCH] Kahan sums 2

Former-commit-id: 0d8c588c1770d76fd7db37e003e357f0b1bc4988
---
 qpms/kahansum.h     |  4 ++--
 qpms/translations.c |  6 +++---
 qpms/vectors.h      | 21 ++++++++++++++-------
 3 files changed, 19 insertions(+), 12 deletions(-)

diff --git a/qpms/kahansum.h b/qpms/kahansum.h
index 466c3dd..66ef1c1 100644
--- a/qpms/kahansum.h
+++ b/qpms/kahansum.h
@@ -6,7 +6,7 @@ static inline void kahaninit(double *sum, double *compensation) {
 	compensation = 0;
 }
 
-static inline void kahaninc(double *sum, double *compensation, double input) {
+static inline void kahanadd(double *sum, double *compensation, double input) {
 	double compensated_input = input - *compensation;
 	double nsum = *sum + compensated_input;
 	*compensation = (nsum - *sum) - compensated_input;
@@ -19,7 +19,7 @@ static inline void ckahaninit(complex double *sum, complex double *compensation)
 	compensation = 0;
 }
 
-static inline void ckahaninc(complex double *sum, complex double *compensation, complex double input) {
+static inline void ckahanadd(complex double *sum, complex double *compensation, complex double input) {
 	complex double compensated_input = input - *compensation;
 	complex double nsum = *sum + compensated_input;
 	*compensation = (nsum - *sum) - compensated_input;
diff --git a/qpms/translations.c b/qpms/translations.c
index c6f2b03..f493368 100644
--- a/qpms/translations.c
+++ b/qpms/translations.c
@@ -820,7 +820,7 @@ static inline complex double qpms_trans_calculator_get_A_precalcbuf(const qpms_t
 		double Pp = legendre_buf[gsl_sf_legendre_array_index(p, abs(mu-m))];
 		complex double zp = bessel_buf[p];
 		complex double multiplier = c->A_multipliers[i][q];
-		ckahaninc(&sum, &kahanc, Pp * zp *  multiplier);
+		ckahanadd(&sum, &kahanc, Pp * zp *  multiplier);
 	}
 	complex double eimf =  cexp(I*(mu-m)*kdlj.phi);
 	return sum * eimf;
@@ -864,13 +864,13 @@ static inline complex double qpms_trans_calculator_get_B_precalcbuf(const qpms_t
 	size_t qmax = c->B_multipliers[i+1] - c->B_multipliers[i] - (1 - BQ_OFFSET);
 	assert(qmax == gauntB_Q_max(-m,n,mu,nu));
 	complex double sum, kahanc;
-	kahaninit(&sum, &kahanc);
+	ckahaninit(&sum, &kahanc);
 	for(int q = BQ_OFFSET; q <= qmax; ++q) {
 		int p = n+nu-2*q;
 		double Pp_ = legendre_buf[gsl_sf_legendre_array_index(p+1, abs(mu-m))];
 		complex double zp_ = bessel_buf[p+1];
 		complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
-		ckahaninc(&sum, &kahanc, Pp_ * zp_ * multiplier);
+		ckahanadd(&sum, &kahanc, Pp_ * zp_ * multiplier);
 	}
 	complex double eimf =  cexp(I*(mu-m)*kdlj.phi);
 	return sum * eimf;
diff --git a/qpms/vectors.h b/qpms/vectors.h
index 88bea09..1fd27cf 100644
--- a/qpms/vectors.h
+++ b/qpms/vectors.h
@@ -49,6 +49,16 @@ static inline ccart3_t ccart3_scale(const complex  double c, const ccart3_t v) {
 	return res;
 }
 
+static inline ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b) {
+	ccart3_t res = {a.x+b.x, a.y+b.y, a.z+b.z};
+	return res;
+}
+
+static inline ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) {
+	ccart3_t res = {a.x-b.x, a.y-b.y, a.z-b.z};
+	return res;
+}
+
 static inline csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) {
 	csphvec_t res = {a.rc + b.rc, a.thetac + b.thetac, a.phic + b.phic};
 	return res;
@@ -131,18 +141,15 @@ static inline void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation) {
 	sum->rc = sum->thetac = sum->phic = compensation->rc = compensation->thetac = compensation->phic = 0;
 }
 
-static inline void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const *ccart3_t input) {
-	ccart3_t comped_input = ccart3_substract(*input, *compensation);
+static inline void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input) {
+	ccart3_t comped_input = ccart3_substract(input, *compensation);
 	ccart3_t nsum = ccart3_add(*sum, comped_input);
 	*compensation = ccart3_substract(ccart3_substract(nsum, *sum), comped_input);
 	*sum = nsum;
 }
 
-
-
-
-static inline void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t *input) {
-	csphvec_t comped_input = csphvec_substract(*input, *compensation);
+static inline void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input) {
+	csphvec_t comped_input = csphvec_substract(input, *compensation);
 	csphvec_t nsum = csphvec_add(*sum, comped_input);
 	*compensation = csphvec_substract(csphvec_substract(nsum, *sum), comped_input);
 	*sum = nsum;