# Help with unstable code

This basically continues from a previous post: generation of certain matrices

I'm trying to draw circles in the plane by using fractional linear transformations.

Here is my code:

K = NumberField(x^2 + 2, 's')
OK = K.ring_of_integers()
def FLT(M,z):
"""takes the fractional linear transformation/Mobius transformation of a z in CC"""
if z == 'infinity' and M[1,0] != 0:
return M[0,0]/M[1,0]
elif z == 'infinity' and M[1,0] == 0:
return 'infinity'
elif M[1,0] != 0 and z == -M[1,1]/M[1,0]:
return 'infinity'
else:
return (M[0,0]*z+M[0,1])/(M[1,0]*z + M[1,1])


Now generate circles based on 3 image points of a matrix M with entries in the ring OK via FLT:

var('a','b','r','x','y')
func_cir=(x-a)**2+(y-b)**2==r**2
circles = set()
j = 0
while j < 10:
M = random_matrix(OK,2,algorithm = 'unimodular')
if FLT(M,-1) != 'infinity' and FLT(M,0) != 'infinity' and FLT(M,1) != 'infinity':
pta=[FLT(M,-1),FLT(M,-1)]
ptb=[FLT(M,0),FLT(M,0)]
ptc=[FLT(M,1),FLT(M,1)]
eq1 = func_cir.subs(x==pta).subs(y==pta)
eq2 = func_cir.subs(x==ptb).subs(y==ptb)
eq3 = func_cir.subs(x==ptc).subs(y==ptc)
sol = solve([eq1,eq2,eq3],a,b,r)
C = circle((sol.rhs(),sol.rhs()),sol.rhs())
j += 1


This second chunk of code throws the error "list index out of range" about 50% of the time. I thought that it may have been due to division by zero, but now I'm not so sure. I've also tried using a for loop with the exact same result. Thank you very much for the help!

edit retag close merge delete

Sort by » oldest newest most voted

The following alternative code is running without errors. The Möbius transform is slightly changed, and the $\sqrt 2=|s|$ was added. (I think it is also wanted.)

K.<s> = QuadraticField( -2 )
OK    = K.ring_of_integers()
IP    = ProjectiveSpace( 1, K )

points = ( -1, 0, 1 )
def lift_point( x ):
return vector( [ x, 1 ] )

def FLT( M, z ):
"""takes the fractional linear transformation/Mobius transformation of a z in CC"""
Mz = IP( ( M * lift_point(z) ).list() )
if Mz == IP( (1,0) ):
return 'infinity'
return Mz

def components( z ):
if z == 'infinity':
return None
return ( z, z*sqrt(2) )

var('a','b','r','x','y');

func_cir = (x-a)**2 + (y-b)**2 == r**2

def get_circles( number_of_circles=10, show=False ):
j = 0
circles  = []
while j < number_of_circles:

M = random_matrix( OK, 2, algorithm = 'unimodular' )
A, B, C = [ components( FLT( M, z ) ) for z in points ]

print ( "\n%s\nM is %s\nA,B,C are %s, %s, %s ."
% ( 30*'=', M.list(), A, B, C ) )
noneOccured = bool( None in ( A, B, C ) )
colinearityMatrix = matrix( 3,3, [ A, A, 1 ,
B, B, 1 ,
C, C, 1 ] ) if not noneOccured else None

if noneOccured:
print "!!! INFINITY... IGNORED !!!\n"

elif bool( colinearityMatrix.det() == 0 ):
# print ( "Are A,B,C colinear? %s" % colinearity )
print "!!! COLINEAR... IGNORED !!!\n"

else:
j += 1

eq1 = func_cir.subs( x==A ).subs( y==A )
eq2 = func_cir.subs( x==B ).subs( y==B )
eq3 = func_cir.subs( x==C ).subs( y==C )

sol = solve( [ eq1, eq2, eq3 ], a,b,r )

print "%s\n[%s]\nSolutions:\n%s" % ( 30*'-', j, sol )
c = circle( ( sol.rhs(),
sol.rhs() ),
abs( sol.rhs() ) )
circles.append( c )

if show:
plot_object = plot( [] )
for c in circles:
plot_object += plot( c )
plot_object.show()
return circles


Now the following calls did the job for me without errors:

get_circles( 1000 );


There is a verbose, aggressive print show in between, but the message is, there are no errors.

Let us also generate three circles, and show them:

get_circles( 3, show=True  );

more

Very nice code. I like how you made use of projective space. In doing so, you remind the user/reader that a FLT is just matrix multiplication. Thanks The error "list index out of range" occurs if solve returns an empty list. You can capture this exception with a try except statement:

    try:
C = circle((sol.rhs(),sol.rhs()),sol.rhs())
except:
pass

more