Compound Means and Fast Computation of Radicals

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 k x with high accuracy, for fixed k ∈  and x + ∈  . The method is based on compound means and Padé approximations.


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 k x with high accuracy, for fixed k ∈  and x + ∈  .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 The convergence has order 2 d + .The order of convergence can be made arbitrarily large by the choice of d .But for larger values of d 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.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 , P Q be two means.Given two positive numbers u t, , put : , , : , , If these two sequences converge to a common limit then this limit is denoted by ( ) , : lim lim .
The function  is called compound mean.The best known application of compound means is Gauss' arithmetic-geometric mean [4] ( ) 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 Q P, be continuous means such that at least one of them is strict.Then the compound mean P Q ◊ exists and is continuous.

Properties
We call a mean P homogeneous if for every , , c t u , , .

P ct cu cP t u =
All power means are homogeneous.If two means are homogeneous then their compound mean is also homogeneous.
On the other hand, there is only one mean R satisfying the functional equation , , , , .

R P t u Q t u R t u =
Easy proofs of these facts can be found in [5].
Example 1.Take the arithmetical mean ( ) and the harmonic mean ( ) 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 The Babylonian method is in fact Newton's method used to solve the equation 2 0 t x − =.Using Newton's method to solve the equation 0 with a quadratical convergence to lim

Our Method
In the present paper we will proceed similarly as in Example 1.For a fixed integer 2 k ≥ and a positive real number x we will find a sequence of approximations converging fast to k x .
We need the following lemma.

Lemma 1. Let function p satisfy
and let ( ) Assume that p and q satisfy (2.1) strictly if 1 t ≠ .Then the function q satisfies ( ) ( ) ( ) ( ) and ( ) 1 w q + is bounded.Let P and Q be the homogeneous means corre- sponding to traces p and q , respectively.Then the compound mean P Q ◊ exists and its convergence has Take the mean ( ) This mean is strict, continuous and homogeneous and it has the property ( ) . We will construct two means Let s ∈  and let ( ) p t be the Padé approximation of the function k t of order [ ] The exact formula for , j j e f will be derived in Lemma 15.In Lemma 20 we will prove that p satisfies ( ) ( ) ( ) for every t + ∈  .Hence it is a trace of a strict homogeneous mean , , and its trace is As in Definition 2, denote the sequences given by the compound mean P Q ◊ by , n n a b , starting with 0 1 a = and 0 b x = .From (3.3) we obtain that ( ) ( ) . So the iterations of the compound mean Note that we don't have to compute the sequence n b .According to (3.1) and Lemma 1 the convergence of the sequence (3.4) to its limit k x has order 2 1 s + .

( )
M N 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 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 N is not very large.The complexity of division of two N-digit numbers differs from the complexity of multiplication only by some multiplicative constant D .Hence the complexity of division is ( ) DM N .Analysis in [9] shows that this constant can be as small as 7 2 D = .
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 r ∈  and that ( ) For every 0 and assume that 1) for every T the image set 2) there is ( ) ≥ for every ( ) From the monotonicity of f we have for every N ( ) ( ) Let 0 δ > .Then (4.1) implies ( ) = ∞ .From this and properties 1 and 2 we deduce that there exists a number ( ) ≥ for every ( ) and every m n ≤ .This implies for every ( ) Passing to the limit 0 δ → implies the result.□ Note that all the above mentioned functions Now we compute the complexity ( ) V N of algorithm computing k x to within N digits.The functions ( ) M N and ( ) V N have asymptotically the same order as N → ∞ , see for instance Theorem 6.3 in [3].Hence all fast algorithms for computing k x differ only in the asymptotic constants.
Let the algorithm for computing k x • performs Z multiplications of two N-digit numbers before the iterations, • performs A multiplications and B divisions of two long numbers in every step, • has order of convergence r .
The accuracy to within N digits is necessary only in the last step.In the previous step we need accuracy only to within N r digits and so on.Hence The error term

( )
O N corresponds to additions of N-digit numbers.This and Lemma 2 imply that 1

Complexity of Newton's Method
Newton's method (2.3) has order 2 and in every step it performs ( ) 1 In the last line we assume the hypothesis that all multiplication algorithms satisfy ( ) ( ) [12]. and 1 division.So the complexity is For the choice of multiplication and division algorithms with 1 ω = and

Complexity of Householder's Method
Consider Householder's method (1.1) applied to the equation ( ) An easy calculation (see [11] for instance) leads to iterations ( ) ( ) ( ) , 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 For the choice of multiplication and division algorithms with 1 The optimal value of d which minimizes the complexity is in this case ( ) ( )

Complexity of Our Method
Given s ∈  , our method (3.4) has order 1 2 + s , before the iterations it performs 1 − s multiplications of N-digit numbers (evaluation of 2 , , s x x  ) and in every step it performs , then evaluations of numerator and denominator by Horner's method 2 , and then the final multiplication) and 1 division.So the complexity of our algorithm is where 2 Not always Horner's method is optimal, see [13].In those cases Householder's and our algorithms are faster.
For the choice of multiplication and division algorithms with 1 The optimal value of s which minimizes the complexity is in this case ( )

Examples
Example 2. Compare the algorithms for computation of 14 x .For 14 k = we have ( ) and, according to (4.4), the optimal value of s for our algorithm is with convergence of order 3.For computation of N digits of 14 x we need to perform ( ) ( ) ( ) Newton's method has order 2 and for computation of N digits it needs ( ) ( ) operations.Hence our method saves 16% of time compared to Newton's method.For Householder's method the optimal value is 1 d = and it leads to the same iterations (5.1) as our method.Example 3. Compare the algorithms for computation of 179 x .For 179 k = we have ( ) and, according to (4.4), the optimal value of s for our algorithm is with convergence of order 5.For computation of N digits of 179 x we need to perform operations.Newton's method This method has order 3 and for computation of N digits it needs Hence our method saves 20% of time compared to Newton's method and saves 0.6% of time compared to Householder's method.
with convergence of order 7.For computation of N digits of 1234567890133 x we need to perform ( ) ( ) ( ) For Householder's method the optimal value is 3 d = and it leads to iterations 0 : 1, a = This method has order 5 and for computation of N digits it needs ( ) ( ) ( ) Newton's method has order 2 and for computation of N 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.

Proofs
In this section we will prove that function p defined by (3.1) satisfies inequalities (3.2).For the sake of brevity, will use the symbol u  for the set

Combinatorial Identities
First we need to prove several combinatorial identities.Our notation will be changed in this subsection.Here, n will be a variable used in mathematical induction, k 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 . This new function is then used for easier evaluation of sums containing ( ) we immediately obtain , , , ) For fixed k n > we have and hence ( ) ( ) Similarly for 0 k = and Then (6.1), (6.2), (6.3) and (6.4) imply ( ) ) and ( ) Similarly for A n < and B n = we obtain From this and (6.6) we obtain that , : , , .
This and the fact that

S n A F n k A F n k A T n k A F n k n A A T n k n A A F n k A T n k A n
, , , .

B n B n S n A B S n A B n
A n This with (6.9) implies the result for every n B > .□ For 0 , , , : , , . (6.10) From this, (6.10) and the polynomial identity

11)
For n = A the sum ( )

S n A contains only one nonzero summand for k A
= and we have ( ) Equation (6.11) implies that 4 S will not change for greater n .
From this, (6.14) and the polynomial identity ! 2 S n B has the same value for greater n . ) From this and (6.16) we obtain ∑ ∑ (6.17) S n B has the same value for greater n .□

Formulas
Now the symbol k again has its original meaning.For See for instance Theorem 159 in [15].□ Now we prove a technical lemma that we need later.Lemma 11.
s j s j hk hk j s j j s j s j i hk hk i j i s j i Proof.On both sides of (6.18) there are polynomials of degree j in variable k .Therefore it suffices to prove the equality for k → ∞ and for j values with 1, , m s j s = − +  we obtain on the left-hand side of (6.18) For such numbers m 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 1 i m ≥ + , the second product on the last line is equal to zero for 1 i m j s ≤ + − − .For other values of i both products contain only positive terms.Hence , for i j > we have ( ) and for i m > we have ( )

RHS . LHS
. 1 This expression, multiplied by j k , is a polynomial of degree The second product is equal to zero for i m j > − , therefore the summations ends for i m j = − .Then Lemma 5 implies    .S n r A i B s s j r s j s i r i j i s j s r s r s i j s j r s j i j s j r s s r s r s i i r i r s j s j h r i j s j s h r i , , , 0    The inequality is strict, since for 1 j ≥ the main bracket in the numerator is greater than the main bracket in the denominator.□ Now we prove that function

0d
∈  be a parameter of the method.When solving the equation ( ) 0 f t = the iterations converging to the solution are means are continuous and strict.There is a known inequality between power means ,

Example 4 .
Compare the algorithms for computation of 1234567890133 x .For 1234567890133 k = the exact value of ( ) k σ is not known.We assume that ( ) 59 k σ = .According to (4.4), the optimal value of s for our algorithm is 3 s = .The iterations of our algorithm are 0

2 k ∈  and 1 ig
are the coefficients of Taylor's polynomial of the function k t .Lemma 10.For 2 k ∈  and

(
(6.22) and (6.21) follow.Putting : k k = − into (6.21) and applying binomial theorem in the same way we obtain (6.20).□On both sides there are polynomials of degree s in variable k .Therefore it suffices to prove the Lemma 10, Lemma 12 and Lemma 13 imply 1