@Rouben Rostamian  

I greatly appreciate your help.

I'll look at the code that you have referenced above and try to understand that.

Thanks once again.

@Rouben Rostamian  

Thank you for your explanation.

With the invocation f([[-1, 0, 0], [2, 1, 0], [2, -1, 0], [3, 0, 2]]); L should be a list of 4 lists with 2 members each. Should that be the case the following assignments would be invalid, right?

T := abs(L[1][3]);
r := L[1][4];

You have correctly pointed out that "You would want J[1][2] to be nonzero, otherwise you will be dividing by zero.  In the example in the last line of your post is J[1][2] zero. "

Let me provide the complete program for you to see if it makes sense. It is supposed to generates Apollonian Gasket.

________________________________________________________________________

f:=proc(J)

L:=map((x)->[x[1], (x[2] + x[3] * I) / x[1] + 50 * (1+I) / abs(J[1][2])], J);

R:=Vector([L]);

T,r:=abs(L[1][3]),L[1][4];

A(L[1][5],L[2][6],L[3][7],L[1][8],L[2][9],L[3][10],R,T,r);

A(L[1][11],L[2][12],L[4][13],L[1][14],L[2][15],L[4][16],R,T,r);

A(L[1][17],L[3][18],L[4][19],L[1][20],L[3][21],L[4][22],R,T,r);

A(L[2][23],L[3][24],L[4][25],L[2][26],L[3][27],L[4][28],R,T,r);

plots[display](seq(plottools[circle]([Re(R[i][29]),Im(R[i][30])],abs(1/R[i][31])),i= 1..numelems(R))):

end proc:

A:=proc(a,b,c,i,j,k,R,E,F)

K:=i+k+j+2*sqrt(i*k+i*j+k*j);

if K>400*E then

   return;

end if;

C:=(a*i+c*k+b*j+2*sqrt(a*c*i*k+b*c*j*k+a*b*i*j))/K;

C2:=(a*i+c*k+b*j-2*sqrt(a*c*i*k+b*c*j*k+a*b*i*j))/K;

if evalf(abs(C-F))<1/E and not member([K,C],R) then

   R(N(R)+1):=[K,C];

   A(a,b,C,i,j,K,R,E,F);

   A(a,c,C,i,k,K,R,E,F);

   A(b,c,C,j,k,K,R,E,F);

end if:

if evalf(abs(C2-F))<1/E and not member([K,C2],R) then

   R(N(R)+1):=[K,C2];

   A(a,b,C2,i,j,K,R,E,F);

   A(a,c,C2,i,k,K,R,E,F);

   A(b,c,C2,j,k,K,R,E,F);

end if:

end proc:

________________________________________________________________________

Call: f([[-1, 0, 0], [2, 1, 0], [2, -1, 0], [3, 0, 2]]);