I adapted billtavis' code to collect bits. First I changed the getGradient function to return the continuous iteration count and the sign of the last iterate:

void getGradient(long double cx, long double cy, long double &xDer,

long double &yDer, long double &gzx, long double &gzy,

long double &length, unsigned int maxIter,

long double &m, bool &b) { // these two are the new ones

...

b = zy < 0 || (zy == 0 && zx < 0);

m = n + 1 - std::log2(std::log(length) / (0.5 * std::log(huge)));

...

}

I used a different formula for m, because I had issues with the original (but they only differ by a constant factor and offset afaict) - the one I use has integer spacing of dwell bands and m = 2 for |c| = sqrt(huge) = escape radius. I also wrote a wrapper without the extra arguments, that just calls the one with extra arguments, so I didn't have to change so much code in externalAngle().

Then I added bit collection in externalAngle():

typedef std::vector<bool> BITS;

BITS externalAngle(long double cx, long double cy, size_t maxIter, long double &e) { // returns a vector of bits, atan2 result in e

...

long double m = maxIter * 2;

long double m_old = m;

BITS bs;

bool b;

while (m_old > 2) { // this replaces the previous while loop start

getGradient(ckx,cky,xDer,yDer,gzxA,gzyA,length,maxIter, m, b);

if (std::ceil(m) < m_old) { // collect a bit if we crossed a dwell band

bs.push_back(b);

m_old = m;

}

...

dist = std::min(std::max(16.0L, 0.01 * std::sqrt(ckx * ckx + cky * cky)),0.5*dist); // replaces the previous clamping, 16 is a tiny step size with very large |c|

... // rest of loop body unchanged

}

e = (std::atan2(cky,ckx));

reverse(bs.begin(), bs.end()); // the bits were collected least-significant first

return bs;

}

Then there's the main program which takes command line arguments, converts the atan2 result from long double to a vector of bits, compares it with the (hopefully more accurate) vector of bits result, and prints them out. See the attached source code file for details.

Here are some tests, which seem to indicate the atan2 result is valid to about 16 bits, even though long double has 64 bits of mantissa precision:

$ # cre cim iter # pre(periodic) angle of nearby landing point

$ ./tavis 0.251 0 1000 # .(0)

.0000000000000000000000000000000000000000000000000000000000000000

.0000000000000000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000

$ ./tavis -2.001 0 1000 # .1(0)

.100000000000 0000000000000000000000000000000000000000000000000000

.100000000000

$ ./tavis -0.75 0.01 1000 # .(01)

.01010101010101001 01000001111000101111100010001010101110101011011

.01010101010101010 1010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101001001111

$ ./tavis 1e-16 1 1000 # .0(01)

.0010101010101001 111110001001100001110001111010000110100110001111

.0010101010101010 101010101010101010101010101100100111011

$ ./tavis -1.749 1e-10 1000 # .(011)

.01101101101100 00100000101001101000000110000110110110011010111011

.01101101101101 101101101101101101101101101101001001001001001001001001001001001100100100101001110001

$ ./tavis -1.01 0.251 1000 # .(01011001)

.010110010101011 1001111010111001000000101111001011001111100001000

.010110010101100 101011001010110010101100101011001010110010101100101011001010110100100111111

$ ./tavis -0.99 0.251 1000 # .(01010110)

.01010110010101011 00111001100111100010010001110011010011101111011

.01010110010101100 1010110010101100101011001010110010101100101011001010110010101010110100011101