bday.c

#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <stdint.h>
#include <pthread.h>
#include <stdlib.h>

typedef struct {
    int N, M;
} bday_arg;

// *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)

typedef struct { uint64_t state;  uint64_t inc; } pcg32_random_t;

uint32_t pcg32_random_r(pcg32_random_t* rng)
{
    uint64_t oldstate = rng->state;
    // Advance internal state
    rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
    // Calculate output function (XSH RR), uses old state for max ILP
    uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
    uint32_t rot = oldstate >> 59u;
    return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}


void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq)
{
    rng->state = 0U;
    rng->inc = (initseq << 1u) | 1u;
    pcg32_random_r(rng);
    rng->state += initstate;
    pcg32_random_r(rng);
}


#define likely(p) __builtin_expect(!!(p), 1)
#define unlikely(p) __builtin_expect(!!(p), 0)

uint64_t nearlydivisionless ( pcg32_random_t* rng, uint64_t s ) {
  uint32_t x = pcg32_random_r (rng) ;
  uint64_t m = ( uint64_t ) x * ( uint64_t ) s;
  uint32_t l = ( uint32_t ) m;
  if (unlikely(l < s)) {
    uint32_t t = -s % s;
    while (l < t) {
      x = pcg32_random_r (rng) ;
      m = ( uint64_t ) x * ( uint64_t ) s;
      l = ( uint32_t ) m;
    }
  }
  return m >> 32;
}

void *bday_thread(void *arg_in) {
    bday_arg *arg = arg_in;
    _Bool days[365];

    pcg32_random_t rng;

    uint64_t seed[2];
    getentropy(seed, sizeof(seed)); 
    pcg32_srandom_r(&rng, seed[0], seed[1]);

    int N = arg->N;
    int M = arg->M;
    int coll = 0;
    for(int i=0; i<N; i++) {
        memset(days, 0, sizeof(days));
        for(int j=0; j<M; j++) {
            int r = nearlydivisionless(&rng, 365);
            if (unlikely(days[r])) { 
                coll ++;
                break;
            }
            days[r] = 1;
        }
    }
    return (void*)(intptr_t)coll;
}

#define NN (399168000)

int main(int argc, char **argv) {
    int coll = 0;
    if (argc > 1) {
        int TT = atoi(argv[1]);
        bday_arg bb[TT];
        pthread_t th[TT];
        for(int i=0; i<TT; i++) {
            bb[i] = (bday_arg){ NN / TT, 24 };
            pthread_create(&th[i], NULL, bday_thread, bb);
        }
        for(int i=0; i<TT; i++) {
            void *res;
            pthread_join(th[i], &res);
            coll += (intptr_t)res;
        }
    } else {
        bday_arg b = {NN, 24 };
        coll = (intptr_t) bday_thread(&b);
    }
    printf("%f\n", coll * 100./ NN);
}