ABSTRACT

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 fast iterative method for computing values with high accuracy, for fixed and. The method is based on compound means and Padé approximations.

Keywords:Compound Means, Padé Approximation, Computation of Radicals, Iteration

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 that¹

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 method², 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.

References

NOTES

¹In the last line we assume the hypothesis that all multiplication algorithms satisfy, see .

²Not always Horner’s method is optimal, see . In those cases Householder’s and our algorithms are faster.