QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Kahan sums 2 

Former-commit-id: 0d8c588c1770d76fd7db37e003e357f0b1bc4988
This commit is contained in:
Marek Nečada 2018-03-13 10:54:32 +02:00
parent cf340560af
commit ad2815434a
3 changed files with 19 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
qpms/kahansum.h
View File

@ -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;

6
qpms/translations.c
View File

@ -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;

21
qpms/vectors.h
View File

@ -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;