this doesn't work

var('c')
mat=matrix([[1,0],[0,c]])
xy0=vector([c,1])
r0=-mat*xy0
p0=r0
for k in range(2):
    la=r0*p0/(p0*mat*p0)
    xy1=xy0+la*p0
    print xy1
    if xy1[0]!=0:
        xy1[0]=xy1[0].factor()
    if (xy1[1])!=0:
        xy1[1]=xy1[1].factor()
    print xy1
    r1=(r0-la*mat*p0)
    p1=r1+r1*r1/(r0.dot_product(r0))*p0
    p1=p1.simplify_full()
    p0=p1
    r0=r1
    xy0=xy1
    print '--------------------'

but this works

var('c')
mat=matrix([[1,0],[0,c]])
xy0=vector([c,1])
r0=-mat*xy0
p0=r0
for k in range(2):
    la=r0*p0/(p0*mat*p0)
    xy1=xy0+la*p0
    print xy1
    xy1=xy1.simplify_full()   # <<<<<<<<< inserted this command
    if xy1[0]!=0:
        xy1[0]=xy1[0].factor()
    if (xy1[1])!=0:
        xy1[1]=xy1[1].factor()
    print xy1
    r1=(r0-la*mat*p0)
    p1=r1+r1*r1/(r0.dot_product(r0))*p0
    p1=p1.simplify_full()
    p0=p1
    r0=r1
    xy0=xy1
    print '--------------------'

v=vec([ax,ay])

should give

a*vec([x,y])

by a command. Which command?

The result equals zero. This ist wrong!