WebSVN – wimsdev – Diff – /trunk/wims/public_html/scripts/gp/hyptiling.gp

 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\
+\\ Homographie (et anti) qui envoie 0 sur a et 1 sur b, avec a \in D 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));\
+lft_fwd(g,beta,l)=my(e=exp(I*beta/2));mob_mul(g,[e,e*l,0]);\
 \
-\\ ------------------------------------------------------------------------------------\
-\\ Dessin d'arcs et de polygones\
-\\ ------------------------------------------------------------------------------------\
-\\ Trace (en tikz) la geodesique entre a et b\
-tikz_arc(o,a,b)={\
+wims_export(pavage)={\
-  my(num=a*b*conj(a-b)+b-a,den=conj(a)*b-conj(b)*a,c,t,arga,argb);\
+  my(res1,res2,res3,res4,[sommets,paves,aretes,daretes]=pavage,n=#paves[1]-1);\
-  if(norm(den)*1e8<=norm(num),\
-    filewrite1(o,strprintf("--(%.4f,%.4f)",real(b),imag(b))),\
+  res1=matrix(#sommets,2,i,j,if(j==1,real(sommets[i]),imag(sommets[i])));\
-    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))))\
+  res2=matrix(#paves,n+2,i,j,if(j==1,n,if(j==n+2,paves[i][n+1][3],paves[i][j-1])));\
-};\
-\\ Trace un polygone hyperbolique, rempli ou non, etiquete ou non\
-tikz_poly(o,v,label,fill)={\
-  my(n=#v,z);\
+  [res1,res2]\
-  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)={\
+hyp_pav(v,d,eps)={\
-  my(n=#v,a,b,s,s1,s2,t,t1,t2,g,gg,cd0,cd1,cd2,r2,f1,v1,w0,z,limit);\
+  my(n=#v,a,b,s,s1,s2,t,t1,t2,g,gg,cd0,cd1,cdd1,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\
+\\ Creation du premier pave 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 arete duale [pave origine, pave extremite, arete croisee]\
+  cdd1=List();\
-  \\ Pour chaque polygone, ses aretes et l'élément du groupe\
+  \\ Pour chaque pave, 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;\
+    w0[s]=a=f1; cd1[a][5]=0; listput(cdd1,[v1[4],#cd0+1,a]);\
     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;\
+      w0[t1]=a=cd2[s1][2]; cd1[a][5]=0; listput(cdd1,[cd1[a][4],#cd0+1,a]);\
       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;\
+      w0[t2]=a=cd2[s2][2]; cd1[a][5]=0; listput(cdd1,[cd1[a][4],#cd0+1,a]);\
       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]);\
   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])));\
+  [vector(#cd2,i,cd2[i][3]),Vec(cd0),Vec(cd1),Vec(cdd1)]\
-  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)={\
+tangentiel(a,centre)={\
   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),\
+    R=sqrt(2/(1+k)-1), g=[1,0,0], h, res=vector(#a));\
-    g=lft_fwd(lft_fwd([1,0,0],-asin(k*cs[1]),R),Pi/2,0),\
+  if(centre,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\
+  [res,R]\
 };\
+\
-\\ ----------------------------------------------------------------------------\
-\\ 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)={\
+catalan(d,limite,centre)={\
-  my(h=tangentiel(apply(x->2*Pi/x,d)));\
+  my([h,R]=tangentiel(apply(x->2*Pi/x,d),centre));\
-  if(depl,h=apply(x->mob(depl,x),h));\
-  hyp_pav(h,d,eps);\
+  [hyp_pav(h,d,limite),R,if(centre,0,R*I)]\
-}\
-\
-\\ 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);\
-*/\

Subversion Repositories wimsdev

(root)/trunk/wims/public_html/scripts/gp/hyptiling.gp – Rev 18506 → 18513