Subversion Repositories wimsdev

!set slib_header_hyptiling=\

\\ -----------------------------------------------------------------------------------\

\\ Action du groupe de Moebius etendu sur le disque de Poincare\

\\ -----------------------------------------------------------------------------------\

\\ Homographie ou antihomographie selon que g[3] vaut 0 ou 1.\

mob(g,z)=if(g[3],z=conj(z));(g[1]*z+g[2])/(conj(g[2])*z+conj(g[1]));\

\\ Produit et inverse dans le groupe de Moebius etendu\

mobx(g)=[g[1],g[2];conj(g[2]),conj(g[1])];\

mob_mul(g1,g2)=my(m=mobx(g1)*mobx(if(g1[3],conj(g2),g2)));[m[1,1],m[1,2],g1[3]!=g2[3]];\

mob_inv(g)=if(g[3],[g[1],-conj(g[2]),1],[conj(g[1]),-g[2],0]);\

\\ Homographie (et anti) qui envoient (0,1) sur (a,b), ou a est un point du disque et b a l'horizon\

mob_ori(a,b,c)=my(e=sqrt(b),g1=(e-conj(e)*a)/(1-norm(a)));[g1,a*conj(g1),c];\

\\ Point d'intersection de la demi-geodesique de a vers b avec l'horizon\

hor(a,b)=my(z=mob([1,-a,0],b));mob([1,a,0],z/sqrt(norm(z)));\

\\ Homographie (et anti) qui envoient (a,b) sur (c,d), sous l'hypothèse d(a,b)=d(c,d)\

mob_abcd(a,b,c,d,e)=mob_mul(mob_ori(c,hor(c,d),e),mob_inv(mob_ori(a,hor(a,b))));\

\\ Homographie qui echange a et b\

mob_exg(a,b)={my(h=[1,a,0],h1=mob_inv(h),c=mob(h1,b));\

  mob_mul(h,mob_mul([I,-I*c,0],h1))};\

\\ Reflection par rapport à la geodesique qui joint a et b\

mob_ref(a,b)={\

  my(num=a*b*conj(a-b)+b-a,den=conj(a)*b-conj(b)*a,c,d,h);\

  d=if(norm(den)*1e8<norm(num),\

    if (norm(b)<norm(a),sqrt(a/conj(a)),sqrt(b/conj(b))),\

    c=num/den;(1+sqrt(1-norm(c)))/conj(c));\

  h=mob_ori(a,d);\

  mob_mul(h,mob_mul([1,0,1],mob_inv(h)))\

  };\

\\ "Distance" entre a et b, comprise entre 0 et 1. La vraie est log((1+d)/(1-d))\

dist(a,b)=sqrt(norm(mob([1,-a,0],b)));\

\\ Angle abc, compris entre -Pi et Pi\

angle(a,b,c)=my(g=[1,-b,0]); arg(mob(g,c)/mob(g,a));\

\

\\ ------------------------------------------------------------------------------------\

\\ Dessin d'arcs et de polygones\

\\ ------------------------------------------------------------------------------------\

\\ Trace (en tikz) la geodesique entre a et b\

tikz_arc(o,a,b)={\

  my(num=a*b*conj(a-b)+b-a,den=conj(a)*b-conj(b)*a,c,t,arga,argb);\

  if(norm(den)*1e8<=norm(num),\

    filewrite1(o,strprintf("--(%.4f,%.4f)",real(b),imag(b))),\

    c=num/den;\

    arga=arg(a-c)/Pi*180;\

    argb=arg(b-c)/Pi*180;\

    if (abs(arga-argb)>180,if(arga<argb,arga+=360,argb+=360));\

    filewrite1(o,strprintf("arc(%.4f:%.4f:%.4f)",arga,argb,sqrt(norm(c)-1))))\

};\

\\ Trace un polygone hyperbolique, rempli ou non, etiquete ou non\

tikz_poly(o,v,label,fill)={\

  my(n=#v,z);\

  if (label > 0,\

    z=vecsum(v)/n;filewrite1(o,strprintf("\\%s(%.4f,%.4f)node{\\tiny$%d$}(%.4f,%.4f)",\

      if(fill,"fill","draw"),real(z),imag(z),label,real(v[1]),imag(v[1]))),\

    filewrite1(o,strprintf("\\%s(%.4f,%.4f)",\

      if(fill,"fill","draw"),real(v[1]),imag(v[1]))));\

  for(i=1,n,tikz_arc(o,v[i],v[i%n+1]));\

  filewrite(o,";")\

};\

\

\\ ------------------------------------------------------------------------------\

\\ Pavage obtenu a partir d'un polygone convexe pavant\

\\ ------------------------------------------------------------------------------\

\\ Entree:  v_i, points du disque de Poincare et d_i>2 des entiers\

\\ On suppose que les v_i forment le bord oriente d'un polygone convexe (pave) P0\

\\ dont l'angle intérieur au point v_i est  2Pi/d_i\

\\ Si deux cotes consecutifs sont inegaux, le d_i entre eux est suppose pair\

\\ Renvoie sous forme de fichier off\

\\ tous les paves dont au moins un sommet est dans le disque euclidien D(0,1-eps)\

\\ Si eps>=1, on s'en sert comme limite sur le nombre de paves.\

\\ Si z0 est présent, c'est un point du pavé de base dont les transformés seront calculés\

\\ Comme base du pavage dual\

\

hyp_pav(v,d,eps,z0)={\

  my(n=#v,a,b,s,s1,s2,t,t1,t2,g,gg,cd0,cd1,cd2,r2,f1,v1,w0,z,limit);\

  if(eps>=1,limit=eps;eps=0,limit=1000);\

  r2=(1-eps)^2;\

  if(!v,return(0));\

\

\\ Les reflexions qui engendrent W\

  g=vector(n,i,mob_ref(v[i],v[i%n+1]));\

\

\\ Creation du premier polygone et de son squelette\

  \\ Pour chaque sommet, [arite restante, derniere arete cree, affixe, dans le disque?]\

  cd2=List(vector(n,s,[d[s]-2,if(s==1,n,s),v[s],norm(v[s])<r2]));\

  \\ Pour chaque arete, [origine, extremite, type, premier polygone borde, active]\

  cd1=List(vector(n,a,[a,a%n+1,a,1,cd2[a][4]||cd2[a%n+1][4]]));\

  \\ Pour chaque polygone, ses aretes et l'élément du groupe\

  cd0=List([vector(n+1,a,if(a>n,[1,0,0],a))]);\

\

  while(#cd0<limit&&f1<#cd1, f1+=1; v1=cd1[f1]; if(v1[5],\

    s=v1[3]; \\ type de l'arete (ambigu)\

\\  Nouveau pave et ses aretes, manquantes ou pas\

    w0=vector(n+1);\

    w0[n+1]=gg=mob_mul(cd0[v1[4]][n+1],g[s]);\

    w0[s]=a=f1; cd1[a][5]=0;\

    s1=cd1[a][1]; s2=cd1[a][2];\

    b=cd0[v1[4]][s%n+1];\

    if(s2!=cd1[b][1]&&s2!=cd1[b][2],t=s1;s1=s2;s2=t); \\ arete mal orientee\

    t1=(s-2)%n+1;\

    while(!cd2[s1][1]&&!w0[t1], \\ aretes precedent f1 qui sont deja la\

      w0[t1]=a=cd2[s1][2]; cd1[a][5]=0;\

      s1=cd1[a][1]+cd1[a][2]-s1; t1=(t1-2)%n+1);\

    t2=s%n+1;\

    while(!cd2[s2][1]&&!w0[t2], \\ aretes suivant f1 qui sont deja la\

      w0[t2]=a=cd2[s2][2]; cd1[a][5]=0;\

      s2=cd1[a][1]+cd1[a][2]-s2; t2=t2%n+1);\

    while(t1!=t2, \\ Il y a au moins deux aretes a creer\

      if(s1<s2,\

        z=mob(gg,v[t1]);\

        listput(cd2,[d[t1]-1,#cd1+1,z,norm(z)<r2]);\

        listput(cd1,[s1,#cd2,t1,#cd0+1,cd2[s1][4]||cd2[#cd2][4]]); w0[t1]=#cd1;\

        cd2[s1][1]-=1; cd2[s1][2]=#cd1; s1=#cd2; t1=(t1-2)%n+1,\

        z=mob(gg,v[t2%n+1]);\

        listput(cd2,[d[t2%n+1]-1,#cd1+1,z,norm(z)<r2]);\

        listput(cd1,[s2,#cd2,t2,#cd0+1,cd2[s2][4]||cd2[#cd2][4]]); w0[t2]=#cd1;\

        cd2[s2][1]-=1; cd2[s2][2]=#cd1; s2=#cd2; t2=t2%n+1));\

      listput(cd1,[s1,s2,t1,#cd0+1,cd2[s1][4]||cd2[s2][4]]); w0[t1]=#cd1;\

      cd2[s1][1]-=1; cd2[s1][2]=#cd1; cd2[s2][1]-=1; cd2[s2][2]=#cd1;\

      listput(cd0,w0)));\

  for(f=1,#cd0,\

    w=cd0[f];a=cd1[w[1]];b=cd1[w[n]];\

    s=a[1];if(s!=b[1]&&s!=b[2],s=a[2]);\

    for(i=1,n,a=cd1[w[i]];w[i]=s;s=a[1]+a[2]-s);\

    cd0[f]=w);\

 \\ Si les d_i ne sont pas tous pairs, la derniere colonne n'a pas de sens\

  res=vector(if(z0,4,2));\

  res[1]=matrix(#cd2,2,i,j,if(j==1,real(cd2[i][3]),imag(cd2[i][3])));\

  res[2]=matrix(#cd0,n+2,i,j,if(j==1,n,if(j<=n+1,cd0[i][j-1],cd0[i][n+1][3])));\

  res\

};\

\

\\ Dessine le pavage en tikz. r est un vecteur a deux composantes "off"\

tikz_off(fname,r,rayon,fill,labels)={\

  my(pts=r[1],pvs=r[2],n=#pvs[1,]-2,o=fileopen(fname,"w"));\

  filewrite(o,strprintf("\\begin{tikzpicture}[scale=5]\n\\draw(0,0)circle(1);\n"));\

  if(rayon,filewrite(o,strprintf("\\draw(0,0)circle(%.4f);\n",rayon)));\

  for(i=1,#pvs[,1],\

    tikz_poly(o,vector(n,j,pts[pvs[i,j+1],1]+I*pts[pvs[i,j+1],2]),if(labels,i),fill&&pvs[i,n+2]));\

  if (#r>2, pts=r[3]; pvs=r[4]; for(i=1,#pvs[,1],\

    tikz_poly(o,vector(pts[i,1],j,pts[pvs[i,j+1],1]+I*pts[pvs[i,j+1],2]),if(labels,i))));\

filewrite(o,strprintf("\\end{tikzpicture}"));\

  fileclose(o);\

};\

\

\\ -----------------------------------------------------------------------------------\

\\ Creation de polygones convexes d'angles et longueurs donnees\

\\ -----------------------------------------------------------------------------------\

\\ S'il existe un triangle de côtés a,b,c et angles A,B,C, les fonctions suivantes\

\\ renvoient les paramètres manquants. Sinon, elles renvoient 0\

abc(A,B,C)={\

  if (A+B+C >= Pi, return(0));\

  my(cha=(cos(A)+cos(B)*cos(C))/(sin(B)*sin(C)),\

    chb=(cos(B)+cos(A)*cos(C))/(sin(A)*sin(C)),\

    chc=(cos(C)+cos(B)*cos(A))/(sin(B)*sin(A)));\

  [sqrt((cha-1)/(cha+1)),sqrt((chb-1)/(chb+1)),sqrt((chc-1)/(chc+1))]\

};\

ABC(a,b,c)={\

  if(c>=a+b-a*b || b>=a+c-a*c || a>=b+c-b*c, return (0));\

  my(cha=(1+a^2)/(1-a^2),chb=(1+b^2)/(1-b^2),chc=(1+c^2)/(1-c^2),\

     sha=2*a/(1-a^2),shb=2*b/(1-b^2),shc=2*c/(1-c^2));\

  [acos((chb*chc-cha)/(shb*shc)),\

   acos((cha*chc-chb)/(sha*shb)),\

   acos((chb*cha-chc)/(shb*sha))]\

};\

AbC(a,B,c)={\

  my(cha=(1+a^2)/(1-a^2),chc=(1+c^2)/(1-c^2),sha=2*a/(1-a^2),shc=2*c/(1-c^2),\

     chb=cha*chc-sha*shc*cos(B),shb=sqrt(chb^2-1));\

  [acos((chb*chc-cha)/(shb*shc)),sqrt((chb-1)/(chb+1)),acos((chb*cha-chc)/(shb*sha))]\

};\

aBc(A,b,C)={\

  my(chb=(1+b^2)/(1-b^2),cB=sin(A)*sin(C)*chb-cos(A)*cos(C),sB,cha,chc);\

  if (abs(cB)>=1,return(0),sB=sqrt(1-cB^2));\

  cha=(cos(A)+cB*cos(C))/(sB*sin(C));\

  chc=(cos(C)+cB*cos(A))/(sB*sin(A));\

  if(cha>1 && chc>1,[sqrt((cha-1)/(cha+1)),acos(cB),sqrt((chc-1)/(chc+1))],0)\

};\

\

\\ S'il existe, renvoie le quadrilatère convexe [0,l,z3,z4]\

\\ d'angles alpha,beta,gamma,delta et longueur(AB)=l, sinon renvoie 0\

quad(al,be,ga,de,l)={\

  my(eps=1e-10,g=[1,l,0],t=aBc(al,l,be),eA=exp(I*al),eB=-exp(-I*be),\

      u,v,lC,lD,minC,maxC=1.);\

  if(t, u=abc(Pi-ga,t[2],Pi-de);\

    if(!u||u[1]>=t[3]||u[3]>=t[1], return(0));\

    lD=(t[3]-u[1])/(1-u[1]); lC=(t[1]-u[3])/(1-u[3]),\

    while(maxC-minC>eps, lC=(maxC+minC)/2; u=AbC(l,be,lC);\

      if(u[1]>=ga, minC=lC, v=aBc(al-u[3],u[2],ga-u[1]);\

      if(!v,maxC=lC, lD=v[3]; if(v[2]<de, maxC=lC, minC=lC)))));\

  [0,l,mob(g,lC*eB),lD*eA]\

};\

\

\\ La tortue hyperbolique, turn left beta, then forward l\

lft_fwd(g,beta,l)=my(e=exp(I*beta/2));mob_mul(g,[e,e*l,0]);\

\

\\ Entree: n angles et n-3 longueurs\

\\ Sortie: [z1=0,z2=l1,z3,....,zn] n points du disque de Poincare qui forment le bord\

\\ oriente d'un polygone convexe tel que... s'il existe. Sinon, 0.\

polygone(a,l)={\

  my(n=#a,res=vector(n),theta,phi,v,g,z,eth,w);\

  if(n==3,v=abc(a[1],a[2],a[3]); return(if(v,[0,v[3],v[2]*exp(I*a[1])])));\

  g=[1,l[1],0]; res[2]=z=l[1];\

  for(i=2,n-3, g=lft_fwd(g,Pi-a[i],l[i]); z=mob(g,0);\

    if(arg(z)<=theta,return(0),theta=arg(z));\

    res[i+1]=z);\

  phi=Pi-angle(mob(g,1),z,0);\

  if(theta>=a[1]||phi>=a[n-2], return(0));\

  v=quad(a[1]-theta,a[n-2]-phi,a[n-1],a[n],sqrt(norm(z)));\

  if(!v,return(0));\

  eth=exp(I*theta); res[n-1]=eth*v[3]; res[n]=eth*v[4];\

  return(res)\

};\

\

\\ Entree: n angles entre 0 et Pi, dont la somme est inferieure a (n-2)*Pi\

\\ Sortie: L'unique polygone convexe dont tous les cotes sont tangents a un meme cercle\

\\   de centre 0,  d'angles interieurs a_i et dont le premier sommet est un reel positif.\

\

tangentiel(a)={\

  my(l,cs=apply(x->cos(x/2),a),\

    k=solve(t=0,1,vecsum(apply(x->asin(t*x),cs))-Pi),\

    R=sqrt(2/(1+k)-1),\

    g=lft_fwd(lft_fwd([1,0,0],-asin(k*cs[1]),R),Pi/2,0),\

\\    g=lft_fwd(lft_fwd([1,0,0],0,R),Pi/2,0),\

    res=vector(#a));\

  for(i=1,#a,\

    l=abc(Pi/2,a[i]/2,asin(k*cs[i]))[3];\

    g=lft_fwd(g,0,l);\

    res[i]=mob(g,0);\

    g=lft_fwd(g,Pi-a[i],l));\

    res\

};\

\

\\ ----------------------------------------------------------------------------\

\\ Quelques paves simples\

\\ ----------------------------------------------------------------------------\

\

\\ Construit un cerf-volant [0,z2,z3,z4] d'angles [a1,a2,a3,a4==a2]\

build_kite(a1,a2,a3)={my(tr=polygone([a1/2,a2,a3/2]));\

  if(tr,[tr[1],tr[2],tr[3],mob(mob_ref(tr[1],tr[3]),tr[2])])};\

\

\\  d2 doit etre pair et 1/d1+2/d2+1/d3 < 1\

kite(d1,d2,d3,eps,depl)={\

  my(p=build_kite(2*Pi/d1,2*Pi/d2,2*Pi/d3));\

  if(!p,return(0));\

  if(depl,p=vector(4,i,mob(depl,p[i])));\

  hyp_pav(p,[d1,d2,d3,d2],eps)};\

\

\\ Polygone a n cotes tous egaux, avec angles tous egaux a 2Pi/d\

\\ On suppose donc (d-2)(n-2)>4\

regular(n,d,eps,depl=[1,0,0])={\

  my(t=polygone([2*Pi/n,Pi/d,Pi/d]));\

  if(t,hyp_pav(vector(n,k,mob(depl,exp(2*I*k*Pi/n)*t[2])),vector(n,k,d),eps))};\

\

\\ 1/d1 + 1/d2+ 1/d3 < 1/2\

\\ Si un des trois est impair, les autres sont supposes egaux\

triangle(d1,d2,d3,eps,depl=[1,0,0])={\

  my(t=polygone([2*Pi/d1,2*Pi/d2,2*Pi/d3]));\

  if(t&&depl,t=vector(3,i,mob(depl,t[i])));\

  if(t,hyp_pav(t,[d1,d2,d3],eps))};\

\

\\ genre de parallelogramme: deux angles distincts\

\\ depend d'un parametre mu continu compris entre 0 et 1\

\\ On doit avoir 1/d1 + 1/d2 < 1/2\

\\ Si d1 ou d2 est impair, on doit aussi avoir mu=1/2.\

parallelogram(d1,d2,mu,eps,depl)={\

  my(t=polygone([2*Pi/d1,2*mu*Pi/d2,2*Pi*(1-mu)/d2]));\

  if(!t,return(0));\

  if(depl,t=vector(3,i,mob(depl,t[i])));\

  hyp_pav([t[1],t[2],mob(mob_exg(t[2],t[3]),t[1]),t[3]],[d1,d2,d1,d2],eps)};\

\

\\ Etant donnes le centre et le rayon (<1) d'un disque hyperbolique,\

\\ renvoie le centre et le rayon du meme disque considere comme euclidien.\

hyp_cercle(c, r)=my(R2=norm(c),den=1-r^2*R2);[(1-r^2)*c/den,(1-R2)*r/den];\

\

\\ On se donne n entiers d[i] >= 3, avec la convention d[] périodique modulo n\

\\ On suppose que si d[i] est impair, alors on a d[i-1]=d[i+1]\

\\ Calcule le pavage associé au polygone tangentiel d'angles interieurs 2*Pi/d[i]\

\\ On peut afficher le pavage dual et les cercles inscrits /* TO DO */\

catalan(d,eps,depl)={\

  my(h=tangentiel(apply(x->2*Pi/x,d)));\

  if(depl,h=apply(x->mob(depl,x),h));\

  hyp_pav(h,d,eps);\

}\

\

\\ Etant donne le "fichier off" d'un pavage issu de catalan\

\\ Calcule un "fichier off dual" dont le premier sommet est 0\

cat_dual(off)={\

  my(cts=vector(#off[2]),fls=vector(#off[1]));\

  for(i=1,#off[2],0);\

/* TO DO */\

\

  [cts,fls]\

};\

\

/*\

tikz_off("8_3_5.tex",kite(8,3,5,0.01))\

regular(7,3,0.1);\

parallelogram(4,6,3,0.1);\

*/\