Учитывая (a, b) вычислить максимальное значение k так, чтобы a ^ {1 / k} и b ^ {1 / k} были целыми числами

Question

Я пишу программу, которая пытается найти минимальное значение k> 1 так, чтобы k-й корень a и b (которые оба заданы) равнялся целому числу.

Вот фрагмент моего кода, который я прокомментировал для пояснения.

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++;

        }

В значительной степени я читаю в (a, b) и начинаю с k = 2 и увеличиваю k до тех пор, пока k-е корни a и b не будут конгруэнтны 0 0 1 (это означает, что они делятся на 1 и, таким образом, целые числа).

Но цикл работает бесконечно. Я пытался исследовать, и я думаю, что это может быть связано с ошибкой точности; однако я не слишком уверен.

Другой подход, который я попробовал, - это изменить условие цикла, чтобы проверить, равен ли пол ^ {1 / k} самому ^ {1 / k}. Но опять же, это происходит бесконечно, вероятно из-за ошибки точности.

Кто-нибудь знает, как я могу решить эту проблему?

РЕДАКТИРОВАТЬ: например, когда (a, b) = (216, 125), я хочу, чтобы k = 3, потому что 216 ^ (1/3) и 125 ^ (1/3) оба целые числа (а именно, 5 и 6).

Answer 1

Как уже упоминалось, сравнение значений с плавающей точкой на равенство проблематично.Если вы найдете способ работать непосредственно с целыми числами, вы можете избежать этой проблемы.Один из способов сделать это - повысить целые числа до степени k вместо взятия корня k.Подробности оставлены в качестве упражнения для читателя.

Answer 2

Это не проблема программирования, а математическая:

, если a - действительное число, а k - положительное целое число, а если a^(1./k) - целое число, то a является целым числом.(в противном случае цель состоит в том, чтобы играть с ошибкой аппроксимации)

Таким образом, самый быстрый подход может состоять в том, чтобы сначала проверить, являются ли a и b целыми числами, затем выполнить простое разложение такой, что a = p ₀ ^{e ₀} * p ₁ ^{e ₁} * ..., где p _i - различные простые числа.

Обратите внимание, что для ^{1 / k} должно быть целое число, каждый e _i также должно делиться на k.Другими словами, k должно быть общим делителем e _i .То же самое должно быть верно для простых степеней b, если b ^{1 / k} должно быть целым числом.

Таким образом, наибольшее k является наибольшим общим делителем всех e _i обоих a и b.

---

При вашем подходе у вас будут проблемы с большими числами.Все двоичные числа с плавающей запятой IIEEE 754 (случай double на x86) имеют 53 значащих бита.Это означает, что все double больше 2 ⁵³ являются целыми числами.

Функция pow(x,1./k) приведет к одному и тому же значению для двух разных x, так что при вашем подходе вам потребуетсяиметь ложный ответ, например, числа 5 ⁵ * 2 ⁹⁰ и 3 ⁵ * 2 ¹²⁰ точно представимы двойным.Результат работы алгоритма k=5.Вы можете найти это значение k с этим числом, но вы также найдете k=5 для 5 ⁵ * 2 ⁹⁰ -2 ⁴⁹ и 3 ⁵ * 2 ¹²⁰ , потому что pow (5 ⁵ * 2 ⁹⁰ -2 ⁴⁹ , 1. / 5) == пау (5 ⁵ * 2 ^{90 * +1087 *).Демонстрация здесь}

С другой стороны, поскольку имеется только 53 значащих бита, разложение по простому числу double тривиально.

Answer 3

Плавающие числа , а не математические действительные числа.Расчет "приблизительный".См. http://floating -point-gui.de /

. Вы можете заменить тест fmod(pow(a, 1.0/k), 1) != 1.0 чем-то вроде fabs(fmod(pow(a, 1.0/k), 1) - 1.0) > 0.0000001 (и поиграть с различными такими ? вместо 0.0000001;см. также std :: numeric_limits :: epsilon , но используйте его осторожно, поскольку pow может привести к некоторой ошибке в его вычислениях, а 1.0/k также добавляет неточности - детали очень сложны, погрузитесь в IEEE754 спецификации).

Конечно, вы можете (и, вероятно, должны) определить вашу bool almost_equal(double x, double y) функцию (и использовать ее вместо ==, и использовать ее отрицание вместо !=).

Как правило, никогда не проверяйте плавающие числа на равенство (т. Е. ==), но вместо этого учитывайте достаточно малое расстояние между ними;то есть замените тест, например, x == y (соответственно x != y) на что-то вроде fabs(x-y) < EPSILON (соответственно fabs(x-y) > EPSILON), где EPSILON - небольшое положительное число, следовательно, тестирование на небольшое L ₁ расстояние (для равенства и достаточно большое расстояние для неравенства).

И избегайте чисел с плавающей запятой в целочисленных задачах.

На самом деле, прогнозировать или оценивать точность с плавающей запятой очень сложно.Возможно, вы захотите рассмотреть такие инструменты, как CADNA .Мой коллега Франк Ведрин является экспертом по статическим анализаторам программ для оценки числовых ошибок (см., Например, его презентацию TERATEC 2017 на Fluctuat ).Это сложная тема для исследования, см. Также статью Д.Моннио , подводные камни проверки вычислений с плавающей запятой и т. Д.

И ошибки с плавающей запятой внекоторые случаи стоили человеческих жизней (или потеря миллиардов долларов).Поиск в Интернете для деталей.В некоторых случаях все цифры вычисленного числа неверны (потому что ошибки могут накапливаться, и окончательный результат был получен путем объединения тысяч операций)!Существует некоторая косвенная связь с теорией хаоса , поскольку многие программы могут иметь некоторую числовую нестабильность .