1. Introduction
In last decades, several algorithms were developed for fast evaluation of some elementary functions with very large arguments, for example for multiplication of million-digit integers. The present paper introduces a new iterative method for computing values 
with high accuracy, for fixed 
and
.
The best-known method used for computing radicals is Newton’s method used to solve the equation
Newton’s method is a general method for numerical solution of equations and for particular choice of the equation it can lead to useful algorithms, for example to algorithm for division of long numbers. This method converges quadratically. Householder [1] found a generalization of this method. Let 
be a parameter of the method. When solving the equation 
the iterations converging to the solution are
(1.1)
The convergence has order
. The order of convergence can be made arbitrarily large by the choice of
. But for larger values of 
it is necessary to perform too many operations in every step and the method gets slower.
The method (3.4) presented in this paper involves compound means. It is proved that this method performs less operations and is faster than the former methods.
Definition 1. A function 
is called mean if for every 
A mean 
is called strict if
A mean 
is called continuous if the function 
is continuous.
A known class of means is the power means defined for 
by
for 
we define
The most used power means are the arithmetic mean
, the geometric mean 
and the harmonic mean
. All power means are continuous and strict. There is a known inequality between power means
see e.g. [2] . From this one directly gets the inequality between arithmetic mean and geometric mean. For other classes of means see e.g. [3] .
Taking two means, one can obtain another mean by composing them by the following procedure.
Definition 2. Let 
be two means. Given two positive numbers
, put
. If these two sequences converge to a common limit then this limit is denoted by
The function 
is called compound mean.
The best known application of compound means is Gauss’ arithmetic-geometric mean [4]
(1.2)
Iterations of the compound mean then give a fast numerical algorithm for computation of the elliptic integral (1.2).
Matkowski [5] proved the following theorem on existence of compound means.
Theorem 1. Let 
be continuous means such that at least one of them is strict. Then the compound mean 
exists and is continuous.
2. Properties
We call a mean 
homogeneous if for every 
All power means are homogeneous. If two means are homogeneous then their compound mean is also homogeneous.
Homogeneous mean 
can be represented by its trace
Conversely, every function 
with property
(2.1)
represents homogeneous mean
Theorem 2. If the compound mean 
exists then it satisfies the functional equation
(2.2)
On the other hand, there is only one mean 
satisfying the functional equation
Easy proofs of these facts can be found in [5] .
Example 1. Take the arithmetical mean 
and the harmonic mean
. The arithmetic-harmonic mean 
exists by Theorem 1 and Theorem 2 implies that
. Hence the iterations of the arithmetic-harmonic mean 
can be used as a fast numerical method of computation of
. This leads to a well known Babylonian method
with a quadratical convergence to
.
The Babylonian method is in fact Newton’s method used to solve the equation
. Using Newton’s method to solve the equation 
leads to iterations
(2.3)
with a quadratical convergence to
.
3. Our Method
In the present paper we will proceed similarly as in Example 1. For a fixed integer 
and a positive real number 
we will find a sequence of approximations converging fast to
.
We need the following lemma.
Lemma 1. Let function 
satisfy 
for 
and let 
be bounded. Let
. Assume that 
and 
satisfy (2.1) strictly if
. Then the function 
satisfies
for 
and 
is bounded. Let 
and 
be the homogeneous means corresponding to traces 
and
, respectively. Then the compound mean 
exists and its convergence has order
.
Proof. The assumption implies that
Then
hence 
for 
and 
is bounded.
The compound mean 
exists by Theorem 1. Let 
and 
be the iterations of
. To find the order of convergence put
. Then
and
. □
Take the mean
. This mean is strict, continuous and homogeneous and it has the property
. We will construct two means 
such that
.
Let 
and let 
be the Padé approximation of the function 
of order 
around
,
(3.1)
The exact formula for 
will be derived in Lemma 15. In Lemma 20 we will prove that 
satisfies
(3.2)
for every
. Hence it is a trace of a strict homogeneous mean
Relation 
and Theorem 2 imply that
(3.3)
and its trace is
Inequalities (3.2) imply that 
is a strict homogeneous mean too.
As in Definition 2, denote the sequences given by the compound mean 
by
, starting with 
and
. From (3.3) we obtain that
hence 
and
. So the iterations of the compound mean 
are
(3.4)
Note that we don’t have to compute the sequence
.
According to (3.1) and Lemma 1 the convergence of the sequence (3.4) to its limit 
has order
.
4. Complexity
Let 
denote the time complexity of multiplication of two N-digit numbers. The classical algorithm of multiplication has asymptotic complexity
. But there are also algorithms with asymptotic complexity
or
see Karatsuba [6] , Schönhage and Strassen [7] or Fürer [8] , respectively. The fastest algorithms have large asymptotic constants, hence it is better to use the former algorithms if the number 
is not very large.
The complexity of division of two N-digit numbers differs from the complexity of multiplication only by some multiplicative constant
. Hence the complexity of division is
. Analysis in [9] shows that this constant can be as small as
.
We will denote by 
the minimal number of multiplications necessary to compute the power
. See [10] for a survey on known results about the function
.
Before the main computation of complexity we need this auxiliary lemma.
Lemma 2. Assume that 
and that 
is a function such that for some 
the function
is nondecreasing with
(4.1)
For every 
put 
and assume that 1) for every 
the image set 
is bounded and 2) there is 
with 
for every
.
Then
Proof. From the monotonicity of 
we have for every 
(4.2)
Let
. Then (4.1) implies
. From this and properties 1 and 2 we deduce that there exists a number 
such that 
for every 
and every
. This implies for every 
(4.3)
Inequalities (4.2) and (4.3) yield
Passing to the limit 
implies the result. □
Note that all the above mentioned functions 
satisfy all assumptions of Lemma 2 with
, 
, 
and
, respectively.
Now we compute the complexity 
of algorithm computing 
to within 
digits. The functions 
and 
have asymptotically the same order as
, see for instance Theorem 6.3 in [3] . Hence all fast algorithms for computing 
differ only in the asymptotic constants.
Let the algorithm for computing 
• performs 
multiplications of two N-digit numbers before the iterations• performs 
multiplications and 
divisions of two long numbers in every step• has order of convergence
.
The accuracy to within 
digits is necessary only in the last step. In the previous step we need accuracy only to within 
digits and so on. Hence
The error term 
corresponds to additions of N-digit numbers. This and Lemma 2 imply that1
4.1. Complexity of Newton’s Method
Newton’s method (2.3) has order 2 and in every step it performs 
multiplications (evaluation of
) and 1 division. So the complexity is
where
For the choice of multiplication and division algorithms with 
and 
we have
4.2. Complexity of Householder’s Method
Consider Householder’s method (1.1) applied to the equation
. Let
. An easy calculation (see [11] for instance) leads to iterations
where 
and 
are suitable constants. The method has order
, before the iterations it performs 
multiplications of N-digit numbers (evaluation of
) and in every step it performs 
multiplications (evaluation of
, then evaluations of numerator and denominator by Horner’s method, and then the final multiplication) and 1 division. So the complexity of Householder’s algorithm is
where
For the choice of multiplication and division algorithms with 
and 
we have
The optimal value of 
which minimizes the complexity is in this case
4.3. Complexity of Our Method
Given
, our method (3.4) has order
, before the iterations it performs 
multiplications of N-digit numbers (evaluation of
) and in every step it performs 
multiplications (evaluation of
, then evaluations of numerator and denominator by Horner’s method2, and then the final multiplication) and 
division. So the complexity of our algorithm is
where
For the choice of multiplication and division algorithms with 
and 
we have
(4.4)
The optimal value of 
which minimizes the complexity is in this case
5. Examples
Example 2. Compare the algorithms for computation of
. For 
we have 
and, according to (4.4), the optimal value of 
for our algorithm is
. Padé approximation of the function 
around 
is
Hence the iterations of our algorithm are
(5.1)
with convergence of order 3. For computation of 
digits of 
we need to perform
operations.
Newton’s method
has order 2 and for computation of 
digits it needs 
operations. Hence our method saves 16% of time compared to Newton’s method.
For Householder’s method the optimal value is 
and it leads to the same iterations (5.1) as our method.
Example 3. Compare the algorithms for computation of
. For 
we have 
and, according to (4.4), the optimal value of 
for our algorithm is
. Padé approximation of the function 
around 
is
Hence the iterations of our algorithm are
with convergence of order 5. For computation of 
digits of 
we need to perform
operations.
Newton’s method
has order 2 and for computation of 
digits it needs
operations.
For Householder’s method the optimal value is 
and it leads to iterations
This method has order 3 and for computation of 
digits it needs
operations.
Hence our method saves 20% of time compared to Newton’s method and saves 0.6% of time compared to Householder’s method.
Example 4. Compare the algorithms for computation of
. For 
the exact value of 
is not known. We assume that
. According to (4.4), the optimal value of 
for our algorithm is
. The iterations of our algorithm are
with convergence of order 7. For computation of 
digits of 
we need to perform
operations.
For Householder’s method the optimal value is 
and it leads to iterations
This method has order 5 and for computation of 
digits it needs
operations.
Newton’s method
has order 2 and for computation of 
digits it needs (assuming that
)
operations.
Hence our method saves 35% of time compared to Newton’s method and saves 7% of time compared to Householder’s method.
6. Proofs
In this section we will prove that function 
defined by (3.1) satisfies inequalities (3.2). For the sake of brevity, will use the symbol 
for the set
.
6.1. Combinatorial Identities
First we need to prove several combinatorial identities. Our notation will be changed in this subsection. Here, 
will be a variable used in mathematical induction, 
will be a summation index, and 
will be additional parameters. The change of notation is made because of easy application of the following methods based on [14] . For a function 
we will denote its differences by
Given a function
, there is some function 
satisfying some relation between 
and
. This new function is then used for easier evaluation of sums containing
. Recall that
for negative integer
.
For 
and 
with 
put
Lemma 3. For every 
satisfying 
we have
Proof. From the polynomial identity
we immediately obtain
(6.1)
For fixed 
we have 
and hence
. Thus
(6.2)
Similarly for 
and 
we have
(6.3)
and
(6.4)
Then (6.1), (6.2), (6.3) and (6.4) imply
(6.5)
If 
and 
then 
for 
and 
for
. Hence the only nonzero summand in 
is the one for 
and
(6.6)
Similarly for 
and 
we obtain
From this and (6.6) we obtain that
Equation (6.5) implies that 
will not change for greater
. □
For 
with 
and 
put
Lemma 4. For every 
with 
Proof. From the polynomial identity
we obtain for 
that
This and the fact that 
imply
□
For 
with 
and 
put
Lemma 5. For every 
with 
Proof. From the polynomial identity
we obtain for 
and 
that
(6.7)
For 
we have
. This and (6.7) imply
(6.8)
Lemma 4 for 
implies
hence
and
(6.9)
For 
we have 
and (6.8) yields
This with (6.9) implies the result for every
. □
For 
with 
put
Lemma 6. For every 
with 
Proof. From the polynomial identity
we obtain for 
and 
that
(6.10)
For 
we have
. From this, (6.10) and the polynomial identity
we obtain
(6.11)
For n = A the sum 
contains only one nonzero summand for 
and we have
. Equation (6.11) implies that 
will not change for greater
. □
For 
with 
put
Lemma 7. For every 
with 
Proof. From the polynomial identity
we obtain for 
and 
(6.12)
For 
and 
we have
and
From this and (6.12) we obtain
(6.13)
Lemma 6 for 
yields
Equation (6.13) implies that 
has the same value for greater
. □
For 
and 
put
Lemma 8. For every 
Proof. From the polynomial identity
we obtain
(6.14)
For 
we have
. From this, (6.14) and the polynomial identity
we obtain
(6.15)
For 
we have
Equation (6.15) implies that 
has the same value for greater
. □
For 
and 
with 
put
Lemma 9. For every 
with 
Proof. From the polynomial identity
we obtain for 
and 
with 
(6.16)
For 
and 
we have
and
From this and (6.16) we obtain
(6.17)
For 
Lemma 8 implies
Equation (6.17) implies that 
has the same value for greater
. □
6.2. Formulas
Now the symbol 
again has its original meaning.
For 
and 
put
The numbers 
are the coefficients of Taylor’s polynomial of the function
.
Lemma 10. For 
and 
Proof. See for instance Theorem 159 in [15] . □
Now we prove a technical lemma that we need later.
Lemma 11. For
, 
and 
(6.18)
Proof. On both sides of (6.18) there are polynomials of degree 
in variable
. Therefore it suffices to prove the equality for 
and for 
values 
with
.
For 
we immediately have
For 
with 
we obtain on the left-hand side of (6.18)
For such numbers 
we have
Therefore we obtain on the right-hand side of (6.18)
The first product on the last line is equal to zero for
, the second product on the last line is equal to zero for
. For other values of 
both products contain only positive terms. Hence
This implies that
For 
we have
, for 
we have 
and for 
we have
. Then Lemma 3 implies
Hence equality (6.18) follows for 
with
. □
For
, 
and 
put
For 
and 
put
(6.19)
We will prove that (6.19) is Padé approximation of
.
Lemma 12. For
, 
the numbers 
and 
satisfy the system of equations
Proof. Lemma 11 implies
□
Lemma 13. For
, 
the numbers 
and 
satisfy the system of equations
Proof. For 
we obtain on the left-hand side
This expression, multiplied by
, is a polynomial of degree 
in variable
. Therefore it suffices to prove the equality for 
values 
with
. For 
we have 
and hence the whole expression is equal to zero. For 
we obtain
The second product is equal to zero for
, therefore the summations ends for
. Then Lemma 5 implies
□
Lemma 14. Function 
is the Padé approximation of the function 
of order 
around
.
Proof. For 
Lemma 10, Lemma 12 and Lemma 13 imply
The result follows. □
Now we find the coefficients 
and 
of the Padé approximation
Lemma 15. For every
, 
and 
(6.20)
(6.21)
Proof. First we prove (6.21). From (6.19) we obtain
Binomial theorem then implies that
Thus equality (6.21) is equivalent to
(6.22)
On both sides of (6.22) there are polynomials of degree 
in variable
. Therefore it suffices to prove the equality for every 
with
,
.
Lemma 7 implies
hence (6.22) and (6.21) follow.
Putting 
into (6.21) and applying binomial theorem in the same way we obtain (6.20). □
6.3. Bounds
By Lemma 14 the function 
is Padé approximation, hence we know its properties in the neighbourhood of 1. Here we find global bounds for 
that are necessary for functionality of our algorithm.
We need another two technical lemmas.
Lemma 16. For 
and 
Proof. On both sides there are polynomials of degree 
in variable
. Therefore it suffices to prove the equality for every 
with
,
.
Lemma 9 implies
This implies the result. □
Lemma 17. For 
and 
with 
Proof. Put
. Then
Lemma 16 implies
□
Now we find lower and upper bounds for the function
.
Lemma 18. For
, 
and 
we have
.
Proof. Put
(6.23)
Lemma 12 and Lemma 13 imply
(6.24)
From (6.23) we obtain for 
We have
. Lemma 17 implies
(6.25)
From (6.23), (6.24) and (6.25) we obtain
Lemma 15 implies that
hence
□
Lemma 19. For every
, 
and 
we have
.
Proof. Directly from the definition of 
and 
we obtain
(6.26)
The inequality is strict, since for 
the main bracket in the numerator is greater than the main bracket in the denominator. □
Now we prove that function 
satisfies inequalities (3.2).
Lemma 20. For every
, 
and 
Proof. The proof splits into four cases.
1) For 
Lemma 18 implies that
.
2) Lemma 15 implies that 
for every
. Hence
(6.27)
Using this and the first case we obtain that for 
we have
.
3) For 
Lemma 19 implies that
.
4) Using this and (6.27) we obtain for 
that
. □
Acknowledgements
The author would like to thank to professor Andrzej Schinzel for recommendation of the paper [14] and also to author’s colleagues Kamil Brezina, Lukáš Novotný and Jan Štĕpnička for checking the results. Publication of this paper was supported by grant P201/12/2351 of the Czech Science Foundation, by grant 01798/2011/RRC of the Moravian-Silesian region and by grant SGS08/PrF/2014 of the University of Ostrava.
NOTES
1In the last line we assume the hypothesis that all multiplication algorithms satisfy
, see .
2Not always Horner’s method is optimal, see . In those cases Householder’s and our algorithms are faster.