Next, in search of palindromes 2

Not so long ago I read the article grey_wolfs “Forward, in search of palindromes” on Habré about solving and optimizing a curious competitive task with a very laconic wording:

“The decimal number 585 is 1001001001 in binary. It is palindromic in both bases. Find n-th palindromic number ". Or, in Russian: “The decimal number 585 in binary number system looks like 1001001001. It is a palindrome in both number systems. Find the nth palindrome like. ”

The author began his fascinating story from the most trivial solution by searching and checking all the numbers from 1 to N and the computational complexity O (N * log (N)) , where N is the maximum checked number. The log (N) factor is necessary, since for each checked number, several actions are performed, depending on the number of its digits.

The first working optimization, namely, replacing the search of all numbers to search only decimal palindromes, dramatically improved the computational complexity to O (N ^1/2 * log (N)) , since the number of tested numbers decreased to O (N ^1/2 ) . And despite some losses due to the complexity of the algorithm, at a sufficiently large N , the execution time predictably improved by orders of magnitude.
')
Having made a few small optimizations, the author improved the execution time by another 3 times, and decided to stop at this good result.

Before I finish reading it, I thought that, with the same computational complexity, O (N ^1/2 * log (N)) would probably work much faster not with strings, but directly with numbers. That is, generate a sequence of palindromes, using not arithmetic operations, but arithmetic operations.

And it turned out to be not difficult at all, because the sequence of numbers of palindromes of the same length resembles an arithmetic progression with a constant step, which must be adjusted every BASE , BASE ² , BASE ³ , ... elements.

For example, for decimal palindromes with a width of 5, the main step would be +100 :

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


The last element of the sequence is not a palindrome and requires a correction of +10 , up to 11011

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


The last element of the sequence element again requires correction +10 , up to 12021

So we reach the 99th and 100th elements: 19991, 20091

The necessary correction for the 100th element + 10-99 = -89 , as a result, we get 20002 and continue until we reach 99999 .

Since the arithmetic generation of palindromes became very fast (only a few assembler commands on average), I decided that the generation of palindromes in one base + checking for palindromes in the other would work more slowly than parallel generation of palindromes in both bases and comparison.

The result is the following C ++ code:

Auxiliary structure with the necessary degrees of the base:

 #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]; }; 

Palindrome iterator:

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

And finally the 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; } 

The 56th palindrome, equal to 18279440804497281, was obtained in 1.85 seconds, which is about 150 times faster than the previous result (assuming that the gray-Wolf computer had a processing power similar to my Intel Core i7-3770 CPU @ 3.40GHz). Of course, a significant part of my advantage was caused by the transition from PHP to C ++, but still I rubbed my hands with satisfaction: the algorithm went through almost 300,000,000 palindromes per second, and I had two more trumps in stock: generate only odd decimal palindromes ( - 20-25 %) and use multithreading ( -70-85 % with 8 threads) ...

And then I saw a comment that simply “killed” me:

@hellman
Not so long ago in one of the codechef contests there was the same problem . Only the bases of the calculus systems are arbitrary from 2 to 16, and the first 1000 palindromes needed less than 2 ^ 60. The time limit is 8 seconds (though there may be 5 tests at the entrance). Editorial is there.

My algorithm found all palindromes up to 2 ⁶⁰ in 15 seconds, i.e. in the worst case, it would not fit into the time limit even using multithreading. But in the editorial there was also a mocking comment “why I’ve seen it?”.

Satisfaction was quickly replaced by disappointment: my implementation of a brute force with computational complexity O (N ^1/2 * log (N)), although it was close to the minimum theoretical limit, was still not enough. I tried to understand the theoretical solution described in the Editorial, I even looked at their code, but from the first time I understood only that they were able to solve the problem with computational complexity around O (N ^1/4 * log (N)) .

At that moment I stopped working on the puzzle, but it did not go out of my head. And a few days later I returned to her.

To be continued .