# Zero-initializing large dynamically allocated arrays of double

Arrays of double can be zero-initialized with the { 0 } initializer, which works even on systems that have a binary representation of 0.0 that is not all bits zero. While the IEEE 754 representation of floating point zero is all bits zero, the C Standard makes no guarantees about the representation of floating point numbers.

Large dynamically allocated arrays of integer types can be zero-initialized by using calloc(), but this method is not portable for floating point types, since calloc() initializes all bits to zero. Assigning a value of 0.0 to each individual double in the segment using a loop seems grossly inefficient.

I wrote the code for the function calloc_d() to solve this problem. The function takes the size_t argument nmemb, and returns a zero-initialized segment of memory large enough to hold nmemb doubles. If there is an allocation error, the function returns NULL. The caller is of course responsible for deallocation.

After using malloc() to allocate the required amount of memory, the first bytes are assigned the value of 0.0. Then memmove() is used to copy the bytes of the first zero, then the bytes of the first two zeros, and so on until at least half of the memory has been initialized. Finally, the remaining memory is initialized to 0.0 by copying bytes from the first half.

The code works, and has no obvious issues, as far as I can tell. I am interested in any comments about shortcomings of this approach, or shortcomings of my code, and suggestions for improvement. I would also be interested in suggestions for alternative, possibly more efficient methods.

I have included the function in a working program below. MAXCOUNT doubles are allocated for, initialized, and displayed. The value of MAXCOUNT is currently set to 100, but I have tested it for values up to 1000000.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAXCOUNT 100      // number of doubles to allocate for

double * calloc_d(size_t nmemb);

int main(void)
{
double *dbl_store = calloc_d(MAXCOUNT);
size_t i;

for(i = 0; i < MAXCOUNT; i++) {
if (i % 10 == 0) putchar('\n');
printf("%-8.4f", dbl_store[i]);
}
putchar('\n');

free(dbl_store);

return 0;
}

double * calloc_d(size_t nmemb)
{
double *ret;
double *next;         // pointer to beginning of uninitialized segment
size_t alloc_bytes = sizeof(*ret) * nmemb;
size_t init_sz;       // size of initialized segment

ret = malloc(alloc_bytes);
next = ret;

if (ret) {
ret[0] = 0.0;
init_sz = sizeof(*ret);
++next;

while (init_sz < (alloc_bytes + sizeof(*ret)) / 2) {
memmove(next, ret, init_sz);
init_sz *= 2;
next = ret + init_sz / sizeof(*ret);
}
memmove(next, ret, alloc_bytes - init_sz);
}

return ret;
}

• How long does a naive solution take? Is it actually making a significant difference compared to the runtime of the rest of your code? These questions should be answered before engaging in a serious attempt to optimise this. – Jack Aidley Dec 19 '16 at 9:45
• I strongly suspect that the compiler will optimise the naive loop to an assembly rep stosq (or its equivalent) on most platforms. "Repeat store string <size>" is a highly optimised in the hardware, and I doubt it can be beaten by anything else for this purpose. It will also never be slower than rep movsq(which I believe is usually how memcpy is implemented) - in fact stosq is almost always faster than movsq because it eliminates the need to read memory while writing to it. – Gavin Lock Dec 19 '16 at 11:06
• Don't use memmove() when you can be confident that memcpy() will suffice. That is indeed the case here. – PellMel Dec 19 '16 at 14:47
• Of course, I had to write the code first before I could test it against the naive solution, and it seemed like an interesting, if unnecessary in the end, problem to solve. But I was surprised to find my solution defeated so resoundingly by compiler optimizations. It looks like, consistently, my code is 18% faster than the naive solution given by @EmilyL with default optimizations, and the naive solution is 38% faster than mine with -O3 optimizations (on gcc). I learned something, which is why I come here :) – ex nihilo Dec 19 '16 at 14:47
• @PellMel-- memcpy() can't be used to copy overlapping arrays, and I was interpreting that to mean that it could not be used to copy from an allocated region of memory to the same region. But I now see that this was wrong, and that specifically memcpy() can't copy overlapping objects, which is not the case here. Thanks for pointing this out. – ex nihilo Dec 19 '16 at 14:57

Okay this guy is making a complicated solution to a simple problem based on a hunch that the trivial solution is "grossly inefficient" without any measurements to back up said claim. I suspect OPs solution is slower than the trivial solution.

So I set out test this, I pit your fancy implementation against my trivial implementation:

double * calloc_d2(size_t nmemb){
double* ret = malloc(sizeof(double)*nmemb);
double* end = ret + nmemb;

if(ret){
double* i = ret;
while(i != end){
*i = 0.0;
++i;
}
}
return ret;
}


I tested this here on ideone.com (I think it is using clang) which shows that this simple bit of code is consistently around 3-5% faster than your fancy code.

## So I have two pieces of advice

1. Express your intent as clearly as possible, and let the compiler worry about generating high performance code. Compilers are pretty smart...
2. Always, always, always measure.

Seeing as this answer has gotten quite some attention I've gone ahead and dug deeper.

Using godbolt compiler explorer I compiled (-O3) and disassembled my trivial loop implementation.

Here is the disassembly on GCC 6.2:

calloc_d2(unsigned long):
push    rbp
push    rbx
lea     rbx, [0+rdi*8]
sub     rsp, 8
mov     rdi, rbx
call    malloc
test    rax, rax
mov     rbp, rax
je      .L1
lea     rax, [rax+rbx]
cmp     rbp, rax
je      .L1
mov     rdx, rbx
xor     esi, esi
mov     rdi, rbp
call    memset
.L1:
mov     rax, rbp
pop     rbx
pop     rbp
ret


Notice how the compiler has just gone "Oh, you want to zero that? Cool, I'll just replace your code with a memset". This goes to show that an optimising compiler is quite smart and is able to transform your code in astounding ways.

I also tried different versions of GCC and clang on ARM and MIPS and PowerPC, they all generate similar instructions which essentially boil down to call malloc then memset.

With the sole exception of ICC 17 (-O3 -march=native) which goes above and beyond:

calloc_d2(unsigned long):
push      rbp                                           #5.1
mov       rbp, rsp                                      #5.1
and       rsp, -32                                      #5.1
push      r15                                           #5.1
push      rbx                                           #5.1
sub       rsp, 16                                       #5.1
mov       rbx, rdi                                      #5.1
lea       rdi, QWORD PTR [rbx*8]                        #6.50
call      malloc
-- snip --
call      __intel_avx_rep_memset                        #12.8
-- snip --


Notice that it calls into __intel_avx_rep_memset instead, this is a 512 bit wide AVX optimized memset that will only function on modern CPUs (it did this because I passed -march=native asking it to throw backward compatibility out the window). There is some branching code at the end which deals with cases where nmemb isn't evenly divisible by 512/64=8.

But here is the key take away, we clearly expressed our intent with the trivial loop. The compiler could then replaced the trivial loop (after deducing it's functionality) with an (512 bit) AVX optimised memset written by Intel's engineers specifically for modern Intel CPUs.

There is this thing called the "as-if rule" in C/C++. Which pretty much says:

As long as the program does the same reads and writes to volatile memory in the same order as-if the program was executed according to the language specification, then the program may do ANY transformation of the code it pleases. Including reordering instructions, reordering stores and loads to non-volatile memory, or removing instructions all together.

Just to show you how awesome some compilers can be look at what clang does here:

auto stupid(){
auto x = new int[30];
for(int i = 0; i < 30; i++)
x[i] = i;

auto sum = 0;
for(int i = 0; i < 30; i++)
sum += x[i]*x[i];
delete [] x;
return sum;
}


gets compiled into:

stupid():                             # @stupid()
mov     eax, 8555
ret


It deduces that the function always returns the same value and it completely removes all code and just returns a constant! It also removes the calls to new and delete as it realises they are pointless and don't affect the output of the function. If that is not impressive, I don't know what is. :)

So I would like to take this opportunity to reiterate my first piece of advice again:

Clearly express what it is you want your code to do. Let the compiler worry about generating good machine code, the compiler is likely much better at it than you are.

• First, thanks for testing this. I did not claim that the trivial solution is inefficient, but that it seems inefficient. I find it a little difficult to believe that the trivial loop-assignment version is faster than my solution, at least for large values of nmemb, which is the case this was designed for. What value did you pass into calloc_d() for the tests? It is certainly possible, and I haven't done any benchmarking yet, but I will look into this further. – ex nihilo Dec 19 '16 at 1:31
• Maybe the loop does get optimized to something like this, but isn't conventional wisdom that it is faster to use memcpy() to initialize some memory than to directly assign values? – ex nihilo Dec 19 '16 at 1:33
• @DavidBowling It doesn't make sense that memcpy() would be fast to initialize memory because memcpy() requires both reads and writes, and the reads are just wasted cycles. What would be fastest would be an optimized "memsetDouble" function that took a double sized bit pattern and set a block of memory to that pattern. The code that Emily L wrote probably is close to optimal if you use an optimizing compiler. But the fastest way would involve the same kinds of target specific optimizations that are used for memset() (which is possibly hand coded). – JS1 Dec 19 '16 at 3:14
• @JS1-- thanks. That makes sense. EmilyL-- I think that the issue was with my timing code, and that once optimized, your code is just flat out faster. – ex nihilo Dec 19 '16 at 4:30
• @EmilyL.-- Thanks for putting so much energy into answering my question. I have learned a lot from the answers and comments on this page, and have a new respect for compiler optimizations. FWIW, the optimizations don't seem to improve my function substantially, and it appears that for your version most of the improvements from optimization are already done in gcc by -01. I have plenty to think about over Christmas :) – ex nihilo Dec 20 '16 at 1:26

In the 99.9% of cases, your architecture will use IEEE floats/doubles, and using calloc() is ideal for both simplicity and speed. calloc() is far more efficent when allocating large zeroed memory since the operating system will often simply allocate a blank zero page of memory from the page table and not need to actually set each word's value.

So you can test when this is safe to do by simply noticing if you're using IEEE floats (std::numeric_limits::is_iec55 will be TRUE) and if so, use calloc(). Otherwise fall back to Emily's method.

• You can even check dynamically whether 0.0 has all bits zero, which may be true even for implementations that don't use IEEE-754 floating point representation. – PellMel Dec 19 '16 at 14:46

This exercise does appear to be trying to optimize a niche part of code and suspect focus on speed improvements are more useful elsewhere. Yet given that, here are some ideas:

1. Using memcpy()/memove() versus a simple loop is certainty only faster when nmemb exceeds some platform dependent threshold n. Consider steering code based on that value. Use a single simple loop to replicate n elements and then after that, double with memcpy()/memove(). This echoes the idea of @Emily L. to measure performance. With this approach, OP may have discovered the value of n is very large rendering the exercise moot.

2. Expand functionality: I'd would expect calloc_d(size_t nmemb) to perform that same, with a ret[0] = 0.0;, ret[0] = 42.0; or ret[0] = DBL_MAX;. So using calloc_d(size_t nmemb, double initial_value) { ... ret[0] = initial_value; makes this function far more useful.

3. memcpy() is usually a bit faster than memmove(), but that difference is more significant with smaller n and #1 above suggest only using memcpy()/memmove() when n is large. Do not expect significant difference in using either of these functions when combined with step #1.

In the end, I would go with usingcalloc(), which zeros the bit data and a while() loop for rare machines that do not encode 0.0 as all zero bits.

• Thanks for the comments, @chux. I considered designing calloc_d() to take an initialization value, but decided not to complicate the code example, and thought that zero-initialization was likely the most useful. I was surprised at the actual performance. For -O0 versions, mine appears a bit faster, but for the -O3 versions the naive approach appears significantly faster. There may be issues with my timing code; links above if you are interested. – ex nihilo Dec 19 '16 at 18:36

### Footnote

After reading all the answers, I decided to see what would be the optimal way to initialize the array and how close the various answers came to the optimal time. Unfortunately, I'm working on a 32-bit target so I used float and uint32_t to test instead of double and uint64_t. But my code can easily be modified if someone wants to test it on a 64-bit target.

Essentially, all I did was to hand code a version that used rep stosl to do the memset. This turned out to be the fastest time:

Time to initialize array of 200 million floats:

memcpy (OP's function)   : 236 msec
loop (Emily L's function): 201 msec (187 msec with -sse3)
rep stosl                : 178 msec


### Update: 64-bit target, sse optimizations

For the above, numbers, I found out that the gcc I was using did not turn on sse optimizations by default. When I turned those on, the loop version dropped to 187 msec.

I then found a 64-bit x86 machine to test on. What I found was that the simple loop actually outperformed the "rep stosq" method. I used gcc -O3 (and sse optimizations were on by default on this machine). Here are the results, using doubles this time:

memcpy (OP's function)   : 845 msec
loop (Emily L's function): 615 msec
rep stosq                : 682 msec


FYI, the assembly code for the simple loop was a 256-byte at a time unrolled loop that looked something like this:

loop_start:
movupd %xmm1, -240(%rdx)
movupd %xmm1, -224(%rdx)
movupd %xmm1, -208(%rdx)
... snip ...
movupd %xmm1, (%rdx)
addq   $256, %rdx addq$-32, %rcx
jne    loop_start


So I guess trusting the compiler is probably the best thing to do.

### Benchmarking code

Here is the program I used to generate the above benchmarks:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <sys/time.h>

void memset_f0(float *array, float val, size_t nmemb)
{
size_t init_sz;       // size of initialized segment
float *next = array; // pointer to beginning of uninitialized segment
size_t numBytes = sizeof(*array) * nmemb;

array[0] = val;
init_sz = sizeof(*array);
++next;

while (init_sz < (numBytes + sizeof(*array)) / 2) {
memmove(next, array, init_sz);
init_sz *= 2;
next = array + init_sz / sizeof(*array);
}
memmove(next, array, numBytes - init_sz);
}

void memset_f1(float *array, float val, size_t nmemb)
{
size_t i;

for (i=0;i<nmemb;i++)
array[i] = val;
}

void memset_f2(float *array, float val, size_t nmemb)
{
union {
uint32_t intVal;
float    fltVal;
} u;

u.fltVal = val;
__asm__ __volatile__ ("cld\n"    \
"rep\n"    \
"stosl\n"  \
: "+a"(u.intVal), "+D"(array), "+c"(nmemb));
}

int getTimeMsec(void)
{
struct timeval tv;

gettimeofday(&tv, NULL);
return (tv.tv_sec * 1000 + tv.tv_usec / 1000);
}

int main(int argc, char *argv[])
{
int size = 200000000;
int method = 0;
float *array;
int begin, end;

if (argc > 1)
method = atoi(argv[1]);

if (argc > 2)
size = atoi(argv[1]);

array = calloc(size, sizeof(float));

begin = getTimeMsec();
switch (method) {
case 0:
memset_f0(array, (float) size, size);
break;
case 1:
memset_f1(array, (float) size, size);
break;
case 2:
memset_f2(array, (float) size, size);
break;
}
end = getTimeMsec();
printf("Elapsed time for size %d method %d: %d msec\n", size, method,
end - begin);
return 0;
}

• NIce. Thanks for this. Did you optimize with -O3 for this? And is there any reason that using floats instead of doubles might change the results? – ex nihilo Dec 20 '16 at 22:50
• @DavidBowling Check the update to this answer. – JS1 Dec 21 '16 at 4:12
• As a note, the code you're testing isn't equivalent to OPs. As val is not a constant 0 the optimisations performed by the compiler will be different. Still an interesting test. An added benefit of the loop is that if the compiler improves so might the performance of the loop while the other two will have the same performance. – Emily L. Dec 21 '16 at 19:32

If you are worried that 0.0 might not be all bits zero, then it's absolutely fine to write a function like calloc_d that will return a double* or float*, with all array elements set to zero. If you are also worried about putting in extra work that isn't needed, you can make that function just call calloc - if it breaks, there is a single function to change, and not possibly 100 calls all over the place, so you have already gained a lot.

If you implement filling memory with memcpy, don't copy too many bytes at a time. I'd say copy no more than 20,000 bytes at a time, so that the bytes you are copying are in the fastest cache.

• Thanks for the reminder about the number of bytes copied in each iteration. – ex nihilo Dec 19 '16 at 13:48