/* Given:

        1) some "small" characteristic exponent h>0, and
        2) (optionally) the number of significative decimals Q>=0,

   we compute the Lazutkin homoclinic invariant lhi of the primary
   homoclinic orbit on the symmetry line {y=0} of the Henon map

          (x,y) -> ( x+y+epsilon*x(1-x), y+epsilon*x*(1-x) )

   up to Q significative decimals. (By default, Q=50.)

   Here, epsilon = 4*sinh^2(h/2) = lambda-2+1/lambda and lambda=exp(h).

   The symmetric homoclinic orbit is found by a Newton's method
   which doubles the precision after each iteration. 

  Variables:

   * lambda=exp(h) is the characteristic multiplier
   * S=2*Pi^2/(h*log(10)) is the number of decimal digits lost by cancellation
   * P_single is the "single" precision of the multiple precision arithmetic
   * P is the "final" precision of the multiple precision arithmetic
   * p is the "variable" precision of the multiple precision arithmetic
   * K is the order of the Taylor expansions
   * N is the number of iterations to reach the symmetry line from the
       fundamental domain
   * [r1/lambda, r1), r1>0, is the fundamental interval
   * m(r)=(X(r),Y(r)) is the natural parameterization of the unstable curve
   * { (x,y) : G(x,y)=y=0 } is the symmetry line
   * Z(r)=G(f^N(m(r))=0 is the one-dimensional equation for the homoclinic
                        point
   * r0 is the solution of the above equation in the fundamental interval
   * [x,y] are the coordinates of the point on which acts the Henon map
   * [dx,dy] is the tangent vector on which acts the differetial of the Henon
             map, or the tangent vector to the unstable curve at the homoclinic
   * [dxs,dys]=[dx-dy,-dy] is the tangent vector to the stable curve
   * lhi is the Lazutkin homoclinic invariant

  User-defined functions (in order of appearance):

   * Taylor() computes the Taylor expansion of m(r)=(X(r),Y(r)) up to order K
   * f() iterates the Henon map
   * df() iterates the differential of the Henon map
   * fN() iterates N times the Henon map
   * dfN() iterates N times the differential of the Henon map
   * XY(r) evaluates m(r)=(X(r),Y(r)) by Horner's rule
   * dXY(r) evaluates m(r)=(X(r),Y(r)) and its derivative by Horner's rule
   * Z(r) evaluates the function Z(r)=G(f^N(m(r))
   * Z_Newton(r) evaluates the "Newton step" Z(r)/Z'(r)

  Some GP-instructions:

   * vector(N,n,expr(n)) creates the vector [expr(1),...,expr(N)]
   * if(cond,"seq1","seq2") executes "seq1" if cond != 0, else "seq2"
   * for(k=a,b,"seq") executes "seq" for k=a,a+1,...,b-1,b
   * while(cond,"seq") executes "seq" while cond != 0
   * sum(k=a,b,"expr") sums "expr" from k=a to k=b
   * default(realprecision,p) sets the precision equal to p
   * default(format,"val") sets the exit format equal to "val"
   * solve(r=a,b,"expr") finds the root of "expr" between a and b
   * x\/y = ceil(x/y)
   * x%y = remainder of the Euclidean division x/y

  To run this program: 1) Open a GP session, 2) Read this file: \r Henon.gp,
   3) Call the main function: LHI_Henon(h,Q), for some real h and integer Q. */
\\ Main function
LHI_Henon(h,Q=50)=
{
  local(K,N,p,P,P_single,S);
  local(epsilon,lambda);
  local(X,Y);
  local(x,y,dx,dy,dxs,dys);
  local(lhi,r0,r1);
  local(time_Taylor,time_homoclinic,time_LHI);

  P_single=28;
  S=ceil( 2*Pi*Pi/(h*log(10)) );
  P=max(ceil(1.1*(Q+S)),P_single);

  h=precision(h,P);
  default(realprecision,P);
  lambda=exp(h);
  epsilon=lambda-2+1/lambda;

  default(format,"e0.10");
  print("The characteristic exponent is ",h,".");
  print1("The multiple-precision works with P=1.1(Q+S)=1.1(");
  print(Q,"+",S,")=",P," decimals.");

  K=max(20, ceil((14*P*log(10)/h)^(1/3)) );
  print1("Computing the Taylor expansions up to order K=",K,"... ");
  gettime();
  X=vector(K,k,0.);
  Y=vector(K,k,0.);
  Taylor();
  time_Taylor=gettime();
  print("Done in ",time_Taylor," ms.");

  print1("Computing the homoclinic orbit ");
  default(realprecision,P_single);
  r1=5*10^(-P/(K+1));
  XY(r1/lambda);
  N=-1;
  while(y>0,N++;f());
  print1("using N=",N," iterations... ");
  r0=solve(r=r1/lambda,r1,Z(r));
  p=P;
  while(p>=P_single-3,p=p\/2);
  while( p<P,
         p=min(2*p,P);
         default(realprecision,p);
         r0=precision(r0,p);
         r0-=Z_Newton(r0);
       );
  time_homoclinic=gettime();
  print("Done in ",time_homoclinic," ms.");

  print1("Computing the Lazutkin homoclinic invariant... ");
  dXY(r0);
  dfN();
  dxs=dx-dy;
  dys=-dy;
  lhi=r0^2*matdet( [dxs,dx; dys,dy] );
  lhi=precision(lhi,Q);
  time_LHI=gettime();
  print("Done in ",time_LHI," ms.");

  print("Total CPU time = ",time_Taylor+time_homoclinic+time_LHI," ms.");
  default(format,"e0.100");
  print("LHI = ",lhi);
}

\\ Computation of the Taylor expansion of m(r)=(X(r),Y(r)) up to order K
Taylor()=
{
  local(l,aux0,aux1,aux2);

  /* Store in vector l the powers of the characteristic multiplier lambda */
  l=vector(K,k,lambda);
  for(k=2,K,l[k]=lambda*l[k-1]);

  /* Normalize & Initialize */
  X[1]=1.;
  Y[1]=(lambda-1)*X[1]/lambda;

  /* The recursive computations */
  for(k=2,K,
      aux0=if( k%2, 2*sum(j=1,(k-1)/2,X[j]*X[k-j]),
                    2*sum(j=1,k/2-1,X[j]*X[k-j]) + X[k/2]^2 );
      aux1=(l[k]-1)/l[k];
      aux2=(l[k]-1)*aux1;
      X[k]=epsilon*aux0/(epsilon-aux2);
      Y[k]=aux1*X[k];
     );
}

\\ Henon map f(x,y)=(x',y')
f()=
{
  local(aux);

  aux=epsilon*x*(1-x);
  y=y+aux;
  x=x+y;
}

\\ Differential of the Henon map: f(x,y)=(x',y') and df(x,y)[dx,dy]=[dx',dy']
df()=
{
  local(aux0,aux1,aux2);

  aux0=epsilon*x;
  aux1=aux0*(1-x);
  aux2=(epsilon-2*aux0)*dx;
  y=y+aux1;
  dy=dy+aux2;
  x=x+y;
  dx=dx+dy;
}

\\ N-th iterate of the Henon map f
fN()=for(n=1,N,f());

\\ N-th iterate of the Differential of the Henon map f
dfN()=for(n=1,N,df());

\\ Evaluation of the natural parameterization m(r)=(X(r),Y(r)) by Horner's rule
XY(r)=
{
  x=X[K]; y=Y[K];
  for(k=1,K-1,
      x=x*r+X[K-k];
      y=y*r+Y[K-k];
     );
  x=x*r; y=y*r;
}

\\ Evaluation of m(r)=(X(r),Y(r)) and its derivative by Horner's rule
dXY(r)=
{
  x=X[K];
  y=Y[K];
  dx=K*X[K];
  dy=K*Y[K];
  for(k=1,K-1,
      x=x*r+X[K-k];
      y=y*r+Y[K-k];
      dx=dx*r+(K-k)*X[K-k];
      dy=dy*r+(K-k)*Y[K-k];
     );
  x=x*r;
  y=y*r;
}

\\ Evaluation of the function Z(r)=Z_N(r)=G(f^N(m(r))), where G(x,y)=y
Z(r)=
{
  XY(r);
  fN();
  return(y);
}

\\ Evaluation of Z(r)/Z'(r), where the function Z(r) is defined above
Z_Newton(r)=
{
  dXY(r);
  dfN();
  return(y/dy);
}


/********************  WRITTEN  BY  RRR (2007, 2010)   ***********************/

/**********  AVAILABLE AT www.ma1.upc.edu/~rafael/research.html  *************/

/*******************  CONTACT:  Rafael.Ramirez@upc.edu  **********************/

/************************   THAT'S  ALL,  FOLKS!   ***************************/