/*

  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 hpxhint9(GEN,GEN,GEN,GEN,GEN,GEN);
GEN hpxhmutable(GEN,GEN,GEN);
void hpxhshowarray(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);
}


#define evalball {\
    ptoradius = gpowgs(p,radius);\
    tsqr = gsqr(t);\
    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;\
    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++;\
    }\
  }


GEN hpxhint9(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);

  PROD = gun;
  musum = 0;
  infty = cgetg(2,t_VEC); infty[1] = un;
  ltop2 = avma; rebound = 1;

  /* balls */
  lim = itos(gsub(gpowgs(p,precs),gun));
  for (j = 0; j <= lim; j++) {
    t = stoi(j); radius = precs;
    evalball;
  }
  
  t = ginv(gpowgs(p,precs+2)); radius = precs;
  evalball;
  
  lim = itos(gsub(gpowgs(p,precs-1),gun));
  for (j = 1; j <= lim; j++) {
    t = ginv(gmul(p,stoi(j))); radius = precs;
    evalball;
  }

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

  return gerepileupto(ltop,PROD);
}



GEN hpxhmutable(GEN manintable, GEN Gp, GEN Gprec) {
  long ltop = avma;
  long ltop2;
  long p = itos(Gp);
  long prec = itos(Gprec);
  GEN Gptoprecminusone = gpowgs(Gp,prec-1);
  GEN Gptoprec = gmul(Gptoprecminusone,Gp);
  long ptoprecminusone = itos(Gptoprecminusone);
  long ptoprec = itos(Gptoprec);
  long uplength = ptoprec;
  long downlength = ptoprecminusone;
  long musum = 0;
  int *uparray,*downarray;
  long j;
  GEN infty,absgamma,gammaz1,gammaz2,temp,outvec;

  uparray = malloc(uplength*sizeof(int));
  downarray = malloc(downlength*sizeof(int));

  infty = cgetg(2,t_VEC); infty[1] = un;
  absgamma = gun; /* THIS NEEDS TO BE CHANGED !!! */

  /* we will assume we want just mu{0,infty}(B((j:1),1/p^prec)) */
  gammaz2 = infty;

  for (j = 0; j < uplength; j++) {
    ltop2 = avma;
    gammaz1 = gdiv(stoi(-j),Gptoprec);
    temp = gmul(absgamma,ellcuspint(manintable,gammaz1,gammaz2));
    uparray[j] = itos((GEN) temp[1]) + itos((GEN) temp[2]);
    musum = musum + uparray[j];
    avma = ltop2;
  }

  /* we will assume we want just mu{0,infty}(B((1:p*j),1/p^prec)) */
  gammaz1 = infty;

  for (j = 0; j < downlength; j++) {
    ltop2 = avma;
    gammaz2 = gdiv(stoi(j),Gptoprecminusone);
    temp = gmul(absgamma,ellcuspint(manintable,gammaz1,gammaz2));
    downarray[j] = itos((GEN) temp[1]) + itos((GEN) temp[2]);
    musum = musum + downarray[j];
    avma = ltop2;
  }

  avma = ltop;

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

  outvec = cgetg(3,t_VEC);
  outvec[1] = (long) stoi((long) uparray);
  outvec[2] = (long) stoi((long) downarray);

  return outvec;
}

void hpxhshowarray(GEN address, GEN start, GEN end) {
  long *intarray = (long *) itos(address);
  long val1 = itos(start);
  long val2 = itos(end);
  long j;

  printf("[ ");
  for (j = val1; j < val2; j++) {
    printf("%i, ",intarray[j]);
  }
  printf("%i ]\n",intarray[val2]);
}

GEN hpxheval(GEN mus, GEN taus) {
  long ltop = avma;
  long i,j,k,up,mu,*muarray,start,end,taupow;
  GEN Gj,mui,tau1,tau2,vl,SUBPROD,PROD;
  long cycles;

  PROD = gun;

  for (i = 1; i < lg(mus); i++) {
    mui = (GEN) mus[i];
    up = itos((GEN) mui[1]);
    muarray = (long *) itos((GEN) mui[2]);
    start = itos((GEN) mui[3]);
    end = itos((GEN) mui[4]);
    cycles = 0;

    if (up) {

      for (j = start; j <= end; j++) {
	SUBPROD = gun;
	mu = muarray[j];
	Gj = stoi(j);

	if (mu != 0) for(k = 1; k < lg(taus); k++) {
	  tau1 = (GEN) ((GEN) taus[k])[1];
	  tau2 = (GEN) ((GEN) taus[k])[2];
	  taupow = itos((GEN) ((GEN) taus[k])[3]);
	  vl = cpdiv(cpsub(Gj,tau2),cpsub(Gj,tau1));
	  SUBPROD = cpmul(cppow(vl,taupow),SUBPROD);
	}  

	PROD = cpmul(cppow(SUBPROD,mu),PROD);

	if (cycles == 1000) { cycles = 0; PROD = gerepileupto(ltop,PROD); }
	cycles++;

      }

    } else {

      for (j = start; j <= end; j++) {
	SUBPROD = gun;
	mu = muarray[j];

	if (mu != 0) for(k = 1; k < lg(taus); k++) {
	  tau1 = (GEN) ((GEN) taus[k])[1];
	  tau2 = (GEN) ((GEN) taus[k])[2];
	  taupow = itos((GEN) ((GEN) taus[k])[3]);
	  vl = cpdiv(cpsub(gun,cpmul(Gj,tau2)),cpsub(gun,cpmul(Gj,tau1)));
	  SUBPROD = cpmul(cppow(vl,taupow),SUBPROD);
	}  

	PROD = cpmul(cppow(SUBPROD,mu),PROD);

	if (cycles == 1000) { cycles = 0; PROD = gerepileupto(ltop,PROD); }
	cycles++;

      }

      PROD = gerepileupto(ltop,PROD);

    }

  }

  return PROD;

}