Stable vectors

Extremely convenient as they are, std::vectors have a limitation that many novice C++ programmers frequently stumble upon: iterators and references to an element of an std:vector are invalidated when a preceding element is erased or when the vector expands and needs to migrate its internal storage to a wider memory region (i.e. when the required size exceeds the vector's capacity). We say then that std::vectors are unstable: by contrast, stable containers are those for which references and iterators to a given element remain valid as long as the element is not erased: examples of stable containers within the C++ standard library are std::list and all the associative containers.

Sometimes stability is too precious a feature to live without. One particular property of std::vectors, namely element contiguity, makes it impossible to add stability to this container. So, provided we sacrifice element contiguity, how much can a stable design approach the behavior of std::vector (random access iterators, amortized constant time end insertion/deletion, minimal memory overhead, etc.)?

The following image describes the layout of a possible data structure upon which to base the design of a stable vector:

Each element is stored in its own separate node. All the nodes are referenced from a contiguous array of pointers, but also every node contains an "up" pointer referring back to the associated array cell. This up pointer is the key element to implementing stability and random accessibility:

1. Iterators point to the nodes rather than to the pointer array. This ensures stability, as it is only the pointer array that needs to be shifted or relocated upon insertion or deletion.
2. Random access operations can be implemented by using the pointer array as a convenient intermediate zone. For instance, if the iterator it holds a node pointer it.p and we want to advance it by n positions, we simply do:

   ‎it.p = *(it.p->up+n);

   That is, we go "up" to the pointer array, add n there and then go "down" to the resulting node.
   

In a later entry we will present a full implementation of an STL-conformant stable vector and compare its characteristics (space/time complexity, exception safety) with those of std::vector.

Cpp bignum arithmetic with memoization

The multiplication algorithm we implemented for C preprocessor bignums incurs a certain amount of repetitive work: multiplying bn1 by bn2 = bn21 bn22···bn2m, involves calculating the set of intermediate products bn1·bn21, ... ,bn1·bn2m, which includes as many repeated values as equal digits are in the decimal expansion of bn2. We can avoid recalculating those by applying memoization to the BN_DIGIT_MUL macro:

‎#define BN_CACHED_DIGIT_MUL(bn,y,cache) \
‎BOOST_PP_IF( \
‎  BOOST_PP_EQUAL( \
‎    1, \
‎    BOOST_PP_SEQ_SIZE(BOOST_PP_SEQ_ELEM(y,cache))), \
‎  BN_CACHED_DIGIT_MUL_CASE1, \
‎  BN_CACHED_DIGIT_MUL_CASE2)(bn,y,cache)
‎
‎#define BN_CACHED_DIGIT_MUL_CASE1(bn,y,cache) \
‎BN_CACHED_DIGIT_MUL_CASE1_IMPL(BN_DIGIT_MUL(bn,y),y,cache)
‎
‎#define BN_CACHED_DIGIT_MUL_CASE1_IMPL(bny,y,cache) \
‎(bny)(BOOST_PP_SEQ_REPLACE(cache,y,(~)bny))
‎
‎#define BN_CACHED_DIGIT_MUL_CASE2(bn,y,cache) \
‎(BOOST_PP_SEQ_TAIL(BOOST_PP_SEQ_ELEM(y,cache)))(cache)

BN_CACHED_DIGIT_MUL(bn,y,cahe) accepts a cache memory in the form of a Boost.Preprocessor sequence storing those of the values bn·0, bn·1, ..., bn·9 calculated so far (with a special format not worth discussing in detail here). If the cache already has the requested multiplication, then this is used, otherwise we calculate the value using BN_DIGIT_MUL and update the cache accordingly. The return value is a pair containing both the result and the updated cache, which can be passed later to subsequent invocations to BN_CACHED_DIGIT_MUL. To make BN_MUL benefit from the usage of BN_CACHED_DIGIT_MUL we have to add the corresponding cache as an additional member of the state structure that constitutes the core of the multiplication algorithm.

‎‎#define BN_MUL(bn1,bn2) \
‎BN_MUL_STATE_RESULT( \
‎  BOOST_PP_SEQ_FOLD_LEFT( \
‎    BN_MUL_OP, \
‎    ((~))(bn2) \
‎    (((~)(0))((~)bn2)((~))((~))((~))((~))((~))((~))((~))((~))) \
‎    ((0)), \
‎    BOOST_PP_SEQ_REVERSE(bn1)))
‎
‎#define BN_MUL_OP(s,state,x) \
‎BN_MUL_OP_IMPL( \
‎  x, \
‎  BN_MUL_STATE_SHIFT(state), \
‎  BN_MUL_STATE_BN2(state), \
‎  BN_MUL_STATE_CACHE(state), \
‎  BN_MUL_STATE_RESULT(state))
‎
‎#define BN_MUL_OP_IMPL(x,shift,bn2,cache,result) \
‎BN_MUL_OP_IMPL_RETURN( \
‎  BN_CACHED_DIGIT_MUL(bn2,x,cache),shift,bn2,result)
‎
‎#define BN_MUL_OP_IMPL_RETURN(bn2x_cache,shift,bn2,result) \
‎(shift(0)) \
‎(bn2) \
‎(BOOST_PP_SEQ_ELEM(1,bn2x_cache)) \
‎(BN_ADD( \
‎   result, \
‎   BOOST_PP_SEQ_ELEM(0,bn2x_cache)BOOST_PP_SEQ_TAIL(shift)))
‎
‎#define BN_MUL_STATE_SHIFT(state) BOOST_PP_SEQ_ELEM(0,state)
‎#define BN_MUL_STATE_BN2(state) BOOST_PP_SEQ_ELEM(1,state)
‎#define BN_MUL_STATE_CACHE(state) BOOST_PP_SEQ_ELEM(2,state)
‎#define BN_MUL_STATE_RESULT(state) BOOST_PP_SEQ_ELEM(3,state)

The initial cache value is marked in red above, note that it already contains precomputed values for the easy cases bn2·0 and bn2·1.

We can also apply memoization to the implementation of BN_DIV: the macro BN_QDIGIT_REMAINDER(bn1,bn2), which returns the largest decimal digit q such that bn2·q ≤ bn1, along with bn1 − bn2·q, can be supplemented with a cache to store the values bn2·0, bn2·1, ..., bn2·9. The implementation details are entirely similar to those involved in memoizing BN_MUL. A full implementation of the memoized bignum arithmetic operations is provided.

Let us measure the improvement brought about by memoization. These 10-, 20- and 30-digit random numbers were selected with the aid of random.org:

n₁ = 8579726378,
n₂ = 88242457695260718048,
n₃ = 7278842148,
n₄ = 17478148223952531376,
n₅ = 22737602045089519821,
n₆ = 859621743146723715335628844223,
n₇ = 97249055978154178611,
n₈ = 502003375419877211491785910042,

and the following operations were timed for the plain and memoized versions of bignum multiplication and division macros; figures are averages of several individual measurements, tests were run with GCC 3.4.4 for Cygwin in a Windows box with a 1.66 GHz Intel Core 2 Duo and 1 GB RAM:

operation	plain version (s)	memoized version (s)	speedup factor
n1 · n2	11.2	9.1	19%
n2 · n1	13.3	9.7	27%
n3 · n4	9.8	7.0	28%
n4 · n3	12.0	9.7	19%
n5 · n6	27.2	20.6	24%
n6 · n5	30.0	22.3	25%
n7 · n8	27.9	19.9	28%
n8 · n7	31.7	22.0	30%
n2 / n1	8.9	7.8	13%
n4 / n3	5.9	4.9	17%
n6 / n5	11.7	9.5	19%
n8 / n7	11.2	9.8	12%

The speedup factor is defined as (p − m)/p, where p is the time taken by the plain version and m the time with memoization. Improvements are noticeable though not dramatic. Multiplication could be further accelerated by using alternative algorithms like for instance Karatsuba's, but we decide to leave it at this for the moment.

Escaping transcendence

When I was a teenager I was fond of 19th century novels such as Flaubert's Sentimental Education or Dostoevsky's Crime and Punishment. In music, my absolute favorite were Pink Floyd's long, discursive compositions: I remember myself those days listening to Comfortably Numb over and over again full of adolescent angst.

There was a uniting theme to these works: they were important, no pastime, bigger than life in some sense. The subjects proposed were deliberately profound: morality, redemption, isolation. There was a search for some sort of transcendence in the literature and music I visited in my youth days.

Twenty years later, I found myself reading Pynchon's Gravity's Rainbow and appreciating the beauty of Portishead's new album, Third. The contrast is stark: some profundity remains, but in a subtler form. Discourse has been supplanted by evocation, and layers of reference take the place of the distinctive stories of the past. Whereas the pieces of my youth sought to prove some thesis or convey important information to a more or less passive audience, Rainbow or Third draw heavily from the reader/listener mental background and cannot exist in isolation from the audience and their implied cultural setting.

I think this is a recurring pattern in the way art has evolved in our days: we have moved from transcendence to self-reference. Art today is less about the world and much more about art itself and its relationship with the viewer. If artistic works of the past projected the viewer into a bigger, external reality, current art constitutes a mere tool for exquisite masturbation.

Indirect survey

The usual way to statistically determine the prevalence of a given trait T in a population involves selecting a representative population sample and having their members answer the question "Do you have trait T?"; the number of affirmative questions divided by the size of the sample is an estimate of the actual prevalence of p = P(T). How can we estimate p if the following question is used instead?

Do you have some friend with trait T?

Let us assume that having a friend with T is independent of the fact that other friends may or may not have T. The probability then that a person has a friend with T is given by:

p' = 1 − ∑f_i(1 − p)ⁱ,
f_i := P(a person has exactly i friends).

So, p' is simply 1 − G_F(1 − p), where G_F is the probability-generating function of F, the random variable expressing the number of friends a person has. If we model F by a Poisson distribution with mean λ, we have

G_F(x)= e^{λ(x − 1)}

and

p' = 1 − e^−λp,
p = −(1/λ) ln(1 − p').

The figure shows p' as a function of p for various values of λ.

There are situations where it can be interesting to estimate p through this indirect survey method:

1. For usual values of λ much greater than 1 , p' is notably larger than p when the latter lies in the vicinity of 0. This implies that indirect estimation of p through p' results in better (narrower) confidence intervals when p is very small. For instance, for p = 1/1000 and λ = 10 a direct survey obtains on average one affirmative answer in a sample of 1000 individuals, with an associated 95% Wilson score confidence interval of (0.0001, 0.0065) (calculated via VassarStats). With an indirect survey we typically get 10 affirmative answers for the same size of population, with a confidence interval (0.0051, 0.019) associated to p'; applying the transformation p = −(1/λ)ln(1 − p'), the latter confidence interval maps to (0.0005, 0.0019), whose width is ~20% of that of the interval for the direct survey.
2. When doing surveys on sensitive issues (illegal drug use, paid sex usage) results can be severely biased due to the reluctance of people to answer thruthfully. In these cases, survey respondents might be more inclined to answer accurately when the questions are indirect: it is easier for a person to recognize that he has friends who had paid for sex than admitting paying for sex oneself.

A serious drawback of this method is its reliance on the a priori unkown parameter λ. One technique for solving this is the following: include into the survey an additional control question for which the associated statistics are already known, so that the survey results can be used to estimate λ. Ideally, the control question's subject matter should be similar to that of the actual question so as to cancel evasive answer bias out.

Formalizing predicate boundedness

Informally, we defined a predicate F to be bounded if the entities for which F is true can be enclosed in a reasonably "small" class. Of course, for this definition to be of any use we have to clearly specify what we mean by "small", or at least use our intuition to propose some axioms for this concept:

• If A and B are small, A U B is small.
• If A and B are small, A ∩ B is small.
• If A is small and B is included in A, B is small.
• If A is small then −A, the class of all entities outside A, is not small.
  

In the definition of −A given above we have assumed that we are working inside some universal class comprising any conceivable entity.

The axioms for "small" are far from categorical. We give a few models for this notion:

• If we decide to work within a mathematical class theory such as NBG, there is an easy definition for "small": A class is small if it is a set (not all classes are sets).
  
• We can take "small" as meaning "finite", provided the universal class is infinite, or as having cardinality less than or equal to some fixed infinite cardinal κ, provided the universal class is larger than κ.
• Suppose we have a fixed set of "primitive" predicates F_i with domains D_i such that UD_i is smaller than the universal class. We can define the class of small sets as the least class S comprising all D_i that is closed under finite union, finite intersection and comprehension: if X belongs to S and Y is a subset of X, Y belongs to S. This definition of "small" can be seen as a formal model for Quine's notion of natural kinds.
  

This is probably as far as we can get with respect to defining what "small" means in the context of predicate boundedness.

Cpp bignum arithmetic: part II

After having studied addition, subtraction and multiplication of C preprocessor bignums, let us now turn to division. To gain some insight into the mechanics of the division algorithm we reproduce a division layout the way we were taught to do at school:

‎         20553
‎8879265 |432
‎ 239
‎ 2392
‎  2326
‎   1665
‎    369

To divide 8879265 by 432, we begin by isolating the first 3-digit group of the dividend, 887, and calculate how many times the divisor 432 makes it into it; we have that 887 = 432·2 + 23, so 2 is the first digit of our result, and we write down the remainder 23. Next, we append the next digit of the dividend, 9 , to the remainder and proceed iteratively:

• 887 = 432·2 + 23,
  
• 239 = 432·0 + 239,
• 2392 = 432·5 + 232,
• 2326 = 432·5 + 166,
• 1665 = 432·3 + 369
  

Let us assume we have a macro BN_QDIGIT_REMAINDER(bn1,bn2) that returns a sequence (q)(r), where q is the greatest digit between 0 and 9 such that bn2·q ≤ bn1 and r is the bignum representing the value bn1 − bn2·q. We can then implement the division of bn1 = (bn11)···(b1n) by bn2 (bn21)···(b2m)by folding over reverse((bn1m)···(b1n)) with the following recursive state definition:

state_i = (lead_i,bn2,res_i);
state₀ = ((bn11)···(b1m-1), bn2,empty sequence),
lead_i+1 = first(BN_QDIGIT_REMAINDER(lead_i,bn2)),
res_i+1 = append(res_i, second(BN_QDIGIT_REMAINDER(lead_i,bn2))),

where the functions first and second simply return the first and second component of a sequence of size 2. This is the complete Boost.Preprocessor-based algorithm:

‎#define BN_DIV(bn1,bn2) \
‎BOOST_PP_IF( \
‎  BN_LESS(bn1,bn2), \
‎  BN_DIV_CASE1, \
‎  BN_DIV_CASE2)(bn1,bn2)
‎
‎#define BN_DIV_CASE1(bn1,bn2) (0)
‎
‎#define BN_DIV_CASE2(bn1,bn2) \
‎BN_DIV_POSTPROCESS( \
‎  BOOST_PP_SEQ_TAIL(BN_DIV_STATE_RESULT(\
‎    BOOST_PP_SEQ_FOLD_LEFT( \
‎      BN_DIV_OP, \
‎      ( \
‎        (0) \
‎        BOOST_PP_SEQ_FIRST_N( \
‎         BOOST_PP_DEC(BOOST_PP_SEQ_SIZE(bn2)),bn1)) \
‎      (bn2) \
‎      ((~)), \
‎      BOOST_PP_SEQ_REST_N( \
‎        BOOST_PP_DEC(BOOST_PP_SEQ_SIZE(bn2)),bn1)))))
‎
‎#define BN_DIV_POSTPROCESS(result) \
‎BOOST_PP_IF( \
‎  BOOST_PP_EQUAL(0,BOOST_PP_SEQ_ELEM(0,result)), \
‎  BOOST_PP_SEQ_TAIL, \
‎  BN_DIV_POSTPROCESS_CASE2)(result)
‎
‎#define BN_DIV_POSTPROCESS_CASE2(result) result
‎
‎#define BN_DIV_OP(s,state,x) \
‎BN_DIV_OP_IMPL( \
‎  BN_QDIGIT_REMAINDER( \
‎    BN_DIV_STATE_LEAD(state)(x), \
‎    BN_DIV_STATE_BN2(state)), \
‎  BN_DIV_STATE_BN2(state), \
‎  BN_DIV_STATE_RESULT(state))
‎
‎#define BN_DIV_OP_IMPL(qr,bn2,result) \
‎(BOOST_PP_SEQ_ELEM(1,qr)) \
‎(bn2) \
‎(result(BOOST_PP_SEQ_ELEM(0,qr)))
‎
‎#define BN_DIV_STATE_LEAD(state) BOOST_PP_SEQ_ELEM(0,state)
‎#define BN_DIV_STATE_BN2(state) BOOST_PP_SEQ_ELEM(1,state)
‎#define BN_DIV_STATE_RESULT(state) BOOST_PP_SEQ_ELEM(2,state)

We are then left with the task of implementing the macro BN_QDIGIT_REMAINDER(bn1,bn2). A brute-force approach is to try the digits 9,8,...,0 until we met the first one satisfying bn2·q ≤ bn1, but this is unnecesarily wasteful: we find in Brinch Hansen's paper "Multiple-length division revisited" the following nice theorem: if bn1 and bn2 are represented in such a way that length(bn1) = length(bn2) +1 (even if that means bn1 has to be padded with some some leading zeros) and we define

qe =min(9, floor((bn11·100 + bn12·10 + bn13)/(bn21·10 + bn22))),

then

q ≤ qe ≤ q + 1,

that is, we can use the first three digits of bn1 and the first two digits of bn2 to estimate q and will be off at most by one:

‎#define BN_QDIGIT_REMAINDER(bn1,bn2) \
‎BN_QDIGIT_REMAINDER_IMPL( \
‎  BN_ZERO_SEQ( \
‎    BOOST_PP_IF( \
‎      BOOST_PP_LESS_EQUAL( \
‎        BOOST_PP_SEQ_SIZE(bn1), \
‎        BOOST_PP_SEQ_SIZE(bn2)), \
‎      BOOST_PP_INC( \
‎        BOOST_PP_SUB( \
‎          BOOST_PP_SEQ_SIZE(bn2), \
‎          BOOST_PP_SEQ_SIZE(bn1))), \
‎      0))bn1, \
‎  bn2)
‎
‎#define BN_ZERO_SEQ(n) \
‎BOOST_PP_REPEAT(n,BN_ZERO_SEQ_OP,~)
‎
‎#define BN_ZERO_SEQ_OP(z,n,data) (0)
‎
‎#define BN_QDIGIT_REMAINDER_IMPL(bn1,bn2) \
‎BOOST_PP_IF( \
‎  BOOST_PP_LESS(BOOST_PP_SEQ_SIZE(bn2),2), \
‎  BN_QDIGIT_REMAINDER_CASE1, \
‎  BN_QDIGIT_REMAINDER_CASE2)(bn1,bn2)
‎
‎#define BN_QDIGIT_REMAINDER_CASE1(bn1,bn2) \
‎BN_QDIGIT_REMAINDER_CASE1_RETURN( \
‎  BOOST_PP_ADD( \
‎    BOOST_PP_MUL(BOOST_PP_SEQ_ELEM(0,bn1),10), \
‎    BOOST_PP_SEQ_ELEM(1,bn1)), \
‎  BOOST_PP_SEQ_ELEM(0,bn2))
‎
‎#define BN_QDIGIT_REMAINDER_CASE1_RETURN(n1,n2) \
‎(BOOST_PP_DIV(n1,n2))((BOOST_PP_MOD(n1,n2)))
‎
‎#define BN_QDIGIT_REMAINDER_CASE2(bn1,bn2) \
‎BN_QDIGIT_REMAINDER_CASE2_TEST( \
‎  bn1,bn2, \
‎  BN_DIV_RRR_DD( \
‎    BOOST_PP_ADD( \
‎      BOOST_PP_MUL(BOOST_PP_SEQ_ELEM(0,bn1),10), \
‎      BOOST_PP_SEQ_ELEM(1,bn1)), \
‎    BOOST_PP_SEQ_ELEM(2,bn1), \
‎    BOOST_PP_SEQ_ELEM(0,bn2), \
‎    BOOST_PP_SEQ_ELEM(1,bn2)))
‎
‎#define BN_QDIGIT_REMAINDER_CASE2_TEST(bn1,bn2,q) \
‎BN_QDIGIT_REMAINDER_CASE2_TEST_IMPL(bn1,bn2,q,BN_DIGIT_MUL(bn2,q))
‎
‎#define BN_QDIGIT_REMAINDER_CASE2_TEST_IMPL(bn1,bn2,q,bn2q) \
‎BOOST_PP_IF( \
‎  BN_LESS_WITH_ZERO(bn1,bn2q), \
‎  BN_QDIGIT_REMAINDER_CASE2_TEST_IMPL_CASE1, \
‎  BN_QDIGIT_REMAINDER_CASE2_TEST_IMPL_CASE2)(bn1,bn2,q,bn2q)
‎
‎#define BN_QDIGIT_REMAINDER_CASE2_TEST_IMPL_CASE1(bn1,bn2,q,bn2q) \
‎(BOOST_PP_DEC(q))(BN_SUB(bn1,BN_SUB(bn2q,bn2)))
‎
‎#define BN_QDIGIT_REMAINDER_CASE2_TEST_IMPL_CASE2(bn1,bn2,q,bn2q) \
‎(q)(BN_SUB(bn1,bn2q))
‎
‎#define BN_LESS_WITH_ZERO(bn1,bn2) \
‎BOOST_PP_IF( \
‎  BOOST_PP_EQUAL(BOOST_PP_SEQ_ELEM(0,bn1),0), \
‎  BN_LESS_WITH_ZERO_CASE1, \
‎  BN_LESS)(bn1,bn2)
‎
‎#define BN_LESS_WITH_ZERO_CASE1(bn1,bn2) \
‎BN_LESS(BOOST_PP_SEQ_TAIL(bn1),bn2)

The code above invokes the macro BN_DIV_RRR_DD(r1r2,r3,d1,d2) whose behavior is defined by:

BN_DIV_RRR_DD(r1r2,r3,d1,d2) = floor((r1r2·10 + r3)/(d1·10 + d2)).

The operation cannot be done directly because Boost.Preprocessor-based arithmetics are limited to numbers between 0 and BOOST_PP_LIMIT_MAG = 256. Instead, floor(r1r2/d1) is used as an initial upper-bound estimation which is then adjusted to the correct value using arithmetic operations that never exceed the BOOST_PP_LIMIT_MAG limit. The implementation of BN_DIV_RRR_DD is somewhat cumbersome but conceptually simple.

The complete implementation of all four arithmetic operations --addition, subtraction, multiplication and division-- is provided. Cpp bignum arithmetic is not very fast; for instance, given the following calculation of 30!:

‎BN_TO_LITERAL(
‎  BN_MUL((3)(0),BN_MUL((2)(9),BN_MUL((2)(8),BN_MUL((2)(7),
‎  BN_MUL((2)(6),BN_MUL((2)(5),BN_MUL((2)(4),BN_MUL((2)(3),
‎  BN_MUL((2)(2),BN_MUL((2)(1),BN_MUL((2)(0),BN_MUL((1)(9),
‎  BN_MUL((1)(8),BN_MUL((1)(7),BN_MUL((1)(6),BN_MUL((1)(5),
‎  BN_MUL((1)(4),BN_MUL((1)(3),BN_MUL((1)(2),BN_MUL((1)(1),
‎  BN_MUL((1)(0),BN_MUL((9)   ,BN_MUL((8)   ,BN_MUL((7),
‎  BN_MUL((6)   ,BN_MUL((5)   ,BN_MUL((4)   ,BN_MUL((3),
‎  (2))))))))))))))))))))))))))))))

it takes around 20 seconds for this expression to expand to

265252859812191058636308480000000

using GCC 3.4.4 for Cygwin in a Windows box with a 1.66 GHz Intel Core 2 Duo and 1 GB RAM. Within the same environment, the long division

‎BN_TO_LITERAL(
‎  BN_DIV(
‎    (2)(6)(5)(2)(5)(2)(8)(5)(9)(8)(1)(2)(1)(9)(1)(0)(5)
‎    (8)(6)(3)(6)(3)(0)(8)(4)(8)(0)(0)(0)(0)(0)(0)(0),
‎    (1)(8)(3)(7)(0)(8)(0)(0)))

takes around 12 seconds. In a later entry, we will explore ways to improve the performance of bignum multiplication and division through memoization techniques.

1015 + 6 = 1000000000000005.9

In a prior entry we saw how the internal representation of numbers in Google Docs spreadsheets as double precision IEEE 754 floating point quantities introduces rounding errors when the stored values get larger than 2⁵³. Yet, this does not explain the fact that 19! is displayed as 121645100408831984, since 19!=121645100408832000 can be represented exactly in double precision IEEE 754.

Let us see how numbers are rounded in the vicinity of 121645100408832000:

Each number on the left column is the least value that gets rounded off to the associated number on the right column: for instance, the numbers 1216···1945 ≤ n ≤ 1216···1959 get mapped to 1216···1952, and so on. Except for the rows in bold, each rounding interval is comprised of 17 or 15 numbers (instead of each interval being of length 16, due to the application of unbiased rounding) around its central value. This is entirely reasonable because 19! is 57 bits wide and double precision IEEE 754 can only represent up to 53 significant bits, hence the hops of size 2⁵⁷⁻⁵³ = 16. On the other hand, we have an anomaly between rows 4 and 5: the values from 1216···1992 to 1216···2008, which should map to 1216···2000, are instead rounded off to the preceding value 1216···1984.

The binary representation of 19! has 57 bits: if we divide by 16 we get the 53 bit wide number 7602818775552000, which ought be exactly representable in double precision IEEE 754. Pasting this into a Google Docs spreadsheet gets us:

which is off by one! If we keep dropping least significant bits from this number (i.e. dividing by 2) we obtain this:

Only the last number is shown correctly, but after that the display is correct all the way down to zero. This suggests that there is a lower limit below which we do not have displaying problems. After probing the number space more or less thoroughly, our findings are the following:

The problems occur from 10¹⁵ up. Displaying errors are extremely frequent and distributed very evenly: if we consider the sequence of exactly representable integers greater than 10¹⁵, roughly 11% of them are displayed incorrectly, with a typical distance between consecutive errors of 9, ocasionally a little more. The magnitud of the error is ~n/2⁵³. The smallest number displayed incorrectly is 10¹⁵ + 6 = 1000000000000006, which Google Docs spreadsheets show as 1000000000000005.9.

The problem seems to be related only to the displaying algorithm; numbers are stored correctly, as suggested by the following: