Transformed Bakshali Square-Root Algorithm with Higher Orders of Convergence

1. Introduction

The time tested method of computing the square root of a number is Newton’s method which is based on the Newton-Raphson iterative scheme,

$$x_{n+1}=x_n - {\displaystyle \frac{\mathrm{f}(x_n)}{\mathrm{f} \text{‘}\left(x_n\right)}}----------$$

(1)

While finding the square root of a number N, this becomes,

$$x_{n+1}={\displaystyle \frac{1}{2}}\left(x_n+{\displaystyle \frac{N}{x_n}}\right)----------$$

(2)

Iteration of this formula will give better and better approximations to$\sqrt{N}$​. The historically important Babylonian and Heron’s methods are identical to Newton’s method as given in eq. (2).[1,2,3,4]

Newton’s method is quadratically converging; after each iteration the number of correct digits doubles. The method is computationally very simple; it has just one addition, one division and one multiplication by half in each iteration. With a suitable choice for initial guess, Newton’s method converges fairly fast to the square root, with minimum utilization of computer resources.[2]

The ancient Bakshali manuscript found in 1881 in the village of Bakshali, approximately 80 km northwest of Peshawar gives an interesting formula for finding the square root of a number. The Bakshali formula is,

$$\sqrt{N}=a+{\displaystyle \frac{r}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2\left(a+\frac{r}{2a}\right)}}----------$$

(3)

In eq.(3), a is such that its square a² is close to N. Then r is the remainder such that N=a²+r.The effectiveness of this formula is demonstrated in the Bakshali manuscript itself by finding the square root of a number with 12 significant digit accuracy. It is truly amazing that calculations with such precision were carried out in ancient times (c.400 CE).[5,6,7,8,9,10]

Bailey and Borwein showed that the Bakshali formula given in eq. (3) can be cast as an iterative formula.[7]

$${\mathrm{P}}_{\mathrm{n}}= {\displaystyle \frac{\mathrm{N}-{\mathrm{x}}_{\mathrm{n}}^2}{2{\mathrm{x}}_{\mathrm{n}}}}$$

$${ x}_{n+1 }= x_n{+\mathrm{P}}_{\mathrm{n}}-{\displaystyle \frac{{\mathrm{P}}_{\mathrm{n}}^2}{2(x_n\mathrm{}+{\mathrm{P}}_{\mathrm{n}\mathrm{}})\mathrm{}\mathrm{}}}$$

They have shown that this is a quartically convergent scheme. Each iteration approximately quadruples the number of correct digits. They have also shown that this follows naturally from the fact that the scheme is equivalent to performing two consecutive iterations of the Newton–Raphson equation, which is quadratically convergent.

Let us look at the process through which the Bakshali formula was most probably arrived at. It is proposed to find the square root of a positive real number N greater than 1. The first step is to guess a square number a² close to N, such that N=a²+r, where r is the remaining part of N once a²is removed from it. The first approximation to the square root of N is,

$$x_1=a+{\displaystyle \frac{\mathrm{r}}{2\mathrm{a}}}------------$$

(4)

This approximation was probably known right from the early part of the 2^nd millennium BCE. The Babylonian tablet YBC 7289 shows the value of$\sqrt{2}$​ with some remarkable accuracy for that time.[3]

Eq.(4) is equivalent to the Newton’s method formula given in eq.(2).

$$x_1=a+{\displaystyle \frac{N-a^2}{2a}} (\because r=N-a^2)$$

$$=a+{\displaystyle \frac{N}{2a}}-{\displaystyle \frac{a}{2}}$$

$$={\displaystyle \frac{1}{2}}(a+{\displaystyle \frac{N}{a}})$$

If we take x₁ as given in eq.(4), as the approximate value of square root of N, the error can be calculated as,

$$e_1=x_1^2-N$$

$$={\left(a+{\displaystyle \frac{r}{2a}}\right)}^2-N$$

$$={ a}^2+r+{\left({\displaystyle \frac{r}{2a}}\right)}^2-N$$

$$={\left({\displaystyle \frac{r}{2a}}\right)}^2$$

Let us define a new variable y₁ which will be applied as a correction to x₁ to reduce the error.

$$y_1={\displaystyle \frac{e_1}{2x_1}}={\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}$$

The new, corrected, approximation to the square root is,

$$x_2=x_1-y_1$$

$$ie x_2= a+{\displaystyle \frac{r}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}------------$$

(5)

When x₂ is considered as the square root of N, the error e₂ can be calculated as,

$${e_{2 =} x}_2^2-\mathrm{N}$$

$${=\left[a+{\displaystyle \frac{r}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}\right]}^2-\mathrm{N}$$

$$={\left(a+{\displaystyle \frac{r}{2a}}\right)}^2-{\left({\displaystyle \frac{r}{2a}}\right)}^2+{\left[{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}\right]}^2-\mathrm{N}$$

$$={a^2+r+\left({\displaystyle \frac{r}{2a}}\right)}^2-{\left({\displaystyle \frac{r}{2a}}\right)}^2+{\left[{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}\right]}^2-\mathrm{N}$$

$$={\left[{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}\right]}^2={\left[{\displaystyle \frac{e_1}{2(a+\frac{r}{2a})}}\right]}^2$$

${ e}_{2 }\cong {\left[{\displaystyle \frac{e_1}{2\sqrt{\mathrm{N}}}}\right]}^2 (a+{\displaystyle \frac{r}{2a}}$ is our first approximation to $\sqrt{N }$)

As N is a positive real number greater than 1, e₂ is always less than e₁. The error decreases, and it justifies the correction term y₁.

Equation (5) is exactly the Bakshali formula given in eq.(3). The process of error reduction we followed here can be extended to obtain a recursive algorithm to find better and better approximations to the square root of any positive real number. Since the process is in the spirit of the Bakshali formula, we call it the Bakshali algorithm.

2. Bakshali Algorithm

To find the square root of N, we first write N=a²+r, where a² is a square number close to N. It is best to choose a² as the greatest square number less than N. At each step, we find the approximation to $\sqrt{N }$ and compute the error. Based on this error, a correction is applied in the next step, to obtain a better approximation to $\sqrt{N } .$The algorithm n proceeds as follows:

The expression for $\sqrt{N }$ obtained in the n^th step is,

$$x_n=a+{\displaystyle \frac{r}{2a}}-y_1-y_2-y_3-\text{- ---------}{-y}_{n-1}$$

The error in the n^th stage can be calculated as,

$${\mathrm{e}}_{\mathrm{n}}=x_n^2-\mathrm{N}$$

$$={\left(x_{n-1}-y_{n-1}\right)}^2$$

$$=x_{n-1}^2-{2x}_{n-1}{y}_{n-1 }+y_{n-1}^2-\mathrm{N}$$

$$=\left(x_{n-1}^2-\mathrm{N}\right)-y_{n-2}^2+y_{n-1}^2$$

$$=\left(e_{n-1}\right)-y_{n-2}^2+y_{n-1}^2$$

$$\mathrm{i}.\mathrm{e}. {\mathrm{e}}_{\mathrm{n}}=y_{n-1}^2..................$$

(6)

But $y_{n-1}={\displaystyle \frac{y_{n-2}^2}{{2x}_{n-1}}} = {\displaystyle \frac{e_{n-1}}{{2x}_{n-1}}}$

$$\mathrm{i}.\mathrm{e}. {e_n=\left[{\displaystyle \frac{e_{n-1}}{{2x}_{n-1}}}\right]}^2$$

As the denominator is $\cong 2\sqrt{\mathrm{N}}$ ,

$${ e}_n\cong {\left({\displaystyle \frac{e_{n-1}}{2\sqrt{N}}}\right)}^2...............$$

(7)

We can see that e_n is always less than e_n-1. The error is progressively decreasing. This ensures the convergence of the algorithm to the true value of $\sqrt{\mathrm{N}}$The first three approximations to $\sqrt{\mathrm{N}}$ in Table1 actually give the following expressions:

$$x_1=a+ {\displaystyle \frac{r}{2a}}$$

$${ x}_2=a+{\displaystyle \frac{r}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}},$$

$$x_3=a+{\displaystyle \frac{r}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}-{\displaystyle \frac{{\left[\frac{{\left(\frac{r}{2a}\right)}^2}{2\left(a+\frac{r}{2a}\right)}\right]}^2}{2 \left[a+\frac{r}{2a}-\frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}\right]}}$$

The following points may be noted about the x_n​ values given in the algorithm.

• Each x_n​ gives a better approximation to the square root than given by the previous expression, x_n−1​.
• Each expression for x_n​ can be used as an iterative formula for getting better and better values of $\sqrt{N}$

3.

Step1 is the same as Newton’s method whose quadratic convergence is well established. Step2 is the Bakshali formula whose quartic(4-fold) convergence is also established previously.

4.

In general the order of convergence gets multiplied by 2 as we go from one step to the next. The order of convergence should be 2ⁿ if we use the formula x_n​ in the n^th step for iteration.

Computationally, using these expressions as iteration formulas don’t give any advantage over the standard Newton’s method. For example, if we use the level 3 formula (x₃) for iteration, we get convergence 2³=8. Each iteration will multiply the number of correct digits by 8. But within each iteration the number of operations will be the same as 3 repeated applications of Newton’s method.

In order to reduce the computational complexity of operations, we now recast the Bakshali algorithm in a different manner. Especially we want to limit the number of divisions to a minimum, as the computational cost of division is much more than other mathematical operations, addition, subtraction, and multiplication.

3. Transformed Bakshali Algorithm

We can recast the first three square root values of the Bakshali algorithm of Table 1 in the following manner.

$$x_1=a+{\displaystyle \frac{r}{2a}}={\displaystyle \frac{{2a}^2+r}{2a}}={\displaystyle \frac{z_1}{2a}}$$

$$x_2={ x}_1-y_1$$

$$={\displaystyle \frac{z_1}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2x_1}}$$

$$={\displaystyle \frac{z_1}{2a}}-{\displaystyle \frac{r^2}{{4a}^2\times 2\frac{z_1}{2a}}}$$

$$={\displaystyle \frac{z_1}{ 2a}}-{\displaystyle \frac{r^2}{2a\times 2z_1}}$$

$$={\displaystyle \frac{{2z}_1^2-r^2}{2a\times 2z_1}}={\displaystyle \frac{z_2}{2a\times 2z_1}}$$

$$x_3={ x}_2-y_2$$

$$= {\displaystyle \frac{z_2}{2a\times {2z}_1}}-{\displaystyle \frac{y_1^2}{2x_2}}$$

$$= {\displaystyle \frac{z_2}{2a\times {2z}_1}}-{\displaystyle \frac{r^4}{{\left(2a\times {2z}_1\right)}^2\times \frac{2z_2}{2a\times {2z}_1}}}$$

$$= {\displaystyle \frac{z_2}{2a\times {2z}_1}}-{\displaystyle \frac{r^4}{2a\times 2z_1{\times 2z}_2}}$$

$$= {\displaystyle \frac{{2 z}_2^2-r^4}{2a\times 2z_1{\times 2z}_2}}= {\displaystyle \frac{z_3}{2a\times 2z_1{\times 2z}_2}}$$

In this manner we can recast the whole Bakshali algorithm. The transformed Bakshali algorithm thus obtained is given below.

This algorithm can be implemented for various orders of convergence. Iteration of operations in step 1 will give x₁,​ with 2-fold (quadratic) convergence. Iterating the operations of steps 1 and 2 together will give x₂ values with 4-fold (quartic) convergence (x₁​ need not be computed). Iteration of steps 1, 2, and 3 together (without computing x₁​ and x₂) will give x₃ values with 8-fold convergence. Putting n steps together in a loop and iterating, we will get 2ⁿ-fold convergence directly for x_n.

4. Implementation of Transformed Bakshali Algorithm

Let us write down the 3-step iteration scheme for 8-fold convergence.

To calculate $\sqrt{N}$​, assumeN=a²+r

Input:   N, a

r(1) = N−a², a(1) = a

For i=1 to n

y(1) = r(i)

z(1) = 2a(i)²+y(1)

t(1) = 2a

y(2) = y(1)²

z(2) = 2z(1)²−y(2)

t(2) = t(1)×2z(1)

y(3) = y(2)²

z(3) = 2z(2)²−y(3)

t(3) = t(2)×2z(2)

x(i) = ${\displaystyle \frac{\mathrm{z}\left(3\right)}{\mathrm{t}\left(3\right)}}$a(i+1) = x(i)

r(i+1) = N−a(i+1)²

Each iteration of this loop ensures 8-fold increase in correct decimal places in the square root of N. There is only one division operation, and the multiplication operations are mostly squaring. This should result in reduced execution times, especially when the number of precision digits becomes large.

The following table gives the results of running the above program implemented in Python language. The results are for a randomly chosen number N=37957. The accuracy of a root obtained is verified by looking at the remainder r=N−a². The number of correct digits after the decimal point in the root is almost the same as the number of leading zeroes in r, the difference being only one digit or two at the most. The correctness of the result was verified by comparing the digits of the square root obtained, with standard library sqrt(N) values.

In just 7 iterations a square root with more than 3.6 million correct digits are extracted by this level-3 Bakshali algorithm. As the iterations progress, the convergence ratio approaches very fast to the number 8. But the ratio is always a little more than 8 due to the presence of few extra digits.

The standard Newton–Raphson method takes 21 iterations to reach the same level of precision. The precision achieved in a single Bakshali iteration will be achieved by 3 iterations of Newton–Raphson. Each iteration of Newton–Raphson implements a division. As the iterations move ahead, these divisions become multi-million digit operations. The Bakshali algorithm of 8-fold convergence condenses the entire algebraic expansion (equivalent to three Newton–Raphson iterations) into a single unified mathematical structure. It achieves 8-fold jump in accuracy using only one final division operation at the very end.

The Bakshali algorithm eliminates the intermediary division operation via the introduction of squaring operations. Modern high performance math libraries use fast multiplication algorithms. Bakshali algorithm is ideally suited to execute these fast squaring operations to achieve high-speed convergence, especially at multi-million digit accuracy.

The remainder term r=N−a² plays a crucial role in the Bakshali algorithm. It is calculated at the end of each iteration and becomes the basis for the calculations in the next iteration. The remainder term provides the exact measurement of the algorithm’s real accuracy. Because the remainder is explicitly calculated at the end of each iteration, it acts as a data anchor which pulls the algorithm back to the correct path if unwanted effects try to limit its efficiency. Normally floating- point calculations across a number of steps cause tiny round-off errors which can compound and degrade the performance of an algorithm. The explicit calculation of the remainder r at the end of each iteration eliminates the possibility of compounding and amplification of errors.

5. Error analysis

Let us see how the 3-stage Bakshali algorithm diminishes error.

Let the initial value of the square root be a.

Let e₀ be the absolute error in our current value of the root.

$$\sqrt{N}=a+e_0$$

$$\mathrm{N}={a^2+2ae}_0+e_0^2$$

The algorithm operates by sequentially transforming this starting error e₀ in to higher powers through three structural squaring levels.

Level 1

Remainder r is calculated.

$$r=N-a^2=2a{e}_0+e_0^2$$

$$y_{1 }=r=2a{e}_0+e_0^2$$

t₁ = 2a

$${\mathrm{z}}_1=2{\mathrm{a}}^2+\mathrm{r}=2{\mathrm{a}}^2+2a{e}_0+e_0^2$$

Level 2

Remainder is squared.

$${\mathrm{y}}_2=y_1^2={\left(2a{e}_0+e_0^2\right)}^2$$

$$=4a^2{e}_0^2+4a{e}_0^3+e_0^4$$

Next, the first correction envelope z₂ is computed.

$${ z}_2=2z_1^2-y_2$$

$$=2{\left(2a^2+2a{e}_0+e_0^2\right)}^2-(4a^2{e}_0^2+4a{e}_0^3+e_0^4$$

)

$$=2(4a^2+4a^2{e}_0^2+e_0^4+8a^3{e}_0+4a{e}_0^3+4a^2{e}_0^2)-\left(4a^2{e}_0^2+4a{e}_0^3+e_0^4\right)$$

$$=8a^4+16a^3{e}_0+12a^2{e}_0^2+4a{e}_0^3+e_0^4$$

t₂ = 2a × 2z₁

$$=2\mathrm{a}\text{×}2\left(2a^2+2a{e}_0+e_0^2\right)$$

$$=8a^3+8a^2{e}_0+4a{e}_0^2$$

Level 3

$${\mathrm{y}}_3=y_2^2={\left(4a^2{e}_0^2+4a{e}_0^3+e_0^4\right)}^2$$

$$=16a^4{e}_0^4+32a^3{e}_0^5+24a^2{e}_0^6+8a{e}_0^7+e_0^8$$

$$z_2^2={\left(8a^4+16a^3{e}_0+12a^2{e}_0^2+4a{e}_0^3+e_0^4\right)}^2$$

$$=64a^8+256a^7{e}_0+448a^6{e}_0^2+448a^5{e}_0^3+288a^4{e}_0^4+128a^3{e}_0^5+ 40a^2{e}_0^6+8a{e}_0^7+e_0^8$$

$${\mathrm{z}}_3=2z_2^2-y_3$$

$${=128 a}^8+512a^7{e}_0+896a^6{e}_0^2+896a^5{e}_0^3+576a^4{e}_0^4+256a^3{e}_{0 }^5+ 80a^2{e}_0^6+16a{e}_0^7+2e_0^8-\left(16a^4{e}_0^4+32a^3{e}_0^5+24a^2{e}_0^6+8a{e}_0^7+ e_0^8\right)$$

$${ =128 a}^8+512a^7{e}_0+896a^6{e}_0^2+896a^5{e}_0^3+560a^4{e}_0^4+224a^3{e}_{0 }^5+ 56a^2{e}_0^6+8a{e}_0^7+e_0^8$$

t₃ = t₂× 2z₂

$$=\left(16a^3+16a^2{e}_0+8a{e}_0^2\right)\times (8a^4+16a^3{e}_0+12a^2{e}_0^2+4a{e}_0^3+e_0^4)$$

$${=128a}^7+384a^6{e}_0+512a^5{e}_0^2+384a^4{e}_0^3+176a^3{e}_0^4+48a^2{e}_{0 }^5+ 8a{e}_0^6$$

Square root R = ${\displaystyle \frac{Z_3}{t_3}}$

$$={\displaystyle \frac{{ 128 a}^8+572a^7{e}_0+896a^6{e}_0^2+896a^5{e}_0^3+560a^4{e}_0^4+224a^3{e}_{0 }^5+ 56a^2{e}_0^6+8a{e}_0^7+e_0^8}{{ 128a}^7+384a^6{e}_0+512a^5{e}_0^2+384a^4{e}_0^3+176a^3{e}_0^4+48a^2{e}_{0 }^5+ 8a{e}_0^6}}$$

$$=a+e_{0 }+{\displaystyle \frac{{ e}_0^8}{{ 128a}^7+384a^6{e}_0+512a^5{e}_0^2+384a^4{e}_0^3+176a^3{e}_0^4+ 48a^2{e}_{0 }^5+ 8a{e}_0^6}}$$

$\cong a+e_{0 }$+${\displaystyle \frac{e_0^8}{{128a}^7}}$(Neglecting the lower order terms in the denominator)

$$=\sqrt{N }$$

The error in this is ${\displaystyle \frac{e_0^8}{{128a}^7}}$There is order of 8 reduction in error. Due to the factor 128a⁷ in the denominator there are some extra digits of precision.

6. Transformed Bakshali Algorithm with 16-fold convergence

To further demonstrate the validity of the transformed Bakshali scheme, we have also implemented a Level-4 Bakshali algorithm, which ensures 16-fold convergence. Here, the first four levels of the transformed Bakshali algorithm of Table 2 are taken as a unit and iterated. The results are shown in Table 4.

Table 2. Transformed Bakshali Algorithm.

Step	y(n)	z(n)	t(n)	x(n)
1	y1= r	z1= 2a2+y1	t1= 2a	x1=${\displaystyle \frac{z_1}{t_1}}$
2	y2= ${\mathrm{y}}_1^2$	z2= $2 {\mathrm{z}}_1^{2 }$-y2	t2= t1×2z1	x2=${\displaystyle \frac{z_2}{t_2}}$
3	y3 = ${\mathrm{y}}_2^2$	z3= $2 {\mathrm{z}}_{2 }^2$-y3	t3= t2×2z2	x3=${\displaystyle \frac{z_3}{t_3}}$
4	y4 =${\mathrm{y}}_3^2$	z4=$2 {\mathrm{z}}_{3 }^2$-y4	t4= t3×2z3	x4=${\displaystyle \frac{z_4}{t_4}}$
.	----------	-------------	-------------	-------
.	----------	-------------	-------------	-------
n	yn=${\mathrm{y}}_{n-1}^2$	zn=$2 {\mathrm{z}}_{n-1}^2$-yn	tn= tn-1×2zn-1	xn=${\displaystyle \frac{z_n}{t_n}}$

Table 2. Transformed Bakshali Algorithm.

Step	y(n)	z(n)	t(n)	x(n)
1	y1= r	z1= 2a2+y1	t1= 2a	x1=${\displaystyle \frac{z_1}{t_1}}$
2	y2= ${\mathrm{y}}_1^2$	z2= $2 {\mathrm{z}}_1^{2 }$-y2	t2= t1×2z1	x2=${\displaystyle \frac{z_2}{t_2}}$
3	y3 = ${\mathrm{y}}_2^2$	z3= $2 {\mathrm{z}}_{2 }^2$-y3	t3= t2×2z2	x3=${\displaystyle \frac{z_3}{t_3}}$
4	y4 =${\mathrm{y}}_3^2$	z4=$2 {\mathrm{z}}_{3 }^2$-y4	t4= t3×2z3	x4=${\displaystyle \frac{z_4}{t_4}}$
.	----------	-------------	-------------	-------
.	----------	-------------	-------------	-------
n	yn=${\mathrm{y}}_{n-1}^2$	zn=$2 {\mathrm{z}}_{n-1}^2$-yn	tn= tn-1×2zn-1	xn=${\displaystyle \frac{z_n}{t_n}}$

Table 3. Square root extraction of 37957 by Level-3 (8-fold convergence) Bakshali algorithm.

Iteration No.	Estimated numberof precision digits in the square root	Magnitudeof remainder(only5 significant digits are shown with normal rounding off)	Convergence ratio
1	P1=15	r =6.2706×10-17	_____
2	P2=115	r =1.1155×10-116	${\displaystyle \frac{\mathrm{P}2}{\mathrm{P}1}}$=7.6666
3	P3=881	r =3.2077×10-882	${\displaystyle \frac{\mathrm{P}3}{\mathrm{P}2}}$ =7.6608
4	P4=7088	r =6.0273×10-7089	${\displaystyle \frac{\mathrm{P}4}{\mathrm{P}3}}$ =8.0454
5	P5=56741	r =9.3654×10-56743	${\displaystyle \frac{\mathrm{P}5}{\mathrm{P}4}}$ =8.0052
6	P6=453972	r =3.1825×10-453973	${\displaystyle \frac{\mathrm{P}6}{\mathrm{P}5}}$ =8.0007
7	P7=3631816	r =5.6596×10-3631817	${\displaystyle \frac{\mathrm{P}7}{\mathrm{P}6}}$ =8.00008

Table 3. Square root extraction of 37957 by Level-3 (8-fold convergence) Bakshali algorithm.

Iteration No.	Estimated numberof precision digits in the square root	Magnitudeof remainder(only5 significant digits are shown with normal rounding off)	Convergence ratio
1	P1=15	r =6.2706×10-17	_____
2	P2=115	r =1.1155×10-116	${\displaystyle \frac{\mathrm{P}2}{\mathrm{P}1}}$=7.6666
3	P3=881	r =3.2077×10-882	${\displaystyle \frac{\mathrm{P}3}{\mathrm{P}2}}$ =7.6608
4	P4=7088	r =6.0273×10-7089	${\displaystyle \frac{\mathrm{P}4}{\mathrm{P}3}}$ =8.0454
5	P5=56741	r =9.3654×10-56743	${\displaystyle \frac{\mathrm{P}5}{\mathrm{P}4}}$ =8.0052
6	P6=453972	r =3.1825×10-453973	${\displaystyle \frac{\mathrm{P}6}{\mathrm{P}5}}$ =8.0007
7	P7=3631816	r =5.6596×10-3631817	${\displaystyle \frac{\mathrm{P}7}{\mathrm{P}6}}$ =8.00008

Table 4. √37957 calculated using Bakshali Level-4 algorithm with 16-fold convergence.

Table 4. √37957 calculated using Bakshali Level-4 algorithm with 16-fold convergence.

The algorithm extracts the square root with more than 2.8 million correctdigits in just five iterations. The convergence ratio approaches 16 very fast.

7. Conclusions

• Following the logic of Bakshali square root formula, a general Bakshali algorithm for square root extraction is presented in Table 1.
• From this general Bakshali algorithm, a transformed Bakshali algorithm is developed in Table 2.
• Different levels of the transformed Bakshali algorithm can be used as iterative schemes for finding the square root of a positive real number.
• The effectiveness of the transformed Bakshali algorithm is demonstrated for Level-3 and Level-4, which shows 8-fold and 16-fold convergences respectively. Millions of digits of precision are achieved in just a small number of iterations.
• It is expected that the error-correcting, self-regulating and fast converging transformed Bakshali algorithm will become the basis for an efficient square-root engine which can produce square roots with arbitrarily large digits of precision with more efficient use of computational resources.

Equations

$x_{n+1}=x_n - {\displaystyle \frac{\mathrm{f}(x_n)}{\mathrm{f} \text{‘}\left(x_n\right)}}$----------------------(1)

$x_{n+1}={\displaystyle \frac{1}{2}}\left(x_n+{\displaystyle \frac{N}{x_n}}\right)$----------------------(2)

$\sqrt{N}=a+{\displaystyle \frac{r}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2\left(a+\frac{r}{2a}\right)}}$ --------------------(3)

$x_1=a+{\displaystyle \frac{\mathrm{r}}{2\mathrm{a}}}$ --------------------(4)

$x_2=a+{\displaystyle \frac{r}{2a}}-{\displaystyle \frac{{\left(\frac{r}{2a}\right)}^2}{2(a+\frac{r}{2a})}}$ --------------------(5)

• e_n = $y_{n-1}^2$

-------------------- (6)

$e_n\cong {\left({\displaystyle \frac{e_{n-1}}{2\sqrt{N}}}\right)}^2$ ----------------------(7)