Routines d'intégration numérique prouvée

Description

Ces fonctions correspondant aux théorèmes et paramètres démontrés dans ma thèse.

Fonctions utiles

/* lambert function Lx */
Lx(y)= {
  my(eps,x,tmp,yexp);
  if(y<exp(1),print("*** Lx : argument lower than e *");return(0););
  eps = 10.^(-default(realprecision)+2);
  /* méthode de Newton */
  x=log(y);
  tmp=0;
  while(abs(x-tmp)>eps,
      tmp = x;
      /* f(x)=x-log(xy), f'(x)=1-1/x */
      x*=(log(x*y)-1)/(x-1);
      );
  return(x);
};
/* upper bound for X=Lαβ(D) s.t.
   \int_X^\infty t^l*e^{-αt^β} \dt < e^{-D}
 */
Lab(D,a=1,b=1,l=1) = {
  \\ solve X^(l+1)e^{-aX^b} = e^{-D}
  \\ ie e^U/U = (l+1)/(ab)e^{bD/(l+1)}, U=ab/(l+1)*X^b
    l++;
  return(l/(a*b)*Lx(l/(a*b)*exp(b*D/l)))^(1/b);
};
/* print informations */
msgdebug(s,header="\t",level=1) = {
  if(default(debug)>=level,printf("%s%s\n",header,s))
};
DBG=1; EXTDBG=3; HEADER="  [intnum] ";
msgprec(e) = {
  if(e,msgdebug(Strprintf("uncertain digits = %i",
          ceil(log(e)/log(10))),HEADER,DBG))
};
msgpoints(npoints,fct) = {
  msgdebug(Strprintf("npoints = %i [%s]",npoints,fct),,DBG)
};

Cadre d'une décroissance DE

méthode des trapèzes avec pas d'intégration donné

f la fonction à intégrer
h le pas d'intégration
npoints le nombre de points
la somme des trapèzes h * sum f(kh)

integration(f,h,npoints) = {
  my(kh,res,a,b,supres);
  kh=0;
  res=f(0);
  supres=0;
  for(k=1,npoints,
      kh+=h;
      \\a = f(-kh);
      \\b = f(kh);
      \\res += a + b
      res+=f(kh)+f(-kh);
      supres=max(abs(res),supres);
     );
  msgprec(abs(supres*npoints/res));
  res*=h;
  return(res);
};
addhelp(integration,"integration(f,h,npoints)");
integration_until(f,h,D=0,moreprec=2) = {
  my(kh,res,tmp,supres,npoints);
  if(!D,D=precision(1.);moreprec=0);
  D = (D+moreprec)*log(10);
  kh=0;
  res=f(0);
  supres=0;
  while(1,
      kh+=h;
      tmp = f(-kh);
      tmp += f(kh);
      res += tmp;
      supres=max(abs(res),supres);
      if(tmp&&log(abs(tmp/res))+D+moreprec<0,
        npoints = round(kh/h);
        msgpoints(npoints,);
        break());
      );
  msgprec(abs(supres*npoints/res));
  res*=h;
  return(res);
};

paramètres pour la méthode des trapèzes

g holomorphe sur Δ_τ
g(±x) ≤ M1 e^(-αe^(βx))
g(z) ≤ M2 e^{λx+Ae^(γx)} sur ±x±it

longueur d'intégration

tmax_DE(D,M1,alpha,beta) = {
  return(log(D+log(4*M1/(alpha*beta)))/beta);
};
addhelp(tmax_DE,"tmax_DE(D,M1,alpha,beta)");

hypothèses du théorème DE 1

norme = \int_{\R+iτ} \abs{g}+\int_{\R-iτ} \abs{g}

pas_h_DE1(D,tau,norme) = {
  my(h);
  if(!D,D=default(realprecision)*log(10));
  h = 2*Pi*tau/(D+log(4*norme+2*exp(-D)));
  msgdebug( Strprintf("h=%1.3e [pas_h_DE1]",h),,DBG);
  return(h);
};
addhelp(pas_h_DE1,"pas_h_DE1(D,tau,norme)");

hypothèses du théorème DE 2

g = O(e^{-αe^β})
g ≤ M2 sur Δ_τ

pas_h_DE2(D,tau,M2=1,alpha=1,beta=1) = {
  my(h,at,bt);
  if(!D,D=default(realprecision)*log(10));
  if(beta*tau<Pi/2,delta=1/tan(beta*tau),delta=0);
  bt = atan((D-beta*delta*tau)/(delta*D+beta*tau));
  at = alpha*(cos(bt)-delta*sin(bt));
  h = 2*Pi*bt/(beta*(D+log(4*M2*exp(-at)/(at*beta)+2*exp(-D))));
  msgdebug( Strprintf("h=%1.3e [pas_h_DE2]",h),,DBG);
  return(h);
};
addhelp(pas_h_DE2,"pas_h_DE2(D,tau,M2,alpha,beta)");

hypothèses du théorème DE 3

g = O(e^{-αe^β})
g ≤ M2e^{λz+Ae^γz} sur Δ_τ

/* premier cas : A = 0 */
pas_h_DE3_noA(D,tau,M2=1,alpha=1,beta=1,lambda=0) = {
  my(Dtilde,a,t,delta,at,Ct,l);
  if(!D,D=default(realprecision)*log(10));
  if(lambda=0,return(pas_h_DE2(D,tau,M2,alpha,beta)));
  a = alpha*beta*(sin(beta*tau)+1/sin(beta*tau));
  l = (lambda+1)/beta;
  Dtilde = (D+log(4*M2+2*exp(-D)))/l+log(l/a)-1;
  /* ne pas descendre en dessous de tau/2 */
  if(tau/2*log(tau/2)+tau-tau/2*(Dtilde+1)>=0,
      t = tau/2,
      \\t = tau-solve(x=exp(-D),tau/2,x*log(x)+tau-x*(Dtilde+1))
      t = tau*(1-1/Lx(exp(Dtilde+1)/tau));
    );
  if(beta*tau>=Pi/2,delta=0,delta=1/tan(beta*tau));
  at = alpha*(cos(beta*t)-delta*sin(beta*t));
  Ct = l*(log(l/at)-1);
  msgdebug( Strprintf("t=tau-%1.2e\n\tCt=%1.2e",tau-t,Ct),,DBG);
  return(pas_h_DE1(D+Ct,t,M2));
};
addhelp(pas_h_DE3_noA,"pas_h_DE3_noA(D,tau,M2,alpha,beta,lambda)");
pas_h_DE3(D,tau,M2=1,alpha=1,beta=1,lambda=0,A=0,gamma_=0) = {
  my(a1,l,U,Dtilde,X,t,at,At,Ct);
  if(!D,D=default(realprecision)*log(10));
  if(A<=0,return(pas_h_DE3_noA(D,tau,M2,alpha,beta,lambda)));
  /* at ~ a1(tau-t) = Ux */
  a1 = alpha*beta*(sin(beta*tau)+1/sin(beta*tau));
  /* Ct ~ U x^(-l) */
  l =  gamma_/(beta-gamma_);
  U = A*((A*gamma_+lambda+1)/(beta*a1))^l-a1*((A*gamma_)/(beta*a1))^(l+1);
  Dtilde = D+log(4*M2+2*exp(-D));
  /* t = tau-X */
  X = solve(x=0,tau,Dtilde*x^(l+1)+U*x-l*U*(tau-x));
  t = max(tau-X,tau/2);
  if(beta*tau>=Pi/2,delta=0,delta=1/tan(beta*tau));
  at = alpha*(cos(beta*t)-delta*sin(beta*t));
  At = A*cos(gamma_*t)/cos(gamma_*tau);
  num = At*gamma_+lambda+1; den = at*beta;
  Ct = At*(num/den)^l -at*(At*gamma_/den)^(l+1)
    +(lambda+1)/(beta-gamma_)*log(num/den);
  msgdebug( Strprintf("\tt=tau-%1.2e\n\tCt=%1.2e",X,Ct),,DBG);
  return(pas_h_DE1(D+Ct,t,M2));
};
addhelp(pas_h_DE3,"pas_h_DE3(D,tau,M2,alpha,beta,lambda,A,gamma)");

routines d'intégration DE

integration_DE1(f,M1,alpha,beta,tau,M2,D) = {
  my(h,npoints);
  if(!D,D=default(realprecision)*log(10));
  h = pas_h_DE1(D,tau,M2);
  npoints = ceil(tmax_DE(D,M1,alpha,beta)/h);
  msgdebug( Strprintf("\tnpoints=%i\n",npoints),,DBG);
  return(integration(f,h,npoints));
};
addhelp(integration_DE1,"integration_DE1(f,M1,alpha,beta,tau,M2,D)");
integration_DE2(f,M1,alpha,beta,tau,M2,D) = {
  my(h,npoints);
  if(!D,D=default(realprecision)*log(10));
  h = pas_h_DE2(D,tau,M2,alpha,beta);
  npoints = ceil(tmax_DE(D,M1,alpha,beta)/h);
  msgdebug( Strprintf("\tnpoints=%i\n",npoints),,DBG);
  return(integration(f,h,npoints));
};
addhelp(integration_DE2,"integration_DE2(f,M1,alpha,beta,tau,M2,D)");
integration_DE3(f,M1,alpha,beta,tau,M2,lambda,A,gamma_,D) = {
  my(h,npoints);
  if(!D,D=default(realprecision)*log(10));
  h = pas_h_DE3(D,tau,M2,alpha,beta,lambda,A,gamma_);
  npoints = ceil(tmax_DE(D,M1,alpha,beta)/h);
  msgdebug( Strprintf("\tnpoints=%i\n",npoints),,DBG);
  return(integration(f,h,npoints));
};
addhelp(integration_DE3,"integration_DE3(f,M1,alpha,beta,tau,M2,lambda,A,gamma_,D)");

Changements de variable

sur R avec décroissance simplement exponentielle

integration_sh(f,h,npoints) = {
  my(res,exph,i_exph,ekh,i_ekh,sh,ch,supres);
  res=f(0);
  exph=exp(h);i_exph=1/exph;
  ekh=i_ekh=1;
  supres=0;
  for(k=1,npoints,
      ekh*=exph;
      i_ekh*=i_exph;
      sh=(ekh-i_ekh)/2;
      ch=(ekh+i_ekh)/2;
      res+=(f(sh)+f(-sh))*ch;
     );
  msgprec(abs(supres*npoints/res));
  res*=h;
  return(res);
};
addhelp(integration_sh,"integration_sh(f,h,npoints)");
integration_sh_until(f,h,D,moreprec=2) = {
  my(res,exph,i_exph,ekh,i_ekh,sh,ch,tmp,supres,npoints);
  if(!D,D=precision(1.));
  D = (D+moreprec)*log(10);
  res=f(0);
  exph=exp(h);i_exph=1/exph;
  ekh=i_ekh=1;
  supres=0;
  while(1,
      ekh*=exph;
      i_ekh*=i_exph;
      sh=(ekh-i_ekh)/2;
      ch=(ekh+i_ekh)/2;
      tmp=(f(sh)+f(-sh))*ch;
      res+=tmp;
      supres=max(abs(res),supres);
      if(tmp&&log(abs(tmp/res))+D+moreprec<0,
        npoints = round(log(ekh)/log(exph));
        msgpoints(npoints,);
        break());
      );
  msgprec(abs(supres*npoints/res));
  res*=h;
  return(res);
};

longueur d'intégration

tmax_sh(D,M1,alpha,beta) = {
  my(a,b,C,X);
  a = 1/beta; b = alpha;
  C = M1/(beta*(D+1)-1);
  X = a/b*Lx(a/b*exp((D+log(C))/a));
  \\msgdebug(Strprintf("\told X=%1.2e\n",X),,DBG);
  \\X = Lab(D+log(M1),a,b,0);
  \\msgdebug(Strprintf("\tnew X=%1.2e\n",X),,DBG);
  return(asinh(X));
};
addhelp(tmax_sh,"tmax_sh(D,M1,alpha,beta)");

pas d'intégration

pas_h_sh(D,tau,M2,alpha=1,beta=1,A=0,gamma_=0) = {
  if(!D,D=default(realprecision)*log(10));
  pas_h_DE3(D,tau,M2,alpha,beta,/*lambda=*/1,A,gamma_);
};
addhelp(pas_h_sh,"pas_h_sh(D,tau,M2,alpha,beta,A,gamma_)");
integration_SE(f,M1,alpha,beta,tau,M2,lambda,A,gamma_,D) = {
  my(h,tmax,npoints);
  if(!D,D=default(realprecision)*log(10));
  h = pas_h_sh(D,tau,M2,alpha,beta,A,gamma_);
  tmax = tmax_sh(D,M1,alpha,beta);
  npoints = ceil(tmax/h);
  msgpoints(npoints,"SE");
  integration_sh(f,h,npoints)
};
addhelp(integration_SE,"integration_SE(f,M1,alpha,beta,tau,M2,lambda,A,gamma_,D)");

intégrale de f sur \R+iθ, par changement exponentiel,

mais en décalant la ligne d'intégration de θ. On doit donc avoir f=O( exp(-alpha x^beta) ) sur un secteur d'ouverture θ. I ~ sum_k [ h * sum f(sh(kh+I*θ))*ch(kh+I*θ) ] NB/ cette fois, ch(kh+I*θ) != ch(-kh+I*θ) f : fonction à intégrer θ : décalage de la droite d'intégration. h : pas d'intégration (cf. SE_int_decale pour valeur usuelle) npoints : nombre de termes sommés

integration_shteta(f,teta,h,npoints) = {
  my(res,ekh,i_ekh,sh,ch,cosd,sind,exph,i_exph,sup,tmp,supres);
  sup=0;
  exph=exp(h);i_exph=1/exph;
  ekh=i_ekh=1;
  cosd=cos(teta);
  sind=sin(teta);
  \\if(flag,
      \\res= sum(k=-n,n,f(sinh(k*h+I*d))*cosh(k*h+I*d))
      \\,
      res=f(0+I*sind)*cosd; /* =f(sinh(I*teta))*cosh(I*teta) */
      supres=0;
      for(k=1,npoints,
        ekh*=exph;
        i_ekh*=i_exph;
        sh=(ekh-i_ekh)/2;
        ch=(ekh+i_ekh)/2;
        /* sinh(kh+I*d)=sh(kh)cos(d)+I*ch(kh)sin(d) */
        /* cosh(kh+I*d)=ch(kh)cos(d)+I*sh(kh)sin(d) */
        shcos=sh*cosd; chsin=ch*sind;
        chcos=ch*cosd; shsin=sh*sind;
        tmp  = f( shcos+I*chsin)*(chcos+I*shsin);
        tmp += f(-shcos+I*chsin)*(chcos-I*shsin);
        if(default(debug)>=EXTDBG && abs(tmp)>sup,
          sup=abs(tmp)); \\ on évalue l'ordre des termes
        res+=tmp;
        supres=max(abs(res),supres);
        /* res = sum_{-k}^k f(sinh(kh+I*d))*cosh(kh+I*d) */
        );
      \\);
      /* res = h sum_{-n}^n f(sinh(kh+I*d))*cosh(kh+I*d) */
      msgdebug(Strprintf("[integration] sup = %1.3e",abs(sup)),,EXTDBG); \\ et on l'indique
        msgprec(abs(supres*npoints/res));
      res*=h;
      return(res);
};
addhelp(integration_shteta,"integration_shteta(f,teta,h,npoints)");
integration_shteta_until(f,teta,h,D,moreprec=2) = {
  my(res,ekh,i_ekh,sh,ch,cosd,sind,exph,i_exph,sup,tmp,supres,npoints);
  if(!D,D=precision(1.);moreprec=0);
  D = (D+moreprec)*log(10);
  sup=0;
  exph=exp(h);i_exph=1/exph;
  ekh=i_ekh=1;
  cosd=cos(teta);
  sind=sin(teta);
  \\if(flag,
      \\res= sum(k=-n,n,f(sinh(k*h+I*d))*cosh(k*h+I*d))
      \\,
      res=f(0+I*sind)*cosd; /* =f(sinh(I*teta))*cosh(I*teta) */
      supres=0;
      while(1,
        ekh*=exph;
        i_ekh*=i_exph;
        sh=(ekh-i_ekh)/2;
        ch=(ekh+i_ekh)/2;
        /* sinh(kh+I*d)=sh(kh)cos(d)+I*ch(kh)sin(d) */
        /* cosh(kh+I*d)=ch(kh)cos(d)+I*sh(kh)sin(d) */
        shcos=sh*cosd; chsin=ch*sind;
        chcos=ch*cosd; shsin=sh*sind;
        tmp = f(shcos+I*chsin)*(chcos+I*shsin);
        tmp+= f(-shcos+I*chsin)*(chcos-I*shsin);
        res+=tmp;
        supres=max(abs(res),supres);
        if(default(debug)>=EXTDBG && abs(tmp)>sup, sup=abs(tmp)); \\ on évalue l'ordre des termes
        if(tmp&&log(abs(tmp/res))+D+moreprec<0,
          npoints = round(log(ekh)/log(exph));
          msgpoints(npoints,"shteta");
          break());
        /* res = sum_{-k}^k f(sinh(kh+I*d))*cosh(kh+I*d) */
        );
  \\);
  /* res = h sum_{-n}^n f(sinh(kh+I*d))*cosh(kh+I*d) */
  msgdebug( Strprintf("sup effectivement rencontré : ",abs(sup)),,EXTDBG); \\ et on l'indique
    res*=h;
  return(res);
};

Changement de variable simplement exponentiel sur [0,oo[ en z=exp(z-α.exp(-βz))

On suppose que f=O(exp(-α.x^β)) sur Delta_tau

integration_exp(f,alpha,beta,h,npoints) = {
  my(ab,kh,res,ebh,ebh_inv,ebkh,ebkh_inv,phikh,supres);
  ebh = exp(-beta*h);
  ebh_inv = 1/ebh;
  ebkh = ebkh_inv = alpha;
  kh=0;
  res=(1+beta*ebkh)*exp(-alpha)*f(exp(-alpha));
  supres=0;
  for(k=1,npoints,
      ebkh *= ebh;
      ebkh_inv *= ebh_inv;
      kh+=h;
      phikh = exp(kh-ebkh);
      res+=(1+beta*ebkh)*phikh*f(phikh);
      phikh = exp(-kh-ebkh_inv);
      res+=(1+beta*ebkh_inv)*phikh*f(phikh);
      supres=max(abs(res),supres);
     );
  msgprec(abs(supres*npoints/res));
  res*=h;
  return(res);
};
addhelp(integration_exp,"integration_exp(f,alpha,beta,h,npoints)");
integration_exp_until(f,alpha,beta,h,D,moreprec=2) = {
  my(ab,kh,res,ebh,ebh_inv,ebkh,ebkh_inv,phikh,supres,npoints);
  if(!D,D=precision(1.);moreprec=0);
  D = (D+moreprec)*log(10);
  ebh = exp(-beta*h);
  ebh_inv = 1/ebh;
  ebkh = alpha*ebh;
  ebkh_inv = alpha*ebh_inv;
  kh=h;
  res=(1+beta*alpha)*exp(-alpha)*f(exp(-alpha));
  supres=0;
  \\wait = 0; \\ prevent premature stopping on oscillating functions
    while(1,
        phikh = exp(kh-ebkh);
        tmp=(1+beta*ebkh)*phikh*f(phikh);
        phikh = exp(-kh-ebkh_inv);
        tmp+=(1+beta*ebkh_inv)*phikh*f(phikh);
        res += tmp;
        supres=max(abs(res),supres);
        if(tmp&&log(abs(tmp/res))+D+moreprec<0,
          \\if(wait++<10,next());
          msgdebug(
            npoints = round(log(ebkh)/log(ebh));
            msgpoints(npoints,"exp");
            );
          break()
          \\,wait = 0
          );
        ebkh *= ebh;
        ebkh_inv *= ebh_inv;
        kh+=h;
        );
  msgprec(abs(supres*npoints/res));
  res*=h;
  return(res);
};
addhelp(integration_exp_until,"integration_exp_until(f,alpha,beta,h,D)");
pas_h_exp(D,tau,M2=1,alpha=1,beta=1,lambda=0,A=0,gamma_=0) = {
  if(!D,D=default(realprecision)*log(10));
  my(lambda=max(1,beta-1));
  return(pas_h_DE3(D,tau,M2*(1+alpha*beta),alpha,beta,lambda,A,gamma_))
};
tmax_exp(D,M1=1,alpha=1,beta=1,lambda=0) = {
  if(!D,D=default(realprecision)*log(10));
  my(tminus = 1/beta*log((D+log(M1))/alpha));
  my(X = Lab(D+log(M1),alpha,beta,lambda));
  my(tplus = log(X*exp(alpha/X^beta)));
  return(max(tplus,tminus));
};
integration_Rplus(f,M1,alpha,beta,tau,M2,lambda,A,gamma_,D) = {
  if(!D,D=default(realprecision)*log(10));
  my(h = pas_h_exp(D,tau,M2,alpha,beta,lambda,A,gamma_));
  my(npoints = ceil(tmax_exp(D,M1,alpha,beta,lambda)/h));
  msgpoints(npoints,"Rplus");
  integration_exp(f,alpha,beta,h,npoints)
};
addhelp(integration_Rplus,"integration_Rplus(f,M1,alpha,beta,tau,M2,lambda,A,gamma_,D)");

Changement de variable pour [-1,1] en z=tanh(Pi/2*sinh(z))

integration_thsh(f,h,npoints,lambda=Pi/2) = {
  my(k,exph,eh_inv,res,ekh,ekh_inv,sh,chi,thsh,supres);
  exph = exp(h);
  eh_inv = 1/exph;
  ekh = ekh_inv = 1;
  res = f(0);
  supres=0;
  for(k=1,npoints,
      ekh *= exph;
      ekh_inv *= eh_inv;
      sh = (ekh-ekh_inv)/2;
      ch2 = (ekh+ekh_inv); /* 2*cosh(kh) */
      esh = exp(lambda*sh);
      esh_inv = 1/esh;
      chsh2_i = 1/(esh+esh_inv); \\ 1/(2*cosh(lambda*sinh(kh)))
      shsh2 = esh-esh_inv; \\ 2*sinh(lambda*sinh(kh))
      thsh = shsh2*chsh2_i;
      res += 2*ch2*(f(thsh)+f(-thsh))*chsh2_i^2;
      supres=max(abs(res),supres);
     );
  msgprec(abs(supres*npoints/res));
  res *= h*lambda;
  return(res);
};
addhelp(integration_thsh,"integration_thsh(f,h,npoints,lambda=Pi/2)");
integration_thsh_until(f,h,D,lambda=Pi/2,moreprec=2) = {
  my(k,exph,eh_inv,res,ekh,ekh_inv,sh,chi,thsh,tmp,supres,npoints);
  if(!D,D=precision(1.);moreprec=0);
  D = (D+moreprec)*log(10);
  exph = exp(h);
  eh_inv = 1/exph;
  ekh = ekh_inv = 1;
  res = f(0);
  supres=0;
  while(1,
      ekh *= exph;
      ekh_inv *= eh_inv;
      sh = (ekh-ekh_inv)/2;
      ch2 = (ekh+ekh_inv); /* 2*cosh(kh) */
      esh = exp(lambda*sh);
      esh_inv = 1/esh;
      chsh2_i = 1/(esh+esh_inv); \\ 1/(2*cosh(lambda*sinh(kh)))
      shsh2 = esh-esh_inv; \\ 2*sinh(lambda*sinh(kh))
      thsh = shsh2*chsh2_i;
      tmp = 2*ch2*(f(thsh)+f(-thsh))*chsh2_i^2;
      res+=tmp;
      supres=max(abs(res),supres);
      if(log(abs(tmp/res))+D+moreprec<0,
        npoints = round(log(ekh)/log(exph));
        msgpoints(npoints,"thsh");
        break());
      );
  msgprec(abs(supres*npoints/res));
  res *= h*lambda;
  return(res);
};

Paramètre tau

best_tau_X(Xm,lambda=Pi/2) = {
  Xm=abs(Xm);
  if(Xm<=1,return(0));
  if(!lambda,lambda=Pi/2);
  /* for Xm */
  Xtau = atanh(1/Xm);
  Ytau = acos(1/Xm);
  imag(asinh((Xtau+I*Ytau)/lambda));
};
best_tau_Y(Ym,lambda=Pi/2) = {
  Ym=abs(Ym);
  if(!Ym,return(0));
  if(!lambda,lambda=Pi/2);
  /* for Ym */
  Xtau = asinh(1/sqrt(2*Ym));
  Ytau = atan(Ym);
  imag(asinh((Xtau+I*Ytau)/lambda));
};
/* tau for z such that z is not
   in Z_τ * mu */
tau_point(a,b,z,mu) = {
  my(x,y,z,tauX,tauY);
  z=(z-(a+b)/2)*2/(b-a);
  x = abs(real(z));
  y = abs(imag(z));
  tauX = best_tau_X((x+mu-1)/mu);
  tauY = best_tau_Y(y/mu);
  max(tauX,tauY)
};
/* find tau and sup M */
tau_vecpoints(a,b,V,mu=1.1) = {
  my(tau,tauM);
  tau=Pi/2;
  B = 0;
  for(i=1,#V,
      tauM = tau_M_point(a,b,V[i],mu);
      tau = min(tau,tauM);
     );
  return(tau);
  /* d >= (mu-1) * (b-a)/2 */
};
boundX(tau,lambda=Pi/2) = {
my(Xtau,Ytau,dec,ftau,eps,f2);
if(!lambda,lambda=Pi/2);
if(tau>=Pi/2,print("*** error boundX : tau>=Pi/2.");return(0));
if(tau<=0,print("*** error boundX : tau<=0.");return(0));
if(2*lambda*sin(tau)>=Pi,print("*** error boundX :
lambda*sin(tau)>=Pi/2.");return(0));
  \\ dichotomy /* method of the secant */
  X1 = 0; Y1 = lambda*sin(tau);
  f1 = -cos(Y1);
  X2 = sqrt(Pi^2/4-lambda^2*sin(tau)^2)/tan(tau); Y2 = Pi/2;
  f2=tanh(X2);
  dec=1; eps = 10^(3-default(realprecision));
  while(1, \\abs(Y2-Y1)>eps,
    Ytau = (Y1+Y2)/2; \\ Y1-f1*(Y2-Y1)/(f2-f1);
    Xtau = sqrt(Ytau^2-lambda^2*sin(tau)^2)/tan(tau);
    ftau = tanh(Xtau)-cos(Ytau);
    if(abs(ftau)<eps,break(1));
    if(ftau<0,Y1=Ytau;f1=ftau,Y2=Ytau;f2=ftau);
    \\print([Y1,f1,Ytau,ftau,Y2,f2]);
  );
  print([1/tanh(Xtau),tanh(Xtau)/cos(Ytau)^2]);
  return(tanh(Xtau)/cos(Ytau)^2);
  return(1/tanh(Xtau));
};
addhelp(boundX,"boundX(tau,lambda) : calcule Xmax");

paramètres d'intégration

tau_thsh(Xm,Ym,lambda=Pi/2) = {
  if(lambda==0,lambda=Pi/2);
  tauX=best_tau_X(Xm,lambda);
  tauY=best_tau_Y(Ym,lambda);
  \\print("tauX=",tauX," ; tauY=",tauY);
  min(tauX,tauY);
  /* by construction, lambda*sin(tau) < Pi/2 */
};
{addhelp(tau_thsh,"tau_thsh(Xmax,Ymax,lambda) : calcule tau
    tel que on reste dans le rectangle fixé")};
/* bound for int_\R |λch(x+iτ)|/|ch(λsh(x+iτ))^(2α)| */
phi_bound(tau,lambda=Pi/2,alpha=1/2) = {
  \\ fixed value of Ytau
    Ytau = sqrt(Pi/2*lambda*sin(tau));
  Xtau = sqrt(Ytau^2-lambda^2*sin(tau)^2)/tan(tau);
  return(2*(1+2*alpha)/(2*alpha*cos(Ytau)^(2*alpha))+2/(2*alpha*sinh(Xtau)^(2*alpha)));
};
pas_h_thsh_Xmax(D,tau,M=1,Xmax=3) = {
  if(!D,D=default(realprecision)*log(10));
  2*Pi*tau/(D+log(16*M*Xmax/cos(tau)+1))
};
addhelp(pas_h_thsh_Xmax,"pas_h_thsh_Xmax(D,tau,M,Xmax)");
pas_h_thsh(D,tau,M=1,alpha=1,lambda=Pi/2) = {
  if(!D,D=default(realprecision)*log(10));
  return(2*Pi*tau/(D+log(1+M*phi_bound(tau,lambda,alpha))));
};
addhelp(pas_h_thsh,"pas_h_thsh(D,tau,M=1,alpha=1,lambda=Pi/2)");
tmax_thsh(D,M,lambda=Pi/2) = {
  if(!D,D=default(realprecision)*log(10));
  asinh((D+log(8*M))/lambda)
};
addhelp(tmax_thsh,"tmax_thsh(D,M,lambda)");
/* integration sur [-1,1] */
integration_segment(f,a=-1,b=1,Xm,Ym,M,lambda=Pi/2,D) = {
my(h,tau,npoints,midpoint,geom_factor);
  if(!D,D=default(realprecision)*log(10));
  tau = tau_thsh(Xm,Ym,lambda);
  h = pas_h_thsh_Xmax(D,tau,M,Xm);
  npoints = ceil(tmax_thsh(D,M,lambda)/h);
  msgpoints(npoints,"thsh");
  midpoint=(a+b)/2; geom_factor=(b-a)/2;
  my(g(t)=f(midpoint+t*geom_factor));
  geom_factor*integration_thsh(g,h,npoints,lambda)
};
integration_segment_tau(f,a=-1,b=1,tau=Pi/4,M=1,lambda=Pi/2,D) = {
  my(h,npoints,midpoint,geom_factor);
  if(!D,D=default(realprecision)*log(10));
  my(alpha = 1);
  if(type(a)==t_VEC,alpha=min(alpha,a[2]);a=a[1]);
  if(type(b)==t_VEC,alpha=min(alpha,b[2]);b=b[1]);
  h = pas_h_thsh(D,tau,M,alpha,lambda);
  npoints = ceil(tmax_thsh(D,M,lambda)/h);
  msgpoints(npoints,"thsh");
  midpoint=(a+b)/2; geom_factor=(b-a)/2;
  my(g(t)=f(midpoint+t*geom_factor));
  geom_factor*integration_thsh(g,h,npoints,lambda)
};
addhelp(integration_segment_tau,"integration_segment_tau(f,a,b,tau,M,lambda=Pi/2)");

Changement de variable pour Rpol en z=sinh(lambda*sinh(z))

integration_shsh(f,h,npoints,lambda=Pi/2) = {
  my(k,exph,eh_inv,res,ekh,ekh_inv,lsh,ch,elsh,elsh_inv,shsh,l,supres);
  l = lambda/2;
  exph = exp(h);
  eh_inv = 1/exph;
  ekh = ekh_inv = 1;
  res = 4*f(0);
  supres=0;
  for(k=1,npoints,
      ekh *= exph;
      ekh_inv *= eh_inv;
      ch = (ekh+ekh_inv); \\ /2 at the end
      lsh = (ekh-ekh_inv)*l;
      elsh = exp(lsh);
      elsh_inv = 1/elsh;
      shsh = (elsh-elsh_inv)/2;
      res += ch*(elsh+elsh_inv)*(f(shsh)+f(-shsh));
      supres=max(abs(res),supres);
     );
  res *= h*l/2;
  msgprec(abs(supres*npoints/res));
  res
};
addhelp(integration_shsh,"integration_shsh(f,h,npoints,lambda=Pi/2)");
pas_h_shsh(D,tau,M2,u) = {
  if(!D,D=default(realprecision)*log(10));
  return(2*Pi*tau/(D+log(2*M2/(u*cos(tau)))));
};
addhelp(pas_h_shsh,"pas_h_shsh(D,tau,M2,u)");
tmax_shsh(D,M1,alpha) = {
  asinh(asinh(exp((D+log(2*M1/(alpha-1))/(alpha-1)))));
};
addhelp(tmax_shsh,"tmax_shsh(D,M1,alpha)");
integration_shsh_until(f,h,D,lambda=Pi/2,moreprec=2) = {
  my(k,exph,eh_inv,res,ekh,ekh_inv,lsh,ch,elsh,elsh_inv,shsh,l,tmp,supres,npoints);
  if(!D,D=precision(1.);moreprec=0);
  D = (D+moreprec)*log(10);
  l = lambda/2;
  exph = exp(h);
  eh_inv = 1/exph;
  ekh = ekh_inv = 1;
  res = 4*f(0);
  supres=0;
  while(1,
      ekh *= exph;
      ekh_inv *= eh_inv;
      ch = (ekh+ekh_inv); \\ /2 at the end
      lsh = (ekh-ekh_inv)*l;
      elsh = exp(lsh);
      elsh_inv = 1/elsh;
      shsh = (elsh-elsh_inv)/2;
      tmp = ch*(elsh+elsh_inv)*(f(shsh)+f(-shsh));
      res +=tmp;
      supres=max(abs(res),supres);
      if(log(abs(tmp/res))+D+moreprec<0,
        npoints = round(log(ekh)/log(exph));
        msgpoints(npoints,"thsh");
        break());
      );
  msgprec(abs(supres*npoints/res));
  res *= h*l/2;
  return(res);
};

Calculs de distances à φ-1(Δ_τ)

dist_sh(s,teta) = {
  my(xp,yp);
  xp = abs(real(s));  yp = abs(imag(s));
  ct = cos(teta); st = sin(teta); tt = st/ct;
  \\x = solve(x=xp, ct^2*(xp+yp*tt), x-xp-x*tt*(yp/sqrt(x^2+ct^2)-tt));
  x = solve(x=xp, sqrt(yp^2-st^2)/tt,x-xp-x*tt*(yp/sqrt(x^2+ct^2)-tt));
  y = tt*sqrt(ct^2+x^2);
  return(sqrt((xp-x)^2+(yp-y)^2));
};
addhelp(dist_sh,"dist_sh(s,teta)");
dist_thsh(p,tau,lambda=Pi/2) = {
  my(p,x0,x1,x);
  \\ take everything in the first quatran
    p = abs(real(p))+I*abs(imag(p));
  \\ restrict curve in the first quatran
    x0 = 0; x1 = acosh(Pi/(2*lambda*sin(tau)));\\ st. λch(x)sin(τ)=Pi/2
    \\ solve orthogonality (p-φ(x)).φ'(x) = 0,
    \\ but with arguments
  my(phi(x)=tanh(lambda*sinh(x+I*tau)));
  my(argphip(x)= arg(cosh(x+I*tau))-2*arg(cosh(lambda*sinh(x+I*tau))));
  my(argphi(x)=arg(p-phi(x)));
  \\print("arg[x0]=",argphi(x0)-argphip(x0)-Pi/2);
  \\print("arg[x1]=",argphi(x1)-argphip(x1)-Pi/2);
  x = solve(x=x0,x1, argphi(x)-argphip(x)-Pi/2);
  abs(p-phi(x));
};
addhelp(dist_thsh,"dist_thsh(p,tau,lambda=Pi/2)");

Détermination des bornes M2 par parcourt des frontières des zones Z_τ

/* TODO : option onlypath pour faire juste le chemin,
   ça permet de voir si on passe dans des pôles.
 */
plot_Ztau_sh(f,tau,tmax=10) = {
  /* run along the hyperbola */
  my(st = sin(tau));
  my(ct2 = (cos(tau))^2);
  my(
      phi_sup(x) = x+I*st*sqrt(x^2+ct2);
    );
  ploth(t=-tmax,tmax,
    my(z=phi_sup(t));
    [abs(f(z)),abs(f(conj(z)))]
  )[4] /* return max */
};
addhelp(plot_Ztau_sh,"plot_Ztau_sh(f,tau,tmax=10)");
plot_Ztau_exp(f,tau,alpha=1,beta=1,tmax=10) = {
  /* run along a truncated half cone */
  my(p2t = Pi/2-tau);
  my(tbt = tan(beta*tau));
  my(rm = -alpha*sin(beta*tau)/p2t*exp(-p2t/tbt));
  my(ym = tbt*(alpha-rm));
  my(
      phi_sup(t) = if(t<1,
        (-rm)+I*t*ym,
        (t*(alpha-rm)-alpha)+I*tbt*(t*(alpha-rm))
        )
    );
  ploth(t=0,tmax,
    my(z=phi_sup(t));
    [abs(f(z)),abs(f(conj(z)))]
  )[4]
};
addhelp(plot_Ztau_exp,"plot_Ztau_exp(f,alpha,beta,tau,tmax=10)");
plot_Ztau_shsh(f,tau,u=1,tmax=10) = {
  my(st = sin(tau));
  my(ct2 = (cos(tau))^2);
  my(
    phi_sup(x) = sinh(x+I*st*sqrt(x^2+ct2))
    );
  ploth(t=-tmax,tmax,
 my(z=phi_sup(t));my(crois=(1+abs(z)^(1+u)));
 [abs(f(z))*crois,abs(f(conj(z)))*crois]
)[4]
 };
addhelp(plot_Ztau_shsh,"plot_Ztau_shsh(f,tau,u=1,tmax=10)");
plot_Ztau_thsh(f,a,b,tau,tmax=10) = {
  my(st = sin(tau));
  my(ct2 = (cos(tau))^2);
  my(midpoint=(a+b)/2);my(geom_factor=(b-a)/2);
  my(
      phi_sup(x) = midpoint+geom_factor*tanh(x+I*st*sqrt(x^2+ct2))
    );
  ploth(t=0,tmax,my(z=phi_sup(t));
      [abs(f(z)),abs(f(conj(z))),abs(f(-z)),abs(f(-conj(z)))]
      )[4]
};

Corrections des résidus

/* Rho vector of poles, Res that of residues */
terme_poles(Rho,Res,h) = {
  my(i,res,rho);
  res=0;
  for(i=1,#Rho,
      rho = Rho[i];
      if(imag(rho)>0,
        \\ then take upper residue, positive loop
        res += Res[i]/(exp(-2*I*Pi*rho/h)-1),
        \\else lower residue, negative loop
        res -= Res[i]/(exp(2*I*Pi*rho/h)-1)
        );
     );
  return(2*I*Pi*res);
};
addhelp(terme_poles,"terme_poles(Rho,Res,h)");
poles_sh(Rho,Res,h) = {
  my(R);
  R=asinh(Rho);
  terme_poles(R,Res,h)
};
addhelp(poles_sh,"poles_SE(Rho,Res,h)");
/* for double sh one must take into account
   many inverse images of sinh :
   sinh(z)=sinh(I*Pi-z)=sin(2I*Pi+z)
 */
poles_shsh(Rho,Res,h,D) = {
  my(res,ash,k,sk,tmp,rho,eps);
  if(!D,D=precision(1.)*log(10));
  res=0;
  ash = asinh(Rho);
  k=0;
  while(1, \\ somme sur k
      if(k%2,eps=-1,eps=1);
      tmp=0;
      for(i=1,#ash, \\ sur les poles
        rho = asinh( I*k*Pi + eps*ash[i] );
        for(sk=0,1, \\ k et -k
          if(imag(rho)>0,
            \\ then take upper residue, positive loop
            tmp += Res[i]/(exp(-2*I*Pi*rho/h)-1),
            \\else lower residue, negative loop
            tmp -= Res[i]/(exp(2*I*Pi*rho/h)-1)
            );
          \\ if k=0 just one term, otherwise take -k
          if(!k,break(),rho = asinh( -I*k*Pi + eps*ash[i] ));
          );
        );
      res+=tmp;
      if(log(abs(tmp))+D<0,break());
      k++;
      );
  return(2*I*Pi*res);
};
addhelp(poles_shsh,"poles_shsh(Rho,Res,h,D)");
poles_thsh(a,b,Rho,Res,h,lambda=Pi/2) = {
  my(i,R,p);
  R=vector(#Rho);
  for(i=1,#Rho,
      p=(2*Rho[i]-a-b)/(b-a);
      R[i]=asinh(atanh(p)/lambda)
     );
  terme_poles(R,Res,h)
};
addhelp(poles_thsh,"poles_thsh(a,b,Rho,Res,h,lambda=Pi/2)");