#include <gmp.h>
#include <mpfr.h>
#include <mpc.h>
#include "params.h"
#include "mpvec.h"
#include "poly.h"
#include "polz.h"

void mpcvec_init(mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_init2(v[i], PREC);
}

void mpcvec_clear(mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_clear(v[i]);
}

void mpcvec_frompoly(mpc_t v[N/2], const poly *a) {
  unsigned int i;

  mpcvec_init(v);

  for(i = 0; i < N/2; ++i)
    mpc_set_si_si(v[i], a->coeffs[i], a->coeffs[i + N/2], MPC_RNDNN);
}

void mpcvec_frompolz(mpc_t v[N/2], const polz *a) {
  unsigned int i;
  __mpfr_struct *real,*imag;

  mpcvec_init(v);

  for(i = 0; i < N/2; ++i) {
    real = mpc_realref(v[i]);
    imag = mpc_imagref(v[i]);
    mpc_set_ui_ui(v[i], a->high[i], a->high[i + N/2], MPC_RNDNN);
    mpc_mul_2ui(v[i], v[i], 31, MPC_RNDNN);
    mpfr_add_ui(real, real, a->mid[i], MPFR_RNDN);
    mpfr_add_ui(imag, imag, a->mid[i + N/2], MPFR_RNDN);
    mpc_mul_2ui(v[i], v[i], 31, MPC_RNDNN);
    mpfr_add_ui(real, real, a->low[i], MPFR_RNDN);
    mpfr_add_ui(imag, imag, a->low[i + N/2], MPFR_RNDN);
  }
}

void mpcvec_topolz(polz *r, const mpc_t v[N/2]) {
  unsigned int i;
  mpz_t q,tmp,rem;

  mpz_init(q);
  mpz_init(tmp);
  mpz_init(rem);

  //FIXME: Barrett reduction!
  mpz_set_ui(q,QH);
  mpz_mul_2exp(q,q,31);
  mpz_add_ui(q,q,QM);
  mpz_mul_2exp(q,q,31);
  mpz_add_ui(q,q,QL);

  for(i = 0; i < N/2; ++i) {
    mpfr_get_z(tmp,mpc_realref(v[i]),MPFR_RNDN);
    mpz_tdiv_r(tmp,tmp,q);
    mpz_fdiv_r_2exp(rem,tmp,31);
    mpz_fdiv_q_2exp(tmp,tmp,31);
    r->low[i] = mpz_get_ui(rem);
    mpz_fdiv_r_2exp(rem,tmp,31);
    mpz_fdiv_q_2exp(tmp,tmp,31);
    r->mid[i] = mpz_get_ui(rem);
    r->high[i] = mpz_get_si(tmp);

    mpfr_get_z(tmp,mpc_imagref(v[i]),MPFR_RNDN);
    mpz_tdiv_r(tmp,tmp,q);
    mpz_fdiv_r_2exp(rem,tmp,31);
    mpz_fdiv_q_2exp(tmp,tmp,31);
    r->low[i+N/2] = mpz_get_ui(rem);
    mpz_fdiv_r_2exp(rem,tmp,31);
    mpz_fdiv_q_2exp(tmp,tmp,31);
    r->mid[i+N/2] = mpz_get_ui(rem);
    r->high[i+N/2] = mpz_get_si(tmp);
  }

  mpz_clear(q);
  mpz_clear(tmp);
  mpz_clear(rem);
}

void mpcvec_set(mpc_t u[N/2], const mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_set(u[i], v[i], MPC_RNDNN);
}

void mpcvec_set_ui(mpc_t u[N/2], unsigned long a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_set_ui(u[i], a, MPC_RNDNN);
}

void mpcvec_set_si(mpc_t u[N/2], long a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_set_si(u[i], a, MPC_RNDNN);
}

void mpcvec_set_fr(mpc_t u[N/2], const mpfr_t a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_set_fr(u[i], a, MPC_RNDNN);
}

void mpcvec_abssq(mpfr_t u[N/2], const mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_norm(u[i], v[i], MPC_RNDNN);
}

void mpcvec_conj(mpc_t u[N/2], const mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_conj(u[i], v[i], MPC_RNDNN);
}

void mpcvec_add(mpc_t w[N/2], const mpc_t u[N/2],  const mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_add(w[i], u[i], v[i], MPC_RNDNN);
}

void mpcvec_add_ui(mpc_t u[N/2], const mpc_t v[N/2], unsigned long a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_add_ui(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_add_fr(mpc_t u[N/2], const mpc_t v[N/2],  const mpfr_t a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_add_fr(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_sub(mpc_t w[N/2], const mpc_t u[N/2],  const mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_sub(w[i], u[i], v[i], MPC_RNDNN);
}

void mpcvec_sub_ui(mpc_t u[N/2], const mpc_t v[N/2], unsigned long a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_sub_ui(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_sub_fr(mpc_t u[N/2], const mpc_t v[N/2],  const mpfr_t a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_sub_fr(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_mul(mpc_t w[N/2], const mpc_t u[N/2],  const mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_mul(w[i], u[i], v[i], MPC_RNDNN);
}

void mpcvec_mul_ui(mpc_t u[N/2], const mpc_t v[N/2], unsigned long a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_mul_ui(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_mul_si(mpc_t u[N/2], const mpc_t v[N/2], long a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_mul_si(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_mul_fr(mpc_t u[N/2], const mpc_t v[N/2],  const mpfr_t a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_mul_fr(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_mul_frvec(mpc_t w[N/2], const mpc_t u[N/2],  const mpfr_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_mul_fr(w[i], u[i], v[i], MPC_RNDNN);
}

void mpcvec_fma(mpc_t w[N/2],
                const mpc_t t[N/2],
                const mpc_t u[N/2],
                const mpc_t v[N/2])
{
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_fma(w[i], t[i], u[i], v[i], MPC_RNDNN);
}

void mpcvec_div(mpc_t w[N/2], const mpc_t u[N/2],  const mpc_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_div(w[i], u[i], v[i], MPC_RNDNN);
}

void mpcvec_div_ui(mpc_t u[N/2], const mpc_t v[N/2], unsigned long a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_div_ui(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_div_fr(mpc_t u[N/2], const mpc_t v[N/2], const mpfr_t a) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_div_fr(u[i], v[i], a, MPC_RNDNN);
}

void mpcvec_div_frvec(mpc_t w[N/2], const mpc_t u[N/2],  const mpfr_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpc_div_fr(w[i], u[i], v[i], MPC_RNDNN);
}

void mpfrvec_init(mpfr_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpfr_init2(v[i], PREC);
}

void mpfrvec_clear(mpfr_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpfr_clear(v[i]);
}

void mpfrvec_add(mpfr_t w[N/2], const mpfr_t u[N/2],  const mpfr_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpfr_add(w[i], u[i], v[i], MPC_RNDNN);
}

void mpfrvec_sub(mpfr_t w[N/2], const mpfr_t u[N/2],  const mpfr_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpfr_sub(w[i], u[i], v[i], MPC_RNDNN);
}

void mpfrvec_mul(mpfr_t w[N/2], const mpfr_t u[N/2],  const mpfr_t v[N/2]) {
  unsigned int i;

  for(i = 0; i < N/2; ++i)
    mpfr_mul(w[i], u[i], v[i], MPC_RNDNN);
}
