Subversion Repositories wimsdev

!set slib_header_coxhyp=\

\\ Action du groupe de Moebius etendu\

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

\

hyp3(p,q,r,eps,depl)={\

  my(m=[1,p,q;p,1,r;q,r,1],nf0,nf1,nf2,f1,f2,nbs,snb,\

  ff1,s,t,u,g,gg,ggg,g1,x,y,z,res1,res2,ses1,r2,l,\

  ep=exp(I*(Pi/p)),\

  eq=exp(I*(Pi/q)),\

  er=exp(I*(Pi/r)),\

  a=0,\

  b=sqrt((cos(Pi/r)+cos(Pi/p+Pi/q))/(cos(Pi/r)+cos(Pi/p-Pi/q))),\

  c=ep*sqrt((cos(Pi/q)+cos(Pi/p+Pi/r))/(cos(Pi/q)+cos(Pi/p-Pi/r))),\

  mm=(b*c*conj(b-c)+c-b)/(conj(b)*c-conj(c)*b));\

  my(N=20); if(eps<1,\

    r2=(1-eps)^2;l=N*round(4*r2/(1-r2)/(1-1/p-1/q-1/r)),\

    l=eps; eps=4*N/(l*(1-1/p-1/q-1/r)+4*N));\

  my(cd0=vector(l), \\ pour chaque f, les numeros des 3 aretes du bord de la chambre f\

  ce0=vector(l), \\ et l'element g du groupe qui envoie la chambre 1 sur la chambre f\

  cd1=vector(2*l+1),\\ pour chaque a, les numeros des extremites de a\

  ce1=vector(2*l+1),\\ et le numero de la premiere chambre qu'elle borde\

  cd2=vector(l+2),\\ pour chaque s, le nombre d'aretes manquantes,\

  ce2=vector(l+2),\\ les deux aretes qui bordent le secteur restant\

  cf2=vector(l+2),\\ coordonnees du sommet\

  cg2=vector(l+2) \\ le sommet est dans la zone centrale\

  );\

\\ Les trois reflexions qui engendrent W\

  g=[[1,0,1],[ep,0,1],[I*mm,-I,1]];\

  if(depl,\

    a=mob(depl,a); b=mob(depl,b); c=mob(depl,c);\

    g=vector(3,i,mob_mul(depl,mob_mul(g[i],mob_inv(depl)))),\

    depl=[1,0,0]);\

\

  \\ Creation de la premiere chambre et de son squelette\

  cd0[nf0=1]=vector(3,s,s);\

  ce0[nf0]=[1,0,0]; \\ Identite\

  for (s = 1, 3,\

    cd1[nf1+=1] = vector(4, t, if(t>3,s,-1));\

    ce1[nf1]=1;\

    for (t=1, s-1,\

      cd2[nf2+=1] = 2*m[s,t]-2;\

      ce2[nf2] = [s,t];\

      snb=6-s-t;\

      cf2[nf2]=if(snb==1,c,if(snb==2,b,a));\

      cg2[nf2]=(sqrt(norm(cf2[nf2]))+eps)<1;\

      cd1[t][s] = cd1[s][t] = nf2));\

  while(f1<nf1 && nf0<l, f1+=1; v1=cd1[f1];\

   if((s=v1[4]) && (cg2[v1[s%3+1]] || cg2[v1[5-s-s%3]]), \\ arete a traiter\

\\ element du groupe associe au triangle precedent\

   gg=ce0[ce1[f1]];\

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

   ce0[nf0+=1]=ggg=mob_mul(gg,g[s]);\

  w0=vector(3); w0[s]=f1; new=vector(3,t,if(t!=s,cd2[v1[t]]));\

    for(t=1, 3, if(t!=s, f2=v1[t]; if(new[t],\

    w1=vector(4); w1[4]=t; w1[s]=f2; w1[t]=-1;\

    cd1[nf1+=1]=w1; ce1[nf1]=nf0; cd2[f2]=new[t]-1; ce2[f2]=w0[t]=nf1,\

    g1=w0[t]=ce2[f2]; cd1[g1][4]=0)));\

\\ Mise a jour (et creation si necessaire) des sommets non encore traites\

    for(t=1,3,if(new[t], g1=w0[t];\

      for(u=1,3,if((u!=t) && (u!=s),\

        if(new[u],\

      if(w0[u]<g1, cd2[nf2+=1]=2*m[t,u]-2; ce2[nf2]=g1;\

      snb=6-t-u;\

      cf2[nf2]=mob(ggg,if(snb==1,c,if(snb==2,b,a)));\

      cg2[nf2]=(sqrt(norm(cf2[nf2]))+eps)<1;\

      cd1[w0[u]][t]=cd1[g1][u]=nf2),\

      cd1[g1][u]=f2=cd1[w0[u]][t]; cd2[f2]-=1; ce2[f2]=g1)))));\

      cd0[nf0]=w0;\

      )\

  );\

  [matrix(nf2,2,i,j,if(j==1,real(cf2[i]),imag(cf2[i]))),\

   matrix(nf0,4,i,j,if(j<4,cd1[cd0[i][j%3+1]][5-j-j%3],ce0[i][3]))]\

}