Projective-Transformation.sagews

A=matrix(3,3,[0,1,1,0,0,1,-1,1,1]);
B=matrix(3,3,[0,2/3,2/3,0,0,2/3,-1/3,2/3,2/3])
show(A)
show(B)
C=matrix(B*(A^(-1)))
show(C)

$\displaystyle \left(\begin{array}{rrr} 0 & 1 & 1 \\ 0 & 0 & 1 \\ -1 & 1 & 1 \end{array}\right)$

$\displaystyle \left(\begin{array}{rrr} 0 & \frac{2}{3} & \frac{2}{3} \\ 0 & 0 & \frac{2}{3} \\ -\frac{1}{3} & \frac{2}{3} & \frac{2}{3} \end{array}\right)$

$\displaystyle \left(\begin{array}{rrr} \frac{2}{3} & 0 & 0 \\ 0 & \frac{2}{3} & 0 \\ \frac{1}{3} & 0 & \frac{1}{3} \end{array}\right)$

#Entonces la transformación proyectiva que manda los puntos (0,0,1), (1,0,1), (1,1,1); (2,1,1)
#a los puntos (0,0,1), (1,0,1),(1,1,0); (4/3,2/3,1) está dada por la matriz C
#La transformación del plano en el plano asociada es la transformación
#f(x,y)=(  ((2/3)x+0y+0)/(1x-(2/3)y-1/3),  (0x+(2/3)y +0)/(1x-(2/3)y-1/3))=(2x/(3x-2y-1),2y/(3x-2y+1))

var('x y')
def myf(x,y):
    return (2*x/(x+1),2*y/(x+1))
show(myf)

(x, y) <function myf at 0x7f1197d1fde8>

for i,j in ((0,0),(1,0),(1,1),(2,1)):
    show(myf(i,j))

($\displaystyle 0$, $\displaystyle 0$)

($\displaystyle 1$, $\displaystyle 0$)

($\displaystyle 1$, $\displaystyle 1$)

($\displaystyle \frac{4}{3}$, $\displaystyle \frac{2}{3}$)

maxiter=30
G=Graphics()
for i in range(maxiter):
    for j in range(i+1):
        if (i,j)==(0,0):
            show("(0,0) should be there")
        G+=point(myf(i,j),size=40/(i+1))
H=Graphics()
for i in range(int(maxiter)):
    H+=line([(2,2),(myf(i,0))],color='red',thickness=.1)
for i in range(int(maxiter)):
    H+=line([(2,0),(myf(i,i))],color='red',thickness=.1)
(H+G).show(aspect_ratio=1)

(0,0) should be there

p=(1,2)
v=vector([10,10])
myf(*v)

(20/11, 20/11)

#As an application, we draw the Farey Tree using this data just to see how it looks
maxlevel=5
show("maxlevel=",maxlevel)
FN=[]
FT=point((0,0))
for i in range(maxlevel):
    FN.append([])
FN.insert(0,[vector((1,0)),vector((1,1))])
for level in range(1,maxlevel+1):
    for i in range(len(FN[level-1])-1):
        (FN[level]).append(FN[level-1][i])
        (FN[level]).append(FN[level-1][i]+FN[level-1][i+1])
    (FN[level]).append(FN[level-1][i+1])
for i in range(maxlevel+1):
    if i<=5:
        show("level=",i,":",FN[i])
def offsprings(point,level):
    i=FN[level+1].index(point)
    return [(FN[level+1][i-1],level+1),(FN[level+1][i+1],level+1)]
FT=point((0,0))
#A point is a vector in the plane
#A pl is a pair of a point and a level in which that point lies.
def drawfareytree(pl):
    """
    Draws the Farey tree in the plane starting at pl;
    pl is pair of a vector in the plane, and the level
    """
    global FT
    p=pl[0];#The point
    l=pl[1];#The level
    #show(p)
    #show("current level=",l,"; maxlevel=",maxlevel)
    if l>=maxlevel:
        #show("level>=maxlevel")
        return
    else:
        #show("se ejecutó esta primero")
        leftchild,rightchild=offsprings(*pl)
        #show("Se ejecutó esta linea")
        FT+= point(myf(*p),color=hue(l/maxlevel),size=.3)
        FT+= point(myf(*leftchild[0]),color=hue((l-1)/(maxlevel-1)),size=.3)+\
            line([myf(*p),myf(*leftchild[0])],color=hue((l-1)/(maxlevel-1)),thickness=.3) #+text(str(leftchild[0]),myf(*leftchild[0]),fontsize=5)
        FT+= point(myf(*rightchild[0]),color=hue((l-1)/(maxlevel-1)),size=.3)+\
            line([myf(*p),myf(*rightchild[0])],color=hue((l-1)/(maxlevel-1)),thickness=.3)#+         text(str(rightchild[0]),myf(*rightchild[0]),fontsize=5)
        drawfareytree(leftchild)
        drawfareytree(rightchild)

maxlevel= $\displaystyle 5$

level= $\displaystyle 0$ : [$\displaystyle \left(1,\,0\right)$, $\displaystyle \left(1,\,1\right)$]

level= $\displaystyle 1$ : [$\displaystyle \left(1,\,0\right)$, $\displaystyle \left(2,\,1\right)$, $\displaystyle \left(1,\,1\right)$]

level= $\displaystyle 2$ : [$\displaystyle \left(1,\,0\right)$, $\displaystyle \left(3,\,1\right)$, $\displaystyle \left(2,\,1\right)$, $\displaystyle \left(3,\,2\right)$, $\displaystyle \left(1,\,1\right)$]

level= $\displaystyle 3$ : [$\displaystyle \left(1,\,0\right)$, $\displaystyle \left(4,\,1\right)$, $\displaystyle \left(3,\,1\right)$, $\displaystyle \left(5,\,2\right)$, $\displaystyle \left(2,\,1\right)$, $\displaystyle \left(5,\,3\right)$, $\displaystyle \left(3,\,2\right)$, $\displaystyle \left(4,\,3\right)$, $\displaystyle \left(1,\,1\right)$]

level= $\displaystyle 4$ : [$\displaystyle \left(1,\,0\right)$, $\displaystyle \left(5,\,1\right)$, $\displaystyle \left(4,\,1\right)$, $\displaystyle \left(7,\,2\right)$, $\displaystyle \left(3,\,1\right)$, $\displaystyle \left(8,\,3\right)$, $\displaystyle \left(5,\,2\right)$, $\displaystyle \left(7,\,3\right)$, $\displaystyle \left(2,\,1\right)$, $\displaystyle \left(7,\,4\right)$, $\displaystyle \left(5,\,3\right)$, $\displaystyle \left(8,\,5\right)$, $\displaystyle \left(3,\,2\right)$, $\displaystyle \left(7,\,5\right)$, $\displaystyle \left(4,\,3\right)$, $\displaystyle \left(5,\,4\right)$, $\displaystyle \left(1,\,1\right)$]

level= $\displaystyle 5$ : [$\displaystyle \left(1,\,0\right)$, $\displaystyle \left(6,\,1\right)$, $\displaystyle \left(5,\,1\right)$, $\displaystyle \left(9,\,2\right)$, $\displaystyle \left(4,\,1\right)$, $\displaystyle \left(11,\,3\right)$, $\displaystyle \left(7,\,2\right)$, $\displaystyle \left(10,\,3\right)$, $\displaystyle \left(3,\,1\right)$, $\displaystyle \left(11,\,4\right)$, $\displaystyle \left(8,\,3\right)$, $\displaystyle \left(13,\,5\right)$, $\displaystyle \left(5,\,2\right)$, $\displaystyle \left(12,\,5\right)$, $\displaystyle \left(7,\,3\right)$, $\displaystyle \left(9,\,4\right)$, $\displaystyle \left(2,\,1\right)$, $\displaystyle \left(9,\,5\right)$, $\displaystyle \left(7,\,4\right)$, $\displaystyle \left(12,\,7\right)$, $\displaystyle \left(5,\,3\right)$, $\displaystyle \left(13,\,8\right)$, $\displaystyle \left(8,\,5\right)$, $\displaystyle \left(11,\,7\right)$, $\displaystyle \left(3,\,2\right)$, $\displaystyle \left(10,\,7\right)$, $\displaystyle \left(7,\,5\right)$, $\displaystyle \left(11,\,8\right)$, $\displaystyle \left(4,\,3\right)$, $\displaystyle \left(9,\,7\right)$, $\displaystyle \left(5,\,4\right)$, $\displaystyle \left(6,\,5\right)$, $\displaystyle \left(1,\,1\right)$]

drawfareytree((vector([2,1]),1))

(H+G+FT).show(aspect_ratio=1,figsize=10)

#Next we will try to have the rerpesentation of the continued fraction drawn as in this projective plane:

def square(origin,size, color=hue(0)):
    origin=vector(origin)
    v1,v2,v3,v4=\
        origin,\
        origin+vector([size,0]),\
        origin+vector([size,size]),\
        origin+vector([0,size])
    return line([myf(*v1),myf(*v2),myf(*v3),myf(*v4),myf(*v1)],color=color)
def squares(origin, size, repeats, direction, color=hue(0)):
    """
    origin is the begining of the squares,
    size is the number of squares to draw
    repeats is the numbers of squares to draw
    direction: 0 means towards the right; 1 upwards
    """
    G=Graphics()
    origin=vector(origin)
    if direction==0:
        dv=vector((size,0))
    else:
        dv=vector((0,size))
    for i in range(repeats):
        G+=square(origin+i*dv, size, color=color)
    return G
def MyQd(number):
    """
    This gives the quotients and the divisors in Euclid's Algorithm.
    It is not as efficient as the continued fraction routine,
    but it gives the devisors.
    """
    d0=number.numerator()
    D0=number.denominator()
    d=[d0]
    D=[D0]
    r=[]
    Q=[]
    while D0>0:
        r0=d0%D0
        r.append(r0)
        Q0=int(d0/D0)
        Q.append(Q0)
        d0=D0
        D0=r0
        D.append(D0)
        d.append(d0)
    D.pop()
    return Q,D
show(MyQd(68/157))

([$\displaystyle 0$, $\displaystyle 2$, $\displaystyle 3$, $\displaystyle 4$, $\displaystyle 5$], [$\displaystyle 157$, $\displaystyle 68$, $\displaystyle 21$, $\displaystyle 5$, $\displaystyle 1$])

def drawcontfrac(number, origin=(1,1),color=hue(0)):
    """
    Draws several rectangles giving a picture of the 
    continuous fraction expantion for number
    Assumes that number<1, and that number is rational.
    A recursive version is possible, but this one is iterative.
    The main idea is that for the fraction 68/157 whose sequence
    of quotients is Q=[0,2,3,4,5], and whose sequence of Divisors
    is D=[157, 68,21,5,1], the rectangle is drawn as follows:
    dv=[vector((1,0)),vector((0,1))]
    O=vector((0,0))
    #The number is 68/157
    G = squares(O,68,2,0)+\
        squares(O+2*68*(dv[0]),21,3,1)+\
        squares(O+2*68*(dv[0])+3*21*(dv[1]),5,4,0)+\
        squares(O+2*68*(dv[0])+3*21*(dv[1])+4*5*(dv[0]),1,5,1)
    G.show(aspect_ratio=1,gridlines=[srange(0,157,10),srange(0,68,10)],figsize=8)
    """
    dv=[vector((1,0)),vector((0,1))]
    Q,D=MyQd(number);#Q is the sequence of Quotients, and D is the sequence of Divisors
    Q.pop(0)
    D.pop(0);#We are assuming that the number is less than one
    G=Graphics()
    O=vector(origin);#A recursive version is also possilbe; but this solution is iterative
    for i in range(len(Q)):
        G+=squares(O,D[i],Q[i],(i%2),color=color)
        O+=D[i]*Q[i]*dv[i%2]
    return G
DC=drawcontfrac(68/157,origin=(1,1))
DC.show(aspect_ratio=1)

SS=squares((2,1),1,5,0)
(H+G+SS).show(aspect_ratio=1)
SS1=square((3,1),1)
(H+G+SS1).show(aspect_ratio=1)

T=Graphics()
for i in range(1,8):
    for j in range(0,i+1):
        T+=text(str((i,j)),myf(i,j),fontsize=20/i)
show(H+G+FT+T+squares((5,2),1,3,1),aspect_ratio=1,figsize=12)

DC=drawcontfrac(1/2,(2,1))+drawcontfrac(1/3,(3,1))+drawcontfrac(2/5,(5,2))+drawcontfrac(2/7,(7,2))
(H+G+FT+T+DC)