ZrHa_update avalanche diagram

ZrHa-avalanche.c

// Copyright (c) 2019 Alexey Tourbin
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

#include <stdio.h>
#include <stdint.h>
#include <emmintrin.h>

// pseudo-random numbers from a counter, by Pelle Evensen
static inline uint64_t rrmxmx(uint64_t x)
{
#define rotr64(x, k) (((x) >> (k)) | ((x) << (64 - (k))))
    x ^= rotr64(x, 49) ^ rotr64(x, 24);
    x *= UINT64_C(0x9FB21C651E98DF25);
    x ^= x >> 28;
    x *= UINT64_C(0x9FB21C651E98DF25);
    x ^= x >> 28;
    return x;
}

static inline void ZrHa_update(void *state_, const void *data)
{
    __m128i state = _mm_loadu_si128(state_);
    __m128i x = _mm_add_epi64(state, _mm_loadu_si128(data));
    __m128i m = _mm_mul_epu32(x, _mm_shuffle_epi32(x, _MM_SHUFFLE(2, 3, 0, 1)));
    state = _mm_add_epi64(m, _mm_shuffle_epi32(x, _MM_SHUFFLE(0, 1, 2, 3)));
    _mm_storeu_si128(state_, state);
}

int main()
{
    // avalanche matrix
    int m[128][128] = { 0, };
#define KLOG 10
    for (int k = 1; k <= (1<<KLOG); k++) {
	__uint128_t state0 = rrmxmx(+k) | 
	       (__uint128_t) rrmxmx(-k) << 64;
	// flip each input bit
	for (int i = 0; i < 128; i++) {
	    __uint128_t state = state0;
	    __uint128_t data = (__uint128_t) 1 << i;
	    ZrHa_update(&state, &data);
	    // see which bits got flipped in the output
	    __uint128_t diff = state0 ^ state;
	    for (int j = 0; j < 128; j++)
		m[i][j] += (diff >> j) & 1;
	}
    }
    // PGM image, XY-plane
    printf("P2\n128 128\n255\n");
    for (int j = 127; j >= 0; j--) {
	for (int i = 0; i < 128; i++)
	    printf("%*d", 3 + (i > 0) , m[i][j] >> (KLOG-8));
	putchar('\n');
    }
    return 0;
}

ZrHa-avalanche.png