/*

  Tools for integration on H_p x H

*/

#include "pari.h"

extern GEN ellcuspint(GEN,GEN,GEN);
extern GEN cpsub(GEN,GEN);
extern GEN cpmul(GEN,GEN);
extern GEN cpdiv(GEN,GEN);
extern GEN cppow(GEN,long);
GEN hpxhint(GEN,GEN,GEN,GEN,GEN,GEN);
GEN hpxhballs(GEN,GEN,GEN);

GEN hpxhint(GEN tau1,GEN tau2,GEN z1,GEN z2,GEN table,GEN prec) {
  long ltop2,ltop = avma;
  GEN a1,a2,b1,b2,p,s,balls,ballsk,PROD,t,val,infty,ptoradius,tsqr,gammaz1,gammaz2,absgamma,temp;
  long k,j,precs,lim,numballs,mu,musum,radius,rebound;

  /* collect info on tau1 and tau2 */
  if (!(typ(tau1) == t_VEC && lg(tau1) == 4 && typ(tau2) == t_VEC && lg(tau2) == 4)) err(talker,"tau1 and tau2 must be in H_p");
  a1 = (GEN) tau1[1]; b1 = (GEN) tau1[2];
  a2 = (GEN) tau2[1]; b2 = (GEN) tau2[2];
  if ((typ(b1) != t_PADIC) || (typ(b2) != t_PADIC)) err(talker,"tau1 and tau2 must be in H_p with explicit precision");
  p = (GEN) b1[2];
  if (cmpii(p,(GEN) b2[2]) != 0) err(talker,"tau1 and tau2 must be in H_p for a fixed p");

  precs = itos(prec);

  /* prepare a list of balls */
  balls = cgetg(itos(gadd(gpowgs(p,precs),gpowgs(p,precs-1)))+1,t_VEC);
  for(j = 1; j < lg(balls); j++) {
    balls[j] = zero;
  }
  k = 1;
  
  lim = itos(gsub(gpowgs(p,precs),gun));
  for (j = 0; j <= lim; j++) {
    ballsk = cgetg(3,t_VEC);
    ballsk[1] = (long) stoi(j);
    ballsk[2] = (long) prec;
    balls[k] = (long) ballsk;
    k++;
  }
  
  ballsk = cgetg(3,t_VEC);
  ballsk[1] = (long) ginv(gpowgs(p,precs+2));
  ballsk[2] = (long) prec;
  balls[k] = (long) ballsk;
  k++;
  
  lim = itos(gsub(gpowgs(p,precs-1),gun));
  for (j = 1; j <= lim; j++) {
    ballsk = cgetg(3,t_VEC);
    ballsk[1] = (long) ginv(gmul(p,stoi(j)));
    ballsk[2] = (long) prec;
    balls[k] = (long) ballsk;
    k++;
  }

  /* main loop */
  PROD = gun;
  musum = 0;
  numballs = lg(balls) - 1;
  infty = cgetg(2,t_VEC); infty[1] = un;
  ltop2 = avma; rebound = 1;

  for (j = 1; j <= numballs; j++) {

    t = (GEN) ((GEN) balls[j])[1]; radius = itos((GEN) ((GEN) balls[j])[2]);
    ptoradius = gpowgs(p,radius);
    tsqr = gsqr(t);

    /* evaluate mu{z1,z2}(U_j) */
    if (ggval(t,p) >= 0) {
      gammaz1 = gdiv(gsub(z1,t),ptoradius);
      gammaz2 = gdiv(gsub(z2,t),ptoradius);
    } else {
      if (gcmp(z1,gzero) == 0) gammaz1 = infty;
      else gammaz1 = gdiv(gsub(ginv(t),ginv(z1)),ptoradius);
      if (gcmp(z2,gzero) == 0) gammaz2 = infty;
      else gammaz2 = gdiv(gsub(ginv(t),ginv(z2)),ptoradius);
    }
    absgamma = gun; /* THIS NEEDS TO BE CORRECTED! */
    temp = gmul(absgamma,ellcuspint(table,gammaz1,gammaz2));
    mu = itos((GEN) temp[1]) + itos((GEN) temp[2]);

    musum = musum + mu;

    /*
      printf("[%i t=",j);
      outbrute(t);
      printf(", gz1=");
      outbrute(gammaz1);
      printf(", gz2=");
      outbrute(gammaz2);
      printf(", radius=%i, mu=%i, musum=%i]\n",radius,mu,musum);
    */

    if (mu != 0) {

      val = cpdiv(cpsub(t,tau2),cpsub(t,tau1));
      
      PROD = cpmul(PROD,cppow(val,mu));
      
      if (rebound == 1000) {
	PROD = gerepileupto(ltop2,PROD);
	rebound = 0;

      }
      
      rebound++;
    }
  }

  printf("musum : %i\n",musum);

  return gerepileupto(ltop,PROD);
}


GEN hpxhballs(GEN table, GEN p, GEN prec) {
  long ltop = avma;
  long ltop2;  
  GEN balls,ballsk,infty,oneoverp,z1,z2,t,ptoradius,gammaz1,gammaz2,absgamma,temp;
  long k,j,lim,numballs,radius,mu,musum;
  long precs = itos(prec);

  /* a little setup */
  infty = cgetg(2,t_VEC); infty[1] = un;
  oneoverp = ginv(p);

  /* prepare a list of balls, with space allocated for measures */
  balls = cgetg(itos(gadd(gpowgs(p,precs),gpowgs(p,precs-1)))+1,t_VEC);
  for(j = 1; j < lg(balls); j++) {
    balls[j] = zero;
  }
  k = 1;

  lim = itos(gsub(gpowgs(p,precs),gun));
  for (j = 0; j <= lim; j++) {
    ballsk = cgetg(3,t_VEC);
    ballsk[1] = (long) stoi(j);
    /* ballsk[2] = lcopy(prec); */
    balls[k] = (long) ballsk;
    k++;
  }
  
  ballsk = cgetg(3,t_VEC);
  ballsk[1] = (long) gpowgs(p,-precs-2);
  /* ballsk[2] = lcopy(prec); */
  balls[k] = (long) ballsk;
  k++;

  lim = itos(gsub(gpowgs(p,precs-1),gun));
  for (j = 1; j <= lim; j++) {
    ballsk = cgetg(3,t_VEC);
    ballsk[1] = (long) gdiv(oneoverp,stoi(j));
    /* ballsk[2] = lcopy(prec); */
    balls[k] = (long) ballsk;
    k++;
  }

  /* main loop */
  musum = 0;
  numballs = lg(balls) - 1;
  z1 = gzero; z2 = infty;

  for (j = 1; j <= numballs; j++) {

    ltop2 = avma;

    t = (GEN) ((GEN) balls[j])[1]; radius = precs /* itos((GEN) ((GEN) balls[j])[2])*/;
    ptoradius = gpowgs(p,radius);

    /* evaluate mu{z1,z2}(U_j) */
    if (ggval(t,p) >= 0) {

      if (typ(z1) == t_VEC) gammaz1 = infty;
      else gammaz1 = gdiv(gsub(z1,t),ptoradius);

      if (typ(z2) == t_VEC) gammaz2 = infty;
      else gammaz2 = gdiv(gsub(z2,t),ptoradius);

    } else {

      if (typ(z1) == t_VEC) gammaz1 = gdiv(ginv(t),ptoradius);
      else if (gcmp0(z1)) gammaz1 = infty;
      else gammaz1 = gdiv(gsub(ginv(t),ginv(z1)),ptoradius);

      if (typ(z2) == t_VEC) gammaz2 = gdiv(ginv(t),ptoradius);
      else if (gcmp0(z2)) gammaz2 = infty;
      else gammaz2 = gdiv(gsub(ginv(t),ginv(z2)),ptoradius);

    }

    absgamma = gun; /* THIS NEEDS TO BE CORRECTED! */

    temp = gmul(absgamma,ellcuspint(table,gammaz1,gammaz2));
    mu = itos((GEN) temp[1]) + itos((GEN) temp[2]);

    musum = musum + mu;

    avma = ltop2;

    ((GEN) balls[j])[2] = (long) stoi(mu);

  }

  printf("musum=%i\n",musum);

  return balls;

}