Calculate the power of a number quickly (MPFR library)

Question

I am using the GMP and MPFR libraries to work with large numbers and I need to calculate the power of a number quickly. The result of the potentiation will always be an integer, but the potency may or may not be a floating number. the GMP library calculates powers very quickly but does not accept floating powers (Using the mpz_pow_ui function), the MPFR library accepts floating powers but is extremely slow as it requires high precision to calculate integers correctly (using the function mpfr_pow). Is there any solution for this? How can GMP accept floating powers, or MPFR calculate whole numbers quickly (and correctly)?

//Ex:

mpz_pow_ui(mpz_power, base, 4790)  //Fast

// power = 4790.60
mpfr_pow(mpfr_power, base, power, MPFR_RNDN)  //Slow

Answer 1

It all depends on the order of magnitude of the numbers you are calculating with. The following assumes both the base and the exponent are positive non-complex numbers.

First note that exponentiation of an integer b by a decimal number like 4790.60 can be rewritten like this (just can't beleive I can't write Latex-style math equations in StackOverflow):

b ^ 4790.60 = b ^ (4790 + 0.60) = (b ^ 4790) * (b ^ 0.60)

Then, the first term (b ^ 4790) can clearly be calculated with GMP and results in an integer value.

The second term has an exponent smaller than 1, so its value will be smaller than b. If b is not a huge integer (say << FLINTMAX, the largest integer that can be represented consecutively in a floating point value, is 2^53 for double), then you can use the native double pow function to calculate it and safely round it to integer as desired, and then multiply by the first term. If b is a huge integer, you have too options: convert it to a double if it's within the double range and then use the double pow function (and perhaps loose some precision in the resulting integer with a speed gain), or you can use the mpfr_pow function to calculate this second term, but noting that it will be an smaller number than b, and so you can adjust MPFR precision in that calculation according to the exponent b: if b is closer to zero use a small precision, if b is close to 1 use a precision not too smaller than that of b.

On the end, all of it is a tradeoff between accuracy and speed, given the machine resource limits you are bound to.

Answer 2

The mpz_pow_ui function is fast since the computation can be done with multiplications. Note that mpfr_pow_ui with MPFR_RNDN should be faster if you don't need the exact integer result (which will be huge if the exponent is large), as the multiplications can be done with a smaller precision and rounded.

Unless the exponent is a not-too-large integer, mpfr_pow will be slower because it needs to compute a logarithm and an exponential, which are more complex than multiplications. I don't think that you can avoid it, even if you know that the result will be an integer. But if you know that the result will be an integer, you can compute it with an error less than 1/2. Thus, with an error analysis, you may be able to reduce the precision of the input and output variables, so that mpfr_pow will be faster.

Note: You are saying

the MPFR library accepts floating powers but is extremely slow as it requires high precision to calculate integers correctly (using the function mpfr_pow)

This is not true. MPFR will carefully chose the intermediate precisions to provide an answer with the required accuracy (though it may sometimes do wrong choices).

However, if the exact result is very close to a "machine number", there may be an issue due to the fact that MPFR needs to return the sign of the error (faithful rounding MPFR_RNDF instead of MPFR_RNDN could help, but it is not optimized yet for mpfr_pow): the Table Maker's Dilemma will occur, requiring internal computations in a higher precision. But if you have carefully chosen the precisions of the inputs, this is unlikely to occur in your case; indeed, reducing the precisions of the input will tend to add some error to the exact result, and this will move the exact result away from the expected integer (by exact result, I mean the exact result with the approximate inputs). You can also use some tricks. For instance, if you know that your integer is not a perfect square, then instead of computing x^y, you can compute x^y/2 (which does not correspond to a machine number, thus should not have an issue due to the Table Maker's Dilemma), then square the result.