The Brouncker algorithm for the continued fractions of a root of a quadratic.

>	with(linalg):

To start, we create a matrix from  

>	A:=matrix(2,4,[1,0,1,0,0,1,0,-7]);

>

The function "fun" produces a polynomial in x,y from a 2 x 2 matrix. We need its values to calculate  steps.

>

>	fun:=m-> unapply(m[1,3]*x^2+2*m[1,4]*x*y+m[2,4]*y^2,x,y);

>	fun(A)(x,y);fun(A)(3,1);

The function ch(m,t,n) changes matrix m as follows:

If n=1 then it add t times row 1 to row 2.

If n=2 it adds t times row 2 to row 1.

Then it makes  a symmetric matrix operation on the last two columns.

>	ch:=proc(m,t,n) local tmp: tmp:=copy(m); if n=1 then  tmp:=addrow(tmp,1,2,t); tmp:=addcol(tmp,3,4,t); elif n=2 then  tmp:=addrow(tmp,2,1,t);tmp:=addcol(tmp,4,3,t);fi;op(tmp);end:

>

We start work on  .

We shall record the t-values as we proceed.

>	A:=matrix(2,4,[1,0,1,0,0,1,0,-7]);

We check the function at (1,1);

If it is less than 1,  then we try values at (t,1) and stop just before the values change sign.

>	fun(A)(1,1); fun(A)(2,1);fun(A)(3,1);

So set t=2 and n=1(because the value at (1,1) was less than 1.)

>	tvals:=[2];

>	A1:=ch(A,2,1);

Repeat with new matrix. Here, the value at (1,1) is bigger than 1 , so we try (1,t), until sign changes.

>	fun(A1)(1,1);fun(A1)(1,2);

So we set t=1 and n=2 (because the value at (1,1) was bigger  than 1.)

>	tvals:=[op(tvals),1];

>	A2:=ch(A1,1,2);

Repeat! Again, we get t=1, n=1.

>	fun(A2)(1,1);fun(A2)(2,1);

>	tvals:=[op(tvals),1];

>	A3:=ch(A2,1,1);

This is special!  

We have the f(1,1) value exactly equal to 1 . In this case, we set t=1 and n equal to the opposite of before, so n=2.

>	fun(A3)(1,1);

>	tvals:=[op(tvals),1];

>	A4:=ch(A3,1,2);

>	fun(A4)(1,1);fun(A4)(2,1);fun(A4)(3,1);fun(A4)(4,1);fun(A4)(5,1);

Now we get t=4, n=1.

>	tvals:=[op(tvals),4];

>	A5:=ch(A4,4,1);

>	seq(cat(A,i)=op(cat(A,i)),i=1..5);

Note that A1 and A4 have a  repeat of the second part .  So, now the pattern repeats!

>	convert(sqrt(7), confrac,8,ans7);

>

ans7;

>

Case of