Given (a, b) compute the maximum value of k such that a^{1/k} and b^{1/k} are whole numbers

Question

I'm writing a program that tries to find the minimum value of k > 1 such that the kth root of a and b (which are both given) equals a whole number.

Here's a snippet of my code, which I've commented for clarification.

int main()
{
    // Declare the variables a and b.
    double a;
    double b;
    // Read in variables a and b.
    while (cin >> a >> b) {

        int k = 2;

        // We require the kth root of a and b to both be whole numbers.
        // "while a^{1/k} and b^{1/k} are not both whole numbers..."
        while ((fmod(pow(a, 1.0/k), 1) != 1.0) || (fmod(pow(b, 1.0/k), 1) != 0)) {

        k++;

        }

Pretty much, I read in (a, b), and I start from k = 2 and increment k until the kth roots of a and b are both congruent to 0 mod 1 (meaning that they are divisible by 1 and thus whole numbers).

But, the loop runs infinitely. I've tried researching, and I think it might have to do with precision error; however, I'm not too sure.

Another approach I've tried is changing the loop condition to check whether the floor of a^{1/k} equals a^{1/k} itself. But again, this runs infinitely, likely due to precision error.

Does anyone know how I can fix this issue?

EDIT: for example, when (a, b) = (216, 125), I want to have k = 3 because 216^(1/3) and 125^(1/3) are both integers (namely, 5 and 6).

Answer 1

That is not a programming problem but a mathematical one:

if a is a real, and k a positive integer, and if a^(1./k) is an integer, then a is an integer. (otherwise the aim is to toy with approximation error)

So the fastest approach may be to first check if a and b are integer, then do a prime decomposition such that a=p₀^e₀ * p₁^e₁ * ..., where p_i are distinct primes.

Notice that, for a^1/k to be an integer, each e_i must also be divisible by k. In other words, k must be a common divisor of the e_i. The same must be true for the prime powers of b if b^1/k is to be an integer.

Thus the largest k is the greatest common divisor of all e_i of both a and b.

---

With your approach you will have problem with large numbers. All IIEEE 754 binary64 floating points (the case of double on x86) have 53 significant bits. That means that all double larger than 2⁵³ are integer.

The function pow(x,1./k) will result in the same value for two different x, so that with your approach you will necessary have false answer, for example the numbers 5⁵*2⁹⁰ and 3⁵*2¹²⁰ are exactly representable with double. The result of the algorithm is k=5. You may find this value of k with these number but you will also find k=5 for 5⁵*2⁹⁰-2⁴⁹ and 3⁵*2¹²⁰, because pow(5⁵*2⁹⁰-2⁴⁹,1./5)==pow(5⁵*2⁹⁰). Demo here

On the other hand, as there are only 53 significant bits, prime number decomposition of double is trivial.

Answer 2

Floating numbers are not mathematical real numbers. The computation is "approximate". See http://floating-point-gui.de/

You could replace the test fmod(pow(a, 1.0/k), 1) != 1.0 with something like fabs(fmod(pow(a, 1.0/k), 1) - 1.0) > 0.0000001 (and play with various such 𝛆 instead of 0.0000001; see also std::numeric_limits::epsilon but use it carefully, since pow might give some error in its computations, and 1.0/k also inject imprecisions - details are very complex, dive into IEEE754 specifications).

Of course, you could (and probably should) define your bool almost_equal(double x, double y) function (and use it instead of ==, and use its negation instead of !=).

As a rule of thumb, never test floating numbers for equality (i.e. ==), but consider instead some small enough distance between them; that is, replace a test like x == y (respectively x != y) with something like fabs(x-y) < EPSILON (respectively fabs(x-y) > EPSILON) where EPSILON is a small positive number, hence testing for a small L₁ distance (for equality, and a large enough distance for inequality).

And avoid floating point in integer problems.

Actually, predicting or estimating floating point accuracy is very difficult. You might want to consider tools like CADNA. My colleague Franck Védrine is an expert on static program analyzers to estimate numerical errors (see e.g. his TERATEC 2017 presentation on Fluctuat). It is a difficult research topic, see also D.Monniaux's paper the pitfalls of verifying floating-point computations etc.

And floating point errors did in some cases cost human lives (or loss of billions of dollars). Search the web for details. There are some cases where all the digits of a computed number are wrong (because the errors may accumulate, and the final result was obtained by combining thousands of operations)! There is some indirect relationship with chaos theory, because many programs might have some numerical instability.

Answer 3

As others have mentioned, comparing floating point values for equality is problematic. If you find a way to work directly with integers, you can avoid this problem. One way to do so is to raise integers to the k power instead of taking the kth root. The details are left as an exercise for the reader.