Вперед, на поиски палиндромов 2

    Не так давно прочитал на хабре статью grey_wolfs «Вперед, на поиски палиндромов» о решении и оптимизации любопытной конкурсной задачки с весьма лаконичной формулировкой:

    «The decimal number 585 is 1001001001 in binary. It is palindromic in both bases. Find n-th palindromic number». Или, по-русски: «Десятичное число 585 в двоичной системе счисления выглядит как 1001001001. Оно является палиндромом в обеих системах счисления. Найдите n-й подобный палиндром».

    Автор начал свой увлекательный рассказ с самого тривиального решения перебором и проверкой всех чисел от 1 до N и вычислительной сложностью O(N * log(N)), где N — максимальное проверяемое число. Множитель log(N) необходим, т.к. для каждого проверяемого числа совершается несколько действий, зависящих от количества его разрядов.

    Первая же рабочая оптимизация, а именно замена перебора всех чисел на перебор только десятичных палиндромов, кардинально улучшила вычислительную сложность до O(N1/2 * log(N)), т. к. количество проверяемых чисел уменьшилось до O(N1/2). И несмотря на некоторые потери из-за усложнения алгоритма, на достаточно больших N время выполнения предсказуемо улучшилось на порядки.

    Сделав еще несколько небольших оптимизаций, автор улучшил время выполнения еще в 3 раза, и на этом неплохом результате решил остановиться.

    Еще не дочитав до конца, я подумал о том, что при той же вычислительной сложности O(N1/2 * log(N)) наверняка будет гораздо быстрее работать не со строками, а напрямую с числами. То есть сгенерировать последовательность палиндромов, используя не операции со строками, а арифметические действия.

    И это оказалось совсем не сложно, ведь последовательность чисел-палиндромов одинаковой длины напоминает арифметическую прогрессию с постоянным шагом, которую надо корректировать каждые BASE, BASE2, BASE3,… элементов.

    Например, для десятичных палиндромов шириной 5 основной шаг будет равн +100:

    10001, 10101, 10201, 10301, 10401, 10501, 10601, 10701, 10801, 10901, 11001


    Последний элемент последовательности не является палиндромом и требует коррекции +10, до 11011

    11011, 11111, 11211, 11311, 11411, 11511, 11611, 11711, 11811, 11911, 12011


    Последний элемент элемент последовательности опять требует коррекции +10, до 12021

    Так доходим то 99-го и 100-го элементов: 19991, 20091

    Необходимая коррекция для 100-го элемента +10-99 = -89, в результате, получаем 20002 и продолжаем до тех пор, пока не дойдем до 99999.

    Так как арифметическая генерация палиндромов стала очень быстрой (в среднем всего несколько команд ассемблера), я решил, что генерация палиндромов в одном основании + проверка на палиндром в другом будет работать медленнее, чем параллельная генерация палиндромов в обоих основаниях и сравнение.

    В итоге, получился следующий код на C++:

    Вспомогательная структура с необходимыми степенями основания:
    #undef POWER
    #define POWER(exp)  m_powers[exp]
    
    template<typename IntT>
    struct BaseData
    {
    	static const size_t MAX_WIDTH = sizeof(IntT) * CHAR_BIT;
    
    	BaseData(
    		size_t base,
    		IntT maxValue)
    		: m_base(base)
    	{
    		POWER(0) = IntT(1);
    		for (size_t i = 1; i < MAX_WIDTH; ++i) {
    			POWER(i) = POWER(i - 1) * IntT(base);
    			if (POWER(i) >= maxValue) {
    				m_maxWidth = i - 1;
    				break;
    			}
    		}
    	}
    
    	size_t  m_base;
    	size_t  m_maxWidth;
    	IntT    m_powers[MAX_WIDTH];
    };
    

    Итератор палиндромов:
    template<typename IntT>
    class Iterator
    {
    	static const size_t MAX_WIDTH = sizeof(IntT) * CHAR_BIT;
    
    private:
    	struct IncrementData
    	{
    		size_t  m_counter;         //	current counter value
    		size_t  m_counterLimit;    //	number of iterations, usually B, but B - 1 for last increment
    		IntT    m_increment;       //	increment value
    	};
    
    public:
    	inline Iterator(
    		const size_t base,
    		const IntT * powers)
    		: m_base(base)
    		, m_powers(powers)
    	{
    		m_value = POWER(0);
    
    		SetWidth(1);
    	}
    
    	inline void operator ++()
    	{
    		//	always increment by basic
    		m_value += m_basicIncrement;
    
    		size_t i = 0;
    		do {
    			if (--m_counters[i].m_counter)
    				return;
    
    			//	reset counter
    			m_counters[i].m_counter = m_counters[i].m_counterLimit;
    			//	correct value on digit overflow
    			m_value += m_counters[i].m_increment;
    		} while (++i < m_countersSize);
    
    		//	prepare to next width
    		SetWidth(m_width + 1);
    	}
    
    	inline const IntT & operator *() const
    	{
    		return m_value;
    	}
    
    private:
    	void SetWidth(
    		size_t width)
    	{
    		m_width = width;
    		m_countersSize = (m_width + 1) / 2;
    
    		m_basicIncrement = POWER(m_countersSize - 1);
    
    		size_t i;
    		for (i = 0; i < m_countersSize - 1; ++i) {
    			m_counters[i].m_counter = m_counters[i].m_counterLimit = m_base;
    			m_counters[i].m_increment = POWER(m_countersSize - i - 2) - POWER(m_countersSize - i - 2) * m_base * m_base;
    		}
    		m_counters[i].m_counter = m_counters[i].m_counterLimit = m_base - 1;
    		m_counters[i].m_increment = POWER(0) - POWER(1);
    
    		if (m_width & 1)
    			m_counters[0].m_increment += POWER(m_countersSize);
    		else
    			m_basicIncrement += POWER(m_countersSize);
    	}
    
    	IntT           m_value;                 //	current value
    	IntT           m_basicIncrement;        //	basic increment (100... for odd width, 1100... for even width
    	size_t         m_countersSize;          //	just greater half of width: (width + 1) / 2
    	IncrementData  m_counters[MAX_WIDTH];
    	size_t         m_width;                 //	current width = number of digits
    	size_t         m_base;
    	const IntT   * m_powers;
    };
    

    И, наконец main:
    int _tmain(int argc, _TCHAR* argv[])
    {
    	int64 limit = 18279440804497281;
    
    	BaseData<int64> base0(2, limit);
    	BaseData<int64> base1(10, limit);
    
    	Iterator<int64> it0(base0.m_base, base0.m_powers);
    	Iterator<int64> it1(base1.m_base, base1.m_powers);
    
    	while (*it0 <= limit) {
    		if (*it0 < *it1)
    			++it0;
    		else if (*it0 > *it1)
    			++it1;
    		else {
    			std::cout << *it0 << std::endl;
    			++it0;
    			++it1;
    		}
    	}
    
    	return 0;
    }
    

    56-й палиндром, равный 18279440804497281 был получен через 1.85 секунды, что примерно в 150 раз быстрее предыдущего результата (предполагая, что компьютер grey_wolfs имел вычислительную мощность, сходную с моим Intel Core i7-3770 CPU @ 3.40GHz). Безусловно, заметная часть моего преимущества была вызвана переходом с PHP на C++, но все равно я удовлетворенно потирал руки: алгоритм уже перебирал почти 300 000 000 палиндромов в секунду, а у меня в запасе было еще два козыря: генерировать только нечетные десятичные палиндромы (-20-25%) и использовать многопоточность (-70-85% при 8 потоках)…

    И тут я увидел комментарий, который меня просто «убил»:

    @hellman
    Не так давно на одном из контестов codechef была такая же задачка. Только базы систем исчисления произвольные от 2 до 16, и нужно было первые 1000 палиндромов меньшие 2^60. Ограничение по времени — 8 сек (правда на входе может быть 5 тестов). Editorial там же есть.

    Мой алгоритм находил все палиндромы до 260 за 15 секунд, т.е. в худшем случае не укладывался бы в ограничение по времени даже используя многопоточность. А ведь в Editorial был еще и издевательский комментарий «why does the time limit for this question is so high 8 sec I haven't seen such high limit in any questions...?».

    Удовлетворение быстро сменилось разочарованием: моя реализация перебора с вычислительной сложностью O(N1/2 * log(N)) хоть и была близка с минимальному теоретическому пределу, но этого все равно было недостаточно. Я попытался разобраться в теоретическом решении, описанном в Editorial, даже посмотрел на их код, но с первого раза понял только то, что они сумели решить задачу с вычислительной сложностью около O(N1/4 * log(N)).

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

    Продолжение следует.
    Поделиться публикацией
    Похожие публикации
    Реклама помогает поддерживать и развивать наши сервисы

    Подробнее
    Реклама
    Комментарии 3
    • 0
      Отличное продолжение!
      К сожалению моя реализация на PHP напрямую с числами не дала большого прироста. Всего раза в 3-5 вроде, уже точно не помню. Даже не стал её развивать и оптимизировать, ведь в задачке, которую упомянул hellman всё на порядки быстрее. И я с удовольствием почитаю разбор решения этой задачки.
      У меня самого, к сожалению, времени пока на это не хватает.
      • +1
        Исправил на «Безусловно, заметная часть моего преимущества была вызвана переходом с PHP на C++, но все равно я удовлетворенно потирал руки» :)

        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.