Compute sqrt(SIZE_MAX+1) using only integer constant expressions, catering for weird ABIs

Question

OpenBSD's C library has an extension called reallocarray(3) which does realloc(array, size*nmemb) without blowing up if the multiplication overflows. The implementation contains this fragment:

/*
 * This is sqrt(SIZE_MAX+1), as s1*s2 <= SIZE_MAX
 * if both s1 < MUL_NO_OVERFLOW and s2 < MUL_NO_OVERFLOW
 */
#define MUL_NO_OVERFLOW (1UL << (sizeof(size_t) * 4))

Over on Programmers.SE a modified version of that calculation got dinged for technical incorrectness. 4 should obviously be CHAR_BIT/2, but that's not the only problem. Suppose an unusual ABI in which size_t has padding bits. (This is not ridiculously implausible: consider a microcontroller with 32-bit registers but a 24-bit address space.) Then SIZE_MAX is less than 1 << (sizeof(size_t)*CHAR_BIT) [in infinite-precision arithmetic] and the calculation is wrong.

So, question: Can you compute floor(sqrt(SIZE_MAX+1)) using only C99 integer-constant-expression arithmetic, making no assumptions whatsoever about the ABI other than what C99 requires plus what you can learn from <limits.h>? Note that SIZE_MAX may equal UINTMAX_MAX, i.e. there may not be any type that can represent SIZE_MAX+1 without overflow.

EDIT: I think SIZE_MAX is required to be 2ⁿ − 1 for some positive integer n, but is not necessarily of the form 2²ⁿ − 1 — consider S/390, one of whose ABIs has a 31-bit address space. Therefore: If sqrt(SIZE_MAX+1) is not an integer, the desired result (given how this constant is used) is floor() of the true value.

Answer 1

The constant SIZE_MAX is non-negative and has type size_t. For shortness, I will define:

  #define S SIZE_MAX  

The mathematical value S+1 is or can be, as you pointed out, out of range for any integer type.
I will write S1 for the mathemtaical value of S+1.
If we consider the logarithm (in base 2, if you want) of S1, then we have:

 logarithm(sqrt(S1)) == (1.0/2.0) logarithm(S1)

On the other hand, in almost-for-sure every realistic situation, we will have that S is represented as a binary number having only 1 bits. The number b of this bits is, in general, the number CHAR_BIT multiplied by a power of two, multiplied by CHAR_BIT: 16, 32, 64, 128... I will denote the exponent of this power by p. Thus, for CHAR_BIT == 8, we have:

16 == CHAR_BIT * 2 ----> p == 1
32 == CHAR_BIT * 4 ----> p == 2
64 == CHAR_BIT * 8 ----> p == 3

Now we have:

logarithm(S1) == b == CHAR_BIT * (2 ** p)   (I am denoting with ** to the "power math. operator").

logarithm(sqrt(S1)) == logaritm(S1) / 2.0 == CHAR_BIT * (2 ** p) / 2.0 == CHAR_BIT * (2 ** (p - 1))

By assuming or knowing that every bit in size_t is used only to represent the bits of an integer number, we have this equality, por some (unknown) value of p:

sizeof(size_t) == b == CHAR_BIT * (2 ** p)

We can assume, for the state-of-the-art in 2014, that the value of p <= 5, say (you can grow this magic number 5 to bigger values in what follows).

Now, consider the following expression, intended to "search and found" the value of b, under the assumption that p <= 5:

#define S_1 ((size_t)1ULL)
#define b (sizeof(size_t))
#define bitexpr(p) ((size_t)(CHAR_BIT * (S_1 << (p))))
#define expr(p) ((size_t) (S_1 << (p)))
#define exp2_expr_1(p) ((size_t)(S_1 << bitexpr(p-1)))
// SRSM() stands for: Square Root SizeMax
#define SRSM  ( \
    (expr(1)==b)? exp2_expr_1(1) :             \
        (expr(2)==b)? exp2_expr_1(2) :         \
            (expr(3)==b)? exp2_expr_1(3) :     \
              (expr(4)==b)? exp2_expr_1(4) :   \
                (expr(5)==b)? exp2_expr_1(5) : \
                   (size_t)0  /* Error! */     \
    ) /* end-of-macro*/

The macro SRSM brings, actually, the square root of S+1, but I suppose you can figure out what to do with this number.

What is important here is that the square root of SIZE_MAX can be obtained by using purely integer constant expressions.

If you want, the "magic" number 5 can be changed by another one.

A more general approach, intended to solve an arbitrary situation, on any possible machine fitting the standard, it would be more complicated.
The method used in this post is independent of the value that has CHAR_BIT, but it uses that the number of bytes is a power of 2.

EDITED: I changed a little the method for "search", starting by 1 and then growing up, to avoid possible "false" matches with << operator and big numbers (one never knows...). Now, the first match is, for sure, the correct.

Answer 2

An alternative to solve the sqrt(SIZE_MAX+1) is to look at the higher level goal shown here and below.

/* s1*s2 <= SIZE_MAX if s1 < K and s2 < K, where K = sqrt(SIZE_MAX+1) */
const size_t MUL_NO_OVERFLOW = ((size_t)1) << (sizeof(size_t) * 4);
if ((nmemb >= MUL_NO_OVERFLOW || size >= MUL_NO_OVERFLOW) &&
    nmemb > 0 && SIZE_MAX / nmemb < size)
  abort();

---

A simple overflow detection could use the following. This unfortunately performs division - an expansive test, yet otherwise "works". No MUL_NO_OVERFLOW needed.

if (nmemb > 0 && (SIZE_MAX / nmemb < size))
  abort();

Instead of dividing, OP's code does 2 simple compares against MUL_NO_OVERFLOW to eliminate most contenders.

if ((nmemb >= MUL_NO_OVERFLOW || size >= MUL_NO_OVERFLOW) && ...

Yet the pre-processor computation of MUL_NO_OVERFLOW is challenging , hence this post.

The alternative is to use an easy to compute MUL_NO_OVERFLOW_ALT to eliminate many size, nmemb pairs from the division test.

// good for up to SIZE_MAX + 1 == 2**128
#define SQN(x) ( !!(x) )
#define SQ0(x) ( (x) >= 0x2u ? 1 + SQN((x)/0x2u) : SQN(x) )
#define SQ1(x) ( (x) >= 0x4u ? 2 + SQ0((x)/0x4u) : SQ0(x) )
#define SQ2(x) ( (x) >= 0x10u ? 4 + SQ1((x)/0x10u) : SQ1(x) )
#define SQ3(x) ( (x) >= 0x100u ? 8 + SQ2((x)/0x100u) : SQ2(x) )
#define SQ4(x) ( (x) >= 0x10000u ? 16 + SQ3((x)/0x10000u) : SQ3(x) )
#define SQ5(x) ( (x) >= 0x100000000u ? 32 + SQ4((x)/0x100000000u) : SQ4(x) )
#define SQ6(x) ( (x)/0x100000000u >= 0x100000000u ? 64 + SQ5((x)/0x1000000000000000u/16) : SQ5(x) )
#define MUL_NO_OVERFLOW_ALT (((size_t)1) << (SQ6(SIZE_MAX)/2))

printf("%zx %zx\n", SIZE_MAX, MUL_NO_OVERFLOW_ALT);
// Output
// ffffffff 10000

---

In OP's case, even if the best (size_t)sqrt(SIZE_MAX + 1.0) was possible at compile time, the division test would still be needed to pass legitimate pairs of size, nmemb like 4, SIZE_MAX/100. So for this task, the best ISQRT(SIZE_MAX + 1.0) is not needed, just something near and not greater than.

With this code the best value is calculated when SIZE_MAX + 1 is an even power of 2 - a common case.

Answer 3

The following works but needs clean-up.

Assumes:
1) SIZE_MAX + 1 is some power of 2, odd or even.
2) unsigned long long range at least size_t.
3) sizeof (SIZE_MAX) <= 128 bits.

Notes:
1) SIZE_MAX must be at least 65535.
2) Way too many ()
3) Needs more testing.

Method:
Find power of 2 of SIZE_MAX/4 + 1 via nested defines. Needs another level for each double of bits. At each bit width level 32, 64, 128, ... test upper bits of SIZE_MAX. Once the maximum power-of-2 bit width is found, recurse down to get bit with of SIZE_MAX.

Form the sqrt(SIZE_MAX + 1) from (1 << Bit_Width(SIZE_MAX/4 + 1)/2)*2. Multiply by a scaled sqrt(2) if bit width was odd.

#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>

int main(void) {
  // Swap 1 for various powers of 2 like 256, 512. etc. to test
#define M  SIZE_MAX/1
  printf("M %zX\n", (size_t) M);
#define N (M/4 + 1)
  printf("N/4 + 1 %zX\n", (size_t) N);

// Obviously this line can be simplified.
#define M0(x)  ((x & ~0x0U) ? 0 : 0)

#define M1(x)  ((x & ~0x1U) ? 1 + M0(x >> 2) : M0(x))

#define M2(x)  ((x & ~0x3U) ? 2 + M1(x >> 2) : M1(x))

#define M3(x)  ((x & ~0xFU) ? 4 + M2(x >> 4) : M2(x))

#define M4(x)  ((x & ~0xFFU) ? 8 + M3(x >> 8) : M3(x))

#define M5(x)  ((x & ~0xFFFFLLU) ? 16 + M4(x >> 16) : M4(x))

#if (M & ~0xFFFFFFFFLLU)
#define M6(x)  ((x & ~0xFFFFFFFFLLU) ? 32 + M5(x >> 32) : M5(x))
#if (M & ~0xFFFFFFFFFFFFFFFFLLU)
#define MN(x)  ((x & ~0xFFFFFFFFFFFFFFFFLLU) ? 64 + M6(x >> 64) : M6(x))
#define MSQ2(odd, x) ((odd&1) ? (x*13043817825332782212llu)>>63 : x)
#else
#define MN(x)  M6(x)
#define MSQ2(odd, x) ((odd&1) ? (x*3037000499llu)>>31 : x)
#endif
#else
#define MN(x)  M5(x)
#define MSQ2(odd, x) ((odd&1) ? (x*46340llu)>>15 : x)
#endif

  printf("MN(N) %d\n", MN(N));
  printf("MN(M) %d\n", MN(M));

#define MUL_NO_OVERFLOW   ((size_t) MSQ2(MN(N), ((((size_t)1) << (MN(N)/2))*2) ))

  printf("MUL_NO_OVERFLOW %zu\n", MUL_NO_OVERFLOW);
  return 0;
}

---

Looks too complicated - will see about simplifying. GTG