/*
Copyright 2016 Ian Jauslin

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 "ss_integral.h"
#include <pthread.h>

// compute the integral
// takes a pointer to an integrand function 'I' as an argument. 'I' must have the following arguments
//   int I(mpfr_t out, mpfr_t p1, mpfr_t p2, mpfr_t q1, mpfr_t q2, mpfr_t F1, mpfr_t F2, mpfr_t t1, mpfr_t t2, mpfr_t phi, mpfr_t W, array_mpfr* tmps, struct ss_cache cache)
int ss_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, unsigned int threads){
  mpfr_t val;
  array_2_mpfr tmpss;
  struct ss_cache cache;
  unsigned int i;

  // compute pi and sqrt3
  mpfr_inits(cache.pi, cache.sqrt3, NULL);
  mpfr_const_pi(cache.pi, MPFR_RNDN);
  mpfr_sqrt_ui(cache.sqrt3, 3, MPFR_RNDN);

  mpfr_init(val);

  // init tmps
  array_2_mpfr_init(&tmpss, 8*threads);
  tmpss.length=8*threads;
  for(i=0;i<threads;i++){
    // 8 array_mpfr's per thread
    // do not allocate too much memory (I'm trying to keep a modicum of optimality)
    array_mpfr_init(tmpss.values+8*i+0, 2);
    array_mpfr_init(tmpss.values+8*i+1, 4);
    array_mpfr_init(tmpss.values+8*i+2, 3);
    array_mpfr_init(tmpss.values+8*i+3, 4);
    array_mpfr_init(tmpss.values+8*i+4, 2);
    array_mpfr_init(tmpss.values+8*i+5, 4);
    array_mpfr_init(tmpss.values+8*i+6, 12);
    array_mpfr_init(tmpss.values+8*i+7, 9);
  }

  // compute A_FF, A_RF and A_RR
  ssFF_integrate(out, I, params, abcissa, weights, cache, tmpss, threads);

  ssRF_integrate(&val, I, params, abcissa, weights, cache, tmpss, threads);
  mpfr_mul_ui(val, val, 2, MPFR_RNDN);
  mpfr_add(*out, *out, val, MPFR_RNDN);

  ssRR_integrate(&val, I, params, abcissa, weights, cache, tmpss, threads);
  mpfr_add(*out, *out, val, MPFR_RNDN);

  mpfr_clears(cache.pi, cache.sqrt3, NULL);
  mpfr_clear(val);

  array_2_mpfr_free(tmpss);
  return(0);
}


// compute the integral A_FF
int ssFF_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads){
  ssFF_argsint_rho args;
  int ret;
  mpfr_t upper, lower;
  array_pthread_t thread_ids;

  // init
  thread_ids.length=0;

  args.params=params;
  args.I=I;
  args.abcissa=abcissa;
  args.weights=weights;
  args.cache=cache;
  args.tmpss=tmpss;
  args.thread_ids=&thread_ids;

  // bounds
  mpfr_inits(upper, lower, NULL);
  mpfr_set_ui(lower, 0, MPFR_RNDN);
  mpfr_set_ui(upper, 1, MPFR_RNDN);

  if(threads==1){
    ret=integrate_gauss_mpfr(out, &ssFF_integrand_rho, lower, upper, abcissa, weights, &args);
  }
  else{
    ret=integrate_gauss_multithread_mpfr(out, &ssFF_integrand_rho, lower, upper, abcissa, weights, &args, threads, args.thread_ids);
    array_pthread_t_free(*(args.thread_ids));
  }

  // extra factor
  mpfr_mul_ui(*out, *out, 6, MPFR_RNDN);

  mpfr_clears(upper, lower, NULL);

  return(ret);
}

// integrand of the integral over rho in A_FF
int ssFF_integrand_rho(mpfr_t* out, mpfr_t rho, void* args){
  ssFF_argsint_theta nargs;
  int ret;
  ssFF_argsint_rho* argument=(ssFF_argsint_rho*)args;
  array_mpfr* tmps1;
  array_mpfr* tmps2;
  int thread_id;

  // the number of the thread
  if(argument->thread_ids->length>1){
    thread_id=array_pthread_t_find(pthread_self(), *(argument->thread_ids));
    if(thread_id<0){
      fprintf(stderr, "error: thread not found\nIf you see this message, then you've uncovered a serious bug...\n");
      return(-1);
    }
  }
  else{
    thread_id=0;
  }

  // tmps used in this function (not passed to integrand_theta)
  tmps1=argument->tmpss.values+8*thread_id;
  tmps2=argument->tmpss.values+8*thread_id+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.abcissa=argument->abcissa;
  nargs.weights=argument->weights;
  nargs.cache=argument->cache;
  nargs.rho=rho;

  // copy tmps
  // only pass tmps reserved for this thread
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+8*thread_id+2;
  // 6 tmps in integrand_theta
  nargs.tmpss.length=6;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_div_si(tmps1->values[0], nargs.cache.pi, -6, MPFR_RNDN);
  mpfr_div_ui(tmps1->values[1], nargs.cache.pi, 2, MPFR_RNDN);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssFF_integrand_theta, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over theta in A_FF
int ssFF_integrand_theta(mpfr_t* out, mpfr_t theta, void* args){
  ssFF_argsint_psi nargs;
  int ret;
  ssFF_argsint_theta* argument=(ssFF_argsint_theta*)args;
  unsigned int a;
  // tmps used in this function (not passed to integrand_psi)
  array_mpfr* tmps1=argument->tmpss.values;
  array_mpfr* tmps2=argument->tmpss.values+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.abcissa=argument->abcissa;
  nargs.weights=argument->weights;
  nargs.cache=argument->cache;
  nargs.rho=argument->rho;
  nargs.theta=theta;

  // copy tmps
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+2;
  // 4 tmps in integrand_psi
  nargs.tmpss.length=4;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(3, tmps1);

  // split the integral into four parts to improve the precision
  mpfr_set_ui(*out, 0, MPFR_RNDN);
  for(a=0;a<4;a++){
    // bounds
    mpfr_set_ui(tmps1->values[0], a, MPFR_RNDN);
    mpfr_sub_ui(tmps1->values[0], tmps1->values[0], 2, MPFR_RNDN);
    mpfr_div_ui(tmps1->values[0], tmps1->values[0], 4, MPFR_RNDN);
    mpfr_mul(tmps1->values[0], tmps1->values[0], nargs.cache.pi, MPFR_RNDN);

    mpfr_set_ui(tmps1->values[1], a, MPFR_RNDN);
    mpfr_sub_ui(tmps1->values[1], tmps1->values[1], 1, MPFR_RNDN);
    mpfr_div_ui(tmps1->values[1], tmps1->values[1], 4, MPFR_RNDN);
    mpfr_mul(tmps1->values[1], tmps1->values[1], nargs.cache.pi, MPFR_RNDN);


    ret=integrate_gauss_smarttmp_mpfr(tmps1->values+2, &ssFF_integrand_psi, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);
    if(ret<0){
      break;
    }
    mpfr_add(*out, *out, tmps1->values[2], MPFR_RNDN);
  }

  return(ret);
}

// integrand of the integral over psi in A_FF
int ssFF_integrand_psi(mpfr_t* out, mpfr_t psi, void* args){
  ssFF_argsint_z nargs;
  int ret;
  ssFF_argsint_psi* argument=(ssFF_argsint_psi*)args;
  // tmps used in this function (not passed to integrand_z)
  array_mpfr* tmps1=argument->tmpss.values;
  array_mpfr* tmps2=argument->tmpss.values+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.cache=argument->cache;
  nargs.rho=argument->rho;
  nargs.theta=argument->theta;
  nargs.psi=psi;

  // copy tmps
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+2;
  // 2 tmps in integrand_z
  nargs.tmpss.length=2;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN);
  mpfr_set_ui(tmps1->values[1], 1, MPFR_RNDN);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssFF_integrand_z, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over z in A_FF
int ssFF_integrand_z(mpfr_t* out, mpfr_t z, void* args){
  ssFF_argsint_z* a=(ssFF_argsint_z*)args;
  mpfr_t* tmps;
  // tmps used in this function
  array_mpfr* tmps1=a->tmpss.values;
  array_mpfr* tmps2=a->tmpss.values+1;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(7, tmps1);

  tmps=tmps1->values;

  // r (store in tmps[0])
  mpfr_cos(tmps[0], a->psi, MPFR_RNDN);
  mpfr_add_ui(tmps[1], tmps[0], 1, MPFR_RNDN);
  mpfr_ui_sub(tmps[0], 1, tmps[0], MPFR_RNDN);
  mpfr_div(tmps[0], tmps[0], tmps[1], MPFR_RNDN);
  mpfr_mul(tmps[0], tmps[0], z, MPFR_RNDN);
  mpfr_add_ui(tmps[0], tmps[0], 1, MPFR_RNDN);
  mpfr_ui_sub(tmps[1], 1, z, MPFR_RNDN);
  mpfr_mul(tmps[0], tmps[0], tmps[1], MPFR_RNDN);

  // varphi (store in tmps[1])
  mpfr_mul_ui(tmps[1], a->psi, 2, MPFR_RNDN);
  mpfr_sin(tmps[2], tmps[1], MPFR_RNDN);
  mpfr_cos(tmps[1], tmps[1], MPFR_RNDN);
  mpfr_pow_ui(tmps[2], tmps[2], 2, MPFR_RNDN);
  mpfr_mul(tmps[2], tmps[2], tmps[0], MPFR_RNDN);
  mpfr_div_ui(tmps[2], tmps[2], 2, MPFR_RNDN);
  mpfr_sub(tmps[1], tmps[1], tmps[2], MPFR_RNDN);
  mpfr_acos(tmps[1], tmps[1], MPFR_RNDN);
  if(mpfr_cmp_ui(a->psi, 0)<0){
    mpfr_neg(tmps[1], tmps[1], MPFR_RNDN);
  }

  // jacobian
  mpfr_mul_ui(*out, a->psi, 2, MPFR_RNDN);
  mpfr_cos(*out, *out, MPFR_RNDN);
  mpfr_mul(*out, *out, tmps[0], MPFR_RNDN);
  mpfr_add_ui(*out, *out, 1, MPFR_RNDN);
  mpfr_cos(tmps[2], a->psi, MPFR_RNDN);
  mpfr_pow_ui(tmps[3], tmps[2], 2, MPFR_RNDN);
  mpfr_mul(tmps[3], tmps[3], tmps[0], MPFR_RNDN);
  mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN);
  mpfr_sqrt(tmps[3], tmps[3], MPFR_RNDN);
  mpfr_div(*out, *out, tmps[3], MPFR_RNDN);
  mpfr_add_ui(tmps[2], tmps[2], 1, MPFR_RNDN);
  mpfr_div(*out, *out, tmps[2], MPFR_RNDN);
  mpfr_mul(*out, *out, tmps[0], MPFR_RNDN);
  mpfr_mul(*out, *out, a->rho, MPFR_RNDN);
  mpfr_mul(*out, *out, a->rho, MPFR_RNDN);
  mpfr_mul(*out, *out, a->rho, MPFR_RNDN);
  mpfr_mul_ui(*out, *out, 4, MPFR_RNDN);

  // cutoffs
  ss_cutoff(tmps[2], a->rho, tmps[3], tmps[4]);
  mpfr_mul(*out, *out, tmps[2], MPFR_RNDN);
  mpfr_mul(tmps[2], a->rho, tmps[0], MPFR_RNDN);
  ss_cutoff(tmps[2], tmps[2], tmps[3], tmps[4]);
  mpfr_mul(*out, *out, tmps[2], MPFR_RNDN);

  // q1 (store in tmps[2])
  mpfr_add(tmps[2], a->theta, tmps[1], MPFR_RNDN);
  mpfr_cos(tmps[2], tmps[2], MPFR_RNDN);
  mpfr_mul(tmps[2], tmps[2], tmps[0], MPFR_RNDN);
  mpfr_mul(tmps[2], tmps[2], a->rho, MPFR_RNDN);
  mpfr_add(tmps[2], tmps[2], a->cache.sqrt3, MPFR_RNDN);
  // q2 (store in tmps[3])
  mpfr_add(tmps[3], a->theta, tmps[1], MPFR_RNDN);
  mpfr_sin(tmps[3], tmps[3], MPFR_RNDN);
  mpfr_mul(tmps[3], tmps[3], tmps[0], MPFR_RNDN);
  mpfr_mul(tmps[3], tmps[3], a->rho, MPFR_RNDN);
  mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN);
  if(a->params.omega==-1){
    mpfr_neg(tmps[3], tmps[3], MPFR_RNDN);
  }
  // p1 (store in tmps[0])
  mpfr_cos(tmps[0], a->theta, MPFR_RNDN);
  mpfr_mul(tmps[0], tmps[0], a->rho, MPFR_RNDN);
  mpfr_add(tmps[0], tmps[0], a->cache.sqrt3, MPFR_RNDN);
  // p2 (store in tmps[1])
  mpfr_sin(tmps[1], a->theta, MPFR_RNDN);
  mpfr_mul(tmps[1], tmps[1], a->rho, MPFR_RNDN);
  mpfr_add_ui(tmps[1], tmps[1], 1, MPFR_RNDN);
  if(a->params.omega==-1){
    mpfr_neg(tmps[1], tmps[1], MPFR_RNDN);
  }
  // F1 (store in tmps[4])
  mpfr_add(tmps[4], tmps[0], tmps[2], MPFR_RNDN);
  mpfr_sub(tmps[4], tmps[4], a->cache.sqrt3, MPFR_RNDN);
  // F2 (store in tmps[5])
  mpfr_add(tmps[5], tmps[1], tmps[3], MPFR_RNDN);
  mpfr_sub_si(tmps[5], tmps[5], a->params.omega, MPFR_RNDN);

  // I
  (*(a->I))(tmps[6], tmps[0], tmps[1], tmps[2], tmps[3], tmps[4], tmps[5], a->params.t1, a->params.t2, a->params.phi, a->params.W, tmps2, a->cache);
  mpfr_mul(*out, *out, tmps[6], MPFR_RNDN);

  return(0);
}


// compute the integral A_RF
int ssRF_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads){
  ssRF_argsint_theta args;
  int ret;
  mpfr_t val;
  mpfr_t lower, upper;
  unsigned int a;
  array_pthread_t thread_ids;

  // init
  thread_ids.length=0;

  args.params=params;
  args.I=I;
  args.abcissa=abcissa;
  args.weights=weights;
  args.cache=cache;
  args.tmpss=tmpss;
  args.thread_ids=&thread_ids;

  mpfr_inits(val, lower, upper, NULL);

  // split the integral in three parts to improve the precision
  mpfr_set_ui(*out, 0, MPFR_RNDN);
  for(a=0;a<3;a++){
    // bounds
    mpfr_set_ui(lower, a, MPFR_RNDN);
    mpfr_mul_ui(lower, lower, 4, MPFR_RNDN);
    mpfr_sub_ui(lower, lower, 1, MPFR_RNDN);
    mpfr_div_ui(lower, lower, 6, MPFR_RNDN);
    mpfr_mul(lower, lower, args.cache.pi, MPFR_RNDN);

    mpfr_set_ui(upper, a, MPFR_RNDN);
    mpfr_mul_ui(upper, upper, 4, MPFR_RNDN);
    mpfr_add_ui(upper, upper, 3, MPFR_RNDN);
    mpfr_div_ui(upper, upper, 6, MPFR_RNDN);
    mpfr_mul(upper, upper, args.cache.pi, MPFR_RNDN);

    if(threads==1){
      ret=integrate_gauss_mpfr(&val, &ssRF_integrand_theta, lower, upper, abcissa, weights, &args);
    }
    else{
      ret=integrate_gauss_multithread_mpfr(&val, &ssRF_integrand_theta, lower, upper, abcissa, weights, &args, threads, args.thread_ids);
      array_pthread_t_free(*(args.thread_ids));
    }
    if(ret<0){
      break;
    }
    
    mpfr_add(*out, *out, val, MPFR_RNDN);
  }

  // extra factor
  mpfr_mul_ui(*out, *out, 3, MPFR_RNDN);

  mpfr_clears(val, lower, upper, NULL);
  return(ret);
}

// integrand of the integral over theta in A_RF
int ssRF_integrand_theta(mpfr_t* out, mpfr_t theta, void* args){
  ssRF_argsint_varphi nargs;
  int ret;
  ssRF_argsint_theta* argument=(ssRF_argsint_theta*)args;
  // tmps used in this function (not passed to integrand_varphi)
  array_mpfr* tmps1;
  array_mpfr* tmps2;
  int thread_id;

  // the number of the thread
  if(argument->thread_ids->length>1){
    thread_id=array_pthread_t_find(pthread_self(), *(argument->thread_ids));
    if(thread_id<0){
      fprintf(stderr, "error: thread not found\nIf you see this message, then you've uncovered a serious bug...\n");
      return(-1);
    }
  }
  else{
    thread_id=0;
  }

  // tmps used in this function (not passed to integrand_varphi)
  tmps1=argument->tmpss.values+8*thread_id;
  tmps2=argument->tmpss.values+8*thread_id+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.abcissa=argument->abcissa;
  nargs.weights=argument->weights;
  nargs.cache=argument->cache;
  nargs.theta=theta;

  // copy tmps
  // only pass tmps reserved for this thread
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+8*thread_id+2;
  // 6 tmps in integrand_varphi
  nargs.tmpss.length=6;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_div_si(tmps1->values[0], nargs.cache.pi, -6, MPFR_RNDN);
  mpfr_div_ui(tmps1->values[1], nargs.cache.pi, 2, MPFR_RNDN);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssRF_integrand_varphi, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over varphi in A_RF
int ssRF_integrand_varphi(mpfr_t* out, mpfr_t varphi, void* args){
  ssRF_argsint_r nargs;
  int ret;
  ssRF_argsint_varphi* argument=(ssRF_argsint_varphi*)args;
  // tmps used in this function (not passed to integrand_r)
  array_mpfr* tmps1=argument->tmpss.values;
  array_mpfr* tmps2=argument->tmpss.values+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.abcissa=argument->abcissa;
  nargs.weights=argument->weights;
  nargs.cache=argument->cache;
  nargs.theta=argument->theta;
  nargs.varphi=varphi;

  // copy tmps
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+2;
  // 4 tmps in integrand_r
  nargs.tmpss.length=4;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN);
  mpfr_set_ui(tmps1->values[1], 1, MPFR_RNDN);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssRF_integrand_r, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over r in A_RF
int ssRF_integrand_r(mpfr_t* out, mpfr_t r, void* args){
  ssRF_argsint_rho nargs;
  int ret;
  ssRF_argsint_r* argument=(ssRF_argsint_r*)args;
  // tmps used in this function (not passed to integrand_rho)
  array_mpfr* tmps1=argument->tmpss.values;
  array_mpfr* tmps2=argument->tmpss.values+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.cache=argument->cache;
  nargs.theta=argument->theta;
  nargs.varphi=argument->varphi;
  nargs.r=r;

  // copy tmps
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+2;
  // 2 tmps in integrand_rho
  nargs.tmpss.length=2;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN);
  ss_R(tmps1->values[1], nargs.varphi, 0, nargs.cache);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssRF_integrand_rho, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over rho in A_RF
int ssRF_integrand_rho(mpfr_t* out, mpfr_t rho, void* args){
  int nu;
  ssRF_argsint_rho* a=(ssRF_argsint_rho*)args;
  mpfr_t* tmps;
  // tmps used in this function
  array_mpfr* tmps1=a->tmpss.values;
  array_mpfr* tmps2=a->tmpss.values+1;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(11, tmps1);

  tmps=tmps1->values;

  mpfr_set_ui(*out, 0, MPFR_RNDN);
  for(nu=-1;nu<=1;nu=nu+2){
    // q1 (store in tmps[2])
    mpfr_cos(tmps[2], a->theta, MPFR_RNDN);
    mpfr_mul(tmps[2], tmps[2], a->r, MPFR_RNDN);
    mpfr_add(tmps[2], tmps[2], a->cache.sqrt3, MPFR_RNDN);
    // q2 (store in tmps[3])
    mpfr_sin(tmps[3], a->theta, MPFR_RNDN);
    mpfr_mul(tmps[3], tmps[3], a->r, MPFR_RNDN);
    mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN);
    if(a->params.omega==-1){
      mpfr_neg(tmps[3], tmps[3], MPFR_RNDN);
    }
    // F1 (store in tmps[4])
    mpfr_cos(tmps[4], a->varphi, MPFR_RNDN);
    mpfr_mul(tmps[4], tmps[4], rho, MPFR_RNDN);
    mpfr_add(tmps[4], tmps[4], a->cache.sqrt3, MPFR_RNDN);
    // F2 (store in tmps[5])
    mpfr_sin(tmps[5], a->varphi, MPFR_RNDN);
    mpfr_mul(tmps[5], tmps[5], rho, MPFR_RNDN);
    mpfr_add_ui(tmps[5], tmps[5], 1, MPFR_RNDN);
    if(nu==-1){
      mpfr_neg(tmps[5], tmps[5], MPFR_RNDN);
    }
    // p1 (store in tmps[0])
    mpfr_sub(tmps[0], tmps[4], tmps[2], MPFR_RNDN);
    mpfr_add(tmps[0], tmps[0], a->cache.sqrt3, MPFR_RNDN);
    // p2 (store in tmps[1])
    mpfr_sub(tmps[1], tmps[5], tmps[3], MPFR_RNDN);
    mpfr_add_si(tmps[1], tmps[1], a->params.omega, MPFR_RNDN);

    // cutoff
    // tmps[6]=(1-ss_cutoff(ss_norm(F1-q1,F2-q2)))*ss_cutoff(r)
    mpfr_sub(tmps[6], tmps[4], tmps[2], MPFR_RNDN);
    mpfr_sub(tmps[7], tmps[5], tmps[3], MPFR_RNDN);
    ss_norm(tmps[8], tmps[6], tmps[7], a->cache, tmps[9], tmps[10]);
    ss_cutoff(tmps[6], tmps[8], tmps[7], tmps[9]);
    mpfr_ui_sub(tmps[6], 1, tmps[6], MPFR_RNDN);
    ss_cutoff(tmps[7], a->r, tmps[8], tmps[9]);
    mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN);

    (*(a->I))(tmps[7], tmps[0], tmps[1], tmps[2], tmps[3], tmps[4], tmps[5], a->params.t1, a->params.t2, a->params.phi, a->params.W, tmps2, a->cache);

    // r*rho*I*cutoffs
    mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN);
    mpfr_mul(tmps[6], tmps[6], a->r, MPFR_RNDN);
    mpfr_mul(tmps[6], tmps[6], rho, MPFR_RNDN);

    // add to *out
    mpfr_add(*out, *out, tmps[6], MPFR_RNDN);
  }

  return(0);
}


// compute the integral A_RR
int ssRR_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads){
  ssRR_argsint_theta args;
  int ret;
  unsigned int a;
  mpfr_t lower, upper;
  mpfr_t val;
  array_pthread_t thread_ids;

  // init
  thread_ids.length=0;

  args.params=params;
  args.I=I;
  args.abcissa=abcissa;
  args.weights=weights;
  args.cache=cache;
  args.tmpss=tmpss;
  args.thread_ids=&thread_ids;

  mpfr_inits(val, lower, upper, NULL);

  // split the integral in three parts
  mpfr_set_ui(*out, 0, MPFR_RNDN);
  for(a=0;a<3;a++){
    // bounds
    mpfr_set_ui(lower, a, MPFR_RNDN);
    mpfr_mul_ui(lower, lower, 4, MPFR_RNDN);
    mpfr_sub_ui(lower, lower, 1, MPFR_RNDN);
    mpfr_div_ui(lower, lower, 6, MPFR_RNDN);
    mpfr_mul(lower, lower, args.cache.pi, MPFR_RNDN);

    mpfr_set_ui(upper, a, MPFR_RNDN);
    mpfr_mul_ui(upper, upper, 4, MPFR_RNDN);
    mpfr_add_ui(upper, upper, 3, MPFR_RNDN);
    mpfr_div_ui(upper, upper, 6, MPFR_RNDN);
    mpfr_mul(upper, upper, args.cache.pi, MPFR_RNDN);

    // save sector
    args.sector_theta=a;

    if(threads==1){
      ret=integrate_gauss_mpfr(&val, &ssRR_integrand_theta, lower, upper, abcissa, weights, &args);
    }
    else{
      ret=integrate_gauss_multithread_mpfr(&val, &ssRR_integrand_theta, lower, upper, abcissa, weights, &args, threads, args.thread_ids);
      array_pthread_t_free(*(args.thread_ids));
    }
    if(ret<0){
      break;
    }
    
    mpfr_add(*out, *out, val, MPFR_RNDN);
  }
  // extra factor
  mpfr_mul_ui(*out, *out, 3, MPFR_RNDN);

  mpfr_clears(val, lower, upper, NULL);

  return(ret);
}

// integrand of the integral over theta in A_RR
int ssRR_integrand_theta(mpfr_t* out, mpfr_t theta, void* args){
  ssRR_argsint_varphi nargs;
  int ret;
  ssRR_argsint_theta* argument=(ssRR_argsint_theta*)args;
  array_mpfr* tmps1;
  array_mpfr* tmps2;
  int thread_id;

  // the number of the thread
  if(argument->thread_ids->length>1){
    thread_id=array_pthread_t_find(pthread_self(), *(argument->thread_ids));
    if(thread_id<0){
      fprintf(stderr, "error: thread not found\nIf you see this message, then you've uncovered a serious bug...\n");
      return(-1);
    }
  }
  else{
    thread_id=0;
  }

  // tmps used in this function (not passed to integrand_varphi)
  tmps1=argument->tmpss.values+8*thread_id;
  tmps2=argument->tmpss.values+8*thread_id+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.abcissa=argument->abcissa;
  nargs.weights=argument->weights;
  nargs.cache=argument->cache;
  nargs.sector_theta=argument->sector_theta;
  nargs.theta=theta;

  // copy tmps
  // only pass tmps reserved for this thread
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+8*thread_id+2;
  // 6 tmps in integrand_varphi
  nargs.tmpss.length=6;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_div_si(tmps1->values[0], nargs.cache.pi, -6, MPFR_RNDN);
  mpfr_div_ui(tmps1->values[1], nargs.cache.pi, 2, MPFR_RNDN);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssRR_integrand_varphi, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over varphi in A_RR
int ssRR_integrand_varphi(mpfr_t* out, mpfr_t varphi, void* args){
  ssRR_argsint_r nargs;
  int ret;
  ssRR_argsint_varphi* argument=(ssRR_argsint_varphi*)args;
  // tmps used in this function (not passed to integrand_r)
  array_mpfr* tmps1=argument->tmpss.values;
  array_mpfr* tmps2=argument->tmpss.values+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.abcissa=argument->abcissa;
  nargs.weights=argument->weights;
  nargs.cache=argument->cache;
  nargs.theta=argument->theta;
  nargs.varphi=varphi;

  // copy tmps
  // only pass tmps reserved for this thread
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+2;
  // 4 tmps in integrand_r
  nargs.tmpss.length=4;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN);
  ss_R(tmps1->values[1], nargs.theta, argument->sector_theta, nargs.cache);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssRR_integrand_r, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over r in A_RR
int ssRR_integrand_r(mpfr_t* out, mpfr_t r, void* args){
  ssRR_argsint_rho nargs;
  int ret;
  ssRR_argsint_r* argument=(ssRR_argsint_r*)args;
  // tmps used in this function (not passed to integrand_rho)
  array_mpfr* tmps1=argument->tmpss.values;
  array_mpfr* tmps2=argument->tmpss.values+1;

  // new argument
  nargs.params=argument->params;
  nargs.I=argument->I;
  nargs.cache=argument->cache;
  nargs.theta=argument->theta;
  nargs.varphi=argument->varphi;
  nargs.r=r;

  // copy tmps
  // only pass tmps reserved for this thread
  // warning: 'nargs.tmpss' must not be freed or resized
  nargs.tmpss.values=argument->tmpss.values+2;
  // 2 tmps in integrand_rho
  nargs.tmpss.length=2;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(2, tmps1);

  // bounds
  mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN);
  ss_R(tmps1->values[1], nargs.varphi, 0, nargs.cache);

  ret=integrate_gauss_smarttmp_mpfr(out, &ssRR_integrand_rho, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs);

  return(ret);
}

// integrand of the integral over rho in A_RR
int ssRR_integrand_rho(mpfr_t* out, mpfr_t rho, void* args){
  int eta, nu;
  ssRR_argsint_rho* a=(ssRR_argsint_rho*)args;
  mpfr_t* tmps;
  // tmps used in this function (not passed to integrand_rho)
  array_mpfr* tmps1=a->tmpss.values;
  array_mpfr* tmps2=a->tmpss.values+1;

  // alloc tmps if needed
  array_mpfr_alloc_tmps(12, tmps1);

  tmps=tmps1->values;

  mpfr_set_ui(*out, 0, MPFR_RNDN);
  for(eta=-1;eta<=1;eta=eta+2){
    for(nu=-1;nu<=1;nu=nu+2){
      // q1 (store in tmps[2])
      mpfr_cos(tmps[2], a->theta, MPFR_RNDN);
      mpfr_mul(tmps[2], tmps[2], a->r, MPFR_RNDN);
      mpfr_add(tmps[2], tmps[2], a->cache.sqrt3, MPFR_RNDN);
      // q2 (store in tmps[3])
      mpfr_sin(tmps[3], a->theta, MPFR_RNDN);
      mpfr_mul(tmps[3], tmps[3], a->r, MPFR_RNDN);
      mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN);
      if(eta==-1){
	mpfr_neg(tmps[3], tmps[3], MPFR_RNDN);
      }
      // F1 (store in tmps[4])
      mpfr_cos(tmps[4], a->varphi, MPFR_RNDN);
      mpfr_mul(tmps[4], tmps[4], rho, MPFR_RNDN);
      mpfr_add(tmps[4], tmps[4], a->cache.sqrt3, MPFR_RNDN);
      // F2 (store in tmps[5])
      mpfr_sin(tmps[5], a->varphi, MPFR_RNDN);
      mpfr_mul(tmps[5], tmps[5], rho, MPFR_RNDN);
      mpfr_add_ui(tmps[5], tmps[5], 1, MPFR_RNDN);
      if(nu==-1){
	mpfr_neg(tmps[5], tmps[5], MPFR_RNDN);
      }
      // p1 (store in tmps[0])
      mpfr_sub(tmps[0], tmps[4], tmps[2], MPFR_RNDN);
      mpfr_add(tmps[0], tmps[0], a->cache.sqrt3, MPFR_RNDN);
      // p2 (store in tmps[1])
      mpfr_sub(tmps[1], tmps[5], tmps[3], MPFR_RNDN);
      mpfr_add_si(tmps[1], tmps[1], a->params.omega, MPFR_RNDN);

      // cutoff
      // tmps[6]=(1-ss_cutoff(ss_norm(F1-q1,F2-q2)))*(1-ss_cutoff(ss_norm(q1-SQRT3,q2-omega)))
      mpfr_sub(tmps[6], tmps[4], tmps[2], MPFR_RNDN);
      mpfr_sub(tmps[7], tmps[5], tmps[3], MPFR_RNDN);
      ss_norm(tmps[8], tmps[6], tmps[7], a->cache, tmps[9], tmps[10]);
      ss_cutoff(tmps[6], tmps[8], tmps[7], tmps[9]);
      mpfr_ui_sub(tmps[6], 1, tmps[6], MPFR_RNDN);
      mpfr_sub(tmps[7], tmps[2], a->cache.sqrt3, MPFR_RNDN);
      mpfr_sub_si(tmps[8], tmps[3], a->params.omega, MPFR_RNDN);
      ss_norm(tmps[9], tmps[7], tmps[8], a->cache, tmps[10], tmps[11]);
      ss_cutoff(tmps[7], tmps[9], tmps[8], tmps[10]);
      mpfr_ui_sub(tmps[7], 1, tmps[7], MPFR_RNDN);
      mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN);

      (*(a->I))(tmps[7], tmps[0], tmps[1], tmps[2], tmps[3], tmps[4], tmps[5], a->params.t1, a->params.t2, a->params.phi, a->params.W, tmps2, a->cache);

      // r*rho*I*cutoffs
      mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN);
      mpfr_mul(tmps[6], tmps[6], a->r, MPFR_RNDN);
      mpfr_mul(tmps[6], tmps[6], rho, MPFR_RNDN);

      // add to *out
      mpfr_add(*out, *out, tmps[6], MPFR_RNDN);
    }
  }

  return(0);
}


// R(theta)
int ss_R(mpfr_t out, mpfr_t theta, unsigned int sector, struct ss_cache cache){
  mpfr_set_ui(out, sector, MPFR_RNDN);
  mpfr_mul_ui(out, out, 8, MPFR_RNDN);
  mpfr_sub_ui(out, out, 1, MPFR_RNDN);
  mpfr_mul(out, out, cache.pi, MPFR_RNDN);
  mpfr_div_ui(out, out, 6, MPFR_RNDN);
  mpfr_add(out, out, theta, MPFR_RNDN);
  mpfr_cos(out, out, MPFR_RNDN);
  mpfr_ui_div(out, 1, out, MPFR_RNDN);
  return(0);
}

// cutoff function
int ss_cutoff(mpfr_t out, mpfr_t x, mpfr_t tmp1, mpfr_t tmp2){

  if(mpfr_cmp_d(x,0.5)<=0){
    mpfr_set_ui(out, 1, MPFR_RNDN);
    return(0);
  }
  if(mpfr_cmp_ui(x,1)>=0){
    mpfr_set_ui(out, 0, MPFR_RNDN);
    return(0);
  }

  mpfr_ui_sub(tmp1, 1, x, MPFR_RNDN);
  mpfr_d_div(tmp1, -0.5, tmp1, MPFR_RNDN);
  mpfr_exp(tmp1, tmp1, MPFR_RNDN);
  
  mpfr_sub_d(tmp2, x, 0.5, MPFR_RNDN);
  mpfr_d_div(tmp2, -0.5, tmp2, MPFR_RNDN);
  mpfr_exp(tmp2, tmp2, MPFR_RNDN);

  mpfr_ui_sub(out, 1, tmp2, MPFR_RNDN);
  mpfr_mul(out, out, tmp1, MPFR_RNDN);
  mpfr_add(out, out, tmp2, MPFR_RNDN);
  mpfr_div(out, tmp1, out, MPFR_RNDN);

  return(0);
}

// periodic norm
int ss_norm(mpfr_t out, mpfr_t k1, mpfr_t k2, struct ss_cache cache, mpfr_t tmp1, mpfr_t tmp2){
  // split R^2 into equilateral triangles and find which triangle (k1,k2) is in
  // tmp1=floor(k1/sqrt3)+1
  // tmp2=floor(k1/2/sqrt3-k2/2)+1
  // out=floor(k1/2/sqrt3+k2/2)+1
  mpfr_div(tmp1, k1, cache.sqrt3, MPFR_RNDN);
  mpfr_sub(tmp2, tmp1, k2, MPFR_RNDN);
  mpfr_add(out, tmp1, k2, MPFR_RNDN);
  mpfr_floor(tmp1, tmp1);
  mpfr_add_ui(tmp1, tmp1, 1, MPFR_RNDN);
  mpfr_div_ui(tmp2, tmp2, 2, MPFR_RNDN);
  mpfr_floor(tmp2, tmp2);
  mpfr_add_ui(tmp2, tmp2, 1, MPFR_RNDN);
  mpfr_div_ui(out, out, 2, MPFR_RNDN);
  mpfr_floor(out, out);
  mpfr_add_ui(out, out, 1, MPFR_RNDN);

  // translate (k1,k2) by -a*(SQRT3,3)-b(SQRT3,-3)
  // a-b
  // tmp2=b-a=floor((tmp2-out+1)/3);
  mpfr_sub(tmp2, tmp2, out, MPFR_RNDN);
  mpfr_add_ui(tmp2, tmp2, 1, MPFR_RNDN);
  mpfr_div_ui(tmp2, tmp2, 3, MPFR_RNDN);
  mpfr_floor(tmp2, tmp2);
  // tmp1=b=floor((tmp1+b-a)/2)
  mpfr_add(tmp1, tmp1, tmp2, MPFR_RNDN);
  mpfr_div_ui(tmp1, tmp1, 2, MPFR_RNDN);
  mpfr_floor(tmp1, tmp1);
  // tmp2=a
  mpfr_sub(tmp2, tmp1, tmp2, MPFR_RNDN);

  // out=(k1-sqrt3*(a+b))**2
  mpfr_add(out, tmp1, tmp2, MPFR_RNDN);
  mpfr_mul(out, out, cache.sqrt3, MPFR_RNDN);
  mpfr_sub(out, k1, out, MPFR_RNDN);
  mpfr_pow_ui(out, out, 2, MPFR_RNDN);
  // tmp1=(k2-3*(a-b))**2
  mpfr_sub(tmp1, tmp2, tmp1, MPFR_RNDN);
  mpfr_mul_ui(tmp1, tmp1, 3, MPFR_RNDN);
  mpfr_sub(tmp1, k2, tmp1, MPFR_RNDN);
  mpfr_pow_ui(tmp1, tmp1, 2, MPFR_RNDN);

  // out=sqrt((k1-sqrt3*(a+b))**2+(k2-3*(a-b))**2)
  mpfr_add(out, out, tmp1, MPFR_RNDN);
  mpfr_sqrt(out, out, MPFR_RNDN);
  
  return(0);
}