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 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];\

\\ Reflexion par rapport a 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)))\

  };\

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

\

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

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

  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 pave et de son squelette\

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

  cd2=List(vector(n,s,[d[s]-2,if(s==1,n,s),v[s],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 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; 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; 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; 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,t1,norm(z)<r2]);\

        listput(cd1,[s1,#cd2,t1,#cd0+1,cd2[s1][5]||cd2[#cd2][5]]); 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,t2%n+1,norm(z)<r2]);\

        listput(cd1,[s2,#cd2,t2,#cd0+1,cd2[s2][5]||cd2[#cd2][5]]); 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][5]||cd2[s2][5]]); w0[t1]=#cd1;\

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

      listput(cd0,w0)));\

  centres=vector(#cd0);\

  for(f=1,#cd0,\

    w=cd0[f]; g=w[#w]; centres[f]=mob(g,z0); w[#w]=g[3];\

    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);\

  [vector(#cd2,i,cd2[i][3..4]),Vec(cd0),Vec(cd1),centres,Vec(cdd1)]\

};\

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))]\

};\

tangentiel(a,auCentre)={\

  my(l,cs=apply(x->cos(x/2),a),rays=vector(1+#a),\

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

    R=sqrt(2/(1+k)-1), g=[1,0,0], h, res=vector(#a));\

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

  for(i=1,#a,\

    rays[i]=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));\

  rays[#rays]=R;\

  [res,rays]\

};\

\

catalan(d,limite,auCentre)={\

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

  [hyp_pav(h,d,if(auCentre,0,rays[#rays]*I),limite),rays];\

};\

\

wims_catalan(d,limite,auCentre,option)={\

  my([pavage,rays]=catalan(d,limite,auCentre));\

  my([sommets,paves,aretes,centres,daretes]=pavage,R,z,p1,p2);\

  my(res0,res1,res2,res3,res4,n=#paves[1]-1,temp);\

  if(option==1 || option==3 ,\

    res1=matrix(#sommets,2,i,j,if(j==1,real(sommets[i][1]),imag(sommets[i][1])));\

    res0=matrix(#paves,n+2,i,j,if(j==1,n,paves[i][j-1]));\

    res2=List(); for(i=1,#paves,listput(res2,Vec(res0[i,])));\

  );\

  if(option==2 || option==3,\

    res3=matrix(#paves,2,i,j,if(j==1,real(centres[i]),imag(centres[i])));\

    temp=vector(#paves,i,List());\

    foreach(daretes,x,listput(temp[x[1]],x[2]));\

    res4=List();\

    for(i=1,#paves,if(#temp[i],listinsert(temp[i],i,1);listput(res4,Vec(temp[i])))));\

  [res1,Vec(res2),res3,Vec(res4)];\

};\