Ошибка в алгоритме

Question

Закодил алгоритм проверки числа на то, является ли оно степенным, работает хорошо, но проверка не засчитывает 1 тест. В первой строке количество чисел которые нужно проверить, от 1 до 10 включительно, потом сами числа от 1 до 10 в 9 степени включительно. Основной принцип проверки - раскладывание числа на простые множители, и если нет ни одного простого множителя, который повторяется, то число не степенное.

#include <iostream>
#include <cmath>

using namespace std;
bool f(int b);

int main() {
    int n, b;
    cin >> n;
    for (int i = 0; i < n; i++)
    {
        cin >> b;
        if (step(b)) cout << "YES" << endl;
        else cout << "NO" << endl;
    }
}

bool f(int b)
{
    long double T=2;
    for (double i = 2; T>1.9; i++)
    {
    T = floor((pow(b, 1.0 / i)*1000000))/1000000;
    if ((double)T == (int)T) return true;
    }
    return false;
}

Answer 1

Если степенные числа это те, которые можно получить из некоторого числа умножением на себя хотя бы один раз.

Вот достаточно быстрый алгоритм основанный на свойствах логарифма:

double log2Values(int x, int y) {
    return log(x) / log(y);
}

bool f(int n)
{
    for (int i = 2; i * i <= n; ++i) {
        if (fabs(log2Values(n, i) - round(log2Values(n, i))) < 1e-7)
            return true;
    }
    return false;
}

А так у вас если ввести число 72, то выведется "YES", что не правильно.

Код выше хоть и простой, но очень долгий и не точный. На числах уже больше 1 миллиона может выдавать неправильные значения.

Пример быстрого и эффективного кода предоставил @StanislavVolodarskiy.

Answer 2

Предупреждение: хотя этот способ работает, есть способ лучше.

Разлагаем на простые множители и подсчитываем НОД их степеней в разложении.

Если НОД не единица, то он есть степень в которую можно возвести (неизвестное нам) число и получить n.

Несколько примеров:

  3 = 3^1       -> НОД(1)    = 1 -> обычное
  4 = 2^2       -> НОД(2)    = 2 -> степенное
  6 = 2^1 * 3^1 -> НОД(1, 1) = 1 -> обычное
 36 = 2^2 * 3^2 -> НОД(2, 2) = 2 -> степенное
 72 = 2^3 * 3^2 -> НОД(3, 2) = 1 -> обычное
144 = 2^4 * 3^2 -> НОД(4, 2) = 2 -> степенное

Разложение сделано перебором делителей i с двойки до тех пор пока i^2 <= n. Если делитель найден (n % i == 0), то делим n на него пока делится. Процесс продолжается дальше. Получается цикл в котором возможный делитель i растёт, а делимое n убывает. Цикл не самый оптимальный: например проверять чётные делители кроме двух не имеет смысла.

Получив разложение (например 2^4 * 3^2) считаем НОД показателей степеней (НОД(4, 2) = 2 в нашем примере). Сами множители нас не интересуют.

// g++ -Wall -Wextra -pedantic -std=c++11 powerful_numbers.cpp

#include <iostream>
#include <numeric>

int gcd(int a, int b) {
    while (b > 0) {
        int t = b;
        b = a % b;
        a = t;
    }
    return a;
}

bool powerful(int n) {
    int g = 0;
    int j = 2;
    while (n > 1) {
        int divisor = n;
        for (; j * j <= n; ++j) {
            if (n % j == 0) {
                divisor = j;
                break;
            }
        }

        int power = 0;
        while (n % divisor == 0) {
            n /= divisor;
            ++power;
        }

        g = gcd(g, power);
        if (g == 1) {
            return false;
        }
    }
    return true;
}


void check(int n) {
    std::cout << n << ' ' << powerful(n) << '\n';
}

int main() {
    for (int i = 0; i < 10; ++i) {
        check(i);
    }
    check(36);
    check(72);
    check(144);
}

Результат:

Answer 3

Ошибка

Ваш способ не работает из-за ошибок округления и неправильного применения конструкции (int). И long double не лучше чем double в этой ситуации.

Пусть n = 125.

Выражение pow(125, 1.0/3.0) вычисляется в 4.999999999999999. Так может быть: для функции возведения в степень нет гарантий точного результата, даже если этот результат представим в виде double без искажений. :(

Выражение floor(4.999999999999999 * 1000000) / 1000000 вычисляется в 4.999999. Здесь сознательно занижается точность, но это не беда.

Выражение (int)4.999999 вычисляется в 4. Ошибка тут. Малая до сих пор ошибка округления стала катастрофой.

Округление можно поправить с (int) на round и этот случай вычислится правильно. Зато перестанут правильно вычислятся такие вещи:

2**29 = 536870912
pow(536870911, 1.0/29.0) = 1.9999999998715416
pow(536870912, 1.0/29.0) = 2.0
pow(536870913, 1.0/29.0) = 2.000000000128458

Исправление

Пусть p - степень которую мы проверяем. Вычислим корень степени p и расстояние от корня до ближайшего целого:

double root = std::pow(n, 1. / p);
double error = std::abs(root - std::round(root));

Даже если у числа есть точный корень степени p ошибка может не быть равной нулю. Если у корня нет точного корня степени p можно надеятся что ошибка будет достаточно велика. Проверить можно на всех числах n. Несколько часов на миллиард тестов и результат готов. Если число степенное, то ошибка не больше 3.41061e-13. Если число не степенное, то ошибки для всех степеней не меньше 1.28458e-10. Ура! Вот работающий вариант:

bool powerful(int n) {
    int powers[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
    for (int p : powers) {
        double root = std::pow(n, 1. / p);
        double error = std::abs(root - std::round(root));
        if (error < 1e-11) { 
            return true;
        }
    }
    return false;
}

Эта функция быстрее всех: десять раз извлекаем корень и проверяем расстояние до целого. Всё. Единственное что трудно объяснить, так это почему она работает и что такое 1e-11.

Другое решение

Чтобы найти ошибки в предыдущем разделе нужно знать какие числа степенные а какие нет. Здесь приводится другой метод выяснения степенности числа.

Из целого числа 0 <= n <= 10^9 извлечём корни простой степени от 2 до 29. Полученные корни округлим и возведём обратно в соответствующие степени. Возводить нужно точно. Результаты сравним с n. Если есть совпадение - n степенное число.

Простых степеней достаточно, так если n - степенное число по составной степени, то n - степенное и по любому делителю степени: n = r^(ab) = (r^a)^b.

bool powerful(int n) {
    int powers[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
    for (int p : powers) {
        double root = std::round(std::pow(n, 1. / p));
        if (pow(root, p) == n) {
            return true;
        }
    }
    return false;
}

Функция возведения в степень самая простая, но точная. double выбран так как для целых возможно переполнение, которое не хочется контролировать (позже выяснилось что переполнения именно в этой задаче не бывает никогда). Вещественные числа в нужном нам диапазоне перемножаются точно. Иначе говоря, что double, что long long будут вести себя одинаково в нашей задаче:

double pow(double k, int p) {
    double r = 1;
    for (int i = 0; i < p; ++i) {
        r *= k;
    }
    return r;
}

Этот способ также проходит проверку на e-olymp. Он проще разложения на простые, работает быстрее.

Answer 4

Просто и незатейливо, и даже не слишком оптимально - но с гарантией и проходимостью :) e-olymp:

#include <iostream>
#include <cmath>

using namespace std;

int qpow(int x, unsigned int e)
{
    int res = 1;
    for(;e;e>>=1)
    {
        if (e&1) res *= x;
        x *= x;
    }
    return res;
}

int main()
{
    int N;
    cin >> N;
    for(int i = 0; i < N; ++i)
    {
        int n;
        cin >> n;
        double x = log(n);
        for (int j = 2; j*j <= n; ++j)
        {
            int p = int(x/log(j)+0.5);
            if (qpow(j,p) == n) { n = 0; break; }
        }
        cout << (n ? "NO" : "YES") << endl;
    }
}

Забавно, что решение побыстрее

#include <iostream>
#include <cmath>

using namespace std;

int qpow(int x, unsigned int e)
{
    int res = 1;
    for(;e;e>>=1)
    {
        if (e&1) res *= x;
        x *= x;
    }
    return res;
}

unsigned int pp[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 };

bool isPow(int n)
{
    double x = log(n);
    for(auto p: pp)
    {
        int a = int(exp(x/p)+0.5);
        if (qpow(a,p) == n) return true;
    }
    return false;
}

int main()
{
    int N;
    cin >> N;
    for(int i = 0; i < N; ++i)
    {
        int n;
        cin >> n;
        cout << (isPow(n) ? "YES" : "NO") << endl;
    }
}

проходит за ту же миллисекунду, что и первое... :)

БЛОГ НА HUSL

Ошибка в алгоритме

Ответы (4 шт):

Ошибка

Исправление

Другое решение