multfly.c

// distributed under public domain
// work in progress

#include <stdint.h>
#include <string.h>

static int multfly_rotl(uint32_t k, int n) {
	return (k << n) | (k >> (32 - n));
}

static const uint32_t multfly_blake2s_iv[8] = {
	UINT32_C(0x6A09E667), UINT32_C(0xBB67AE85), UINT32_C(0x3C6EF372), UINT32_C(0xA54FF53A),
	UINT32_C(0x510E527F), UINT32_C(0x9B05688C), UINT32_C(0x1F83D9AB), UINT32_C(0x5BE0CD19)
};

static const uint8_t multfly_blake2s_sigma[20][8] = {
	{ 0,  1,  2,  3,  4,  5,  6,  7},
	{ 8,  9, 10, 11, 12, 13, 14, 15},
	{14, 10,  4,  8,  9, 15, 13,  6},
	{ 1, 12,  0,  2, 11,  7,  5,  3},
	{11,  8, 12,  0,  5,  2, 15, 13},
	{10, 14,  3,  6,  7,  1,  9,  4},
	{ 7,  9,  3,  1, 13, 12, 11, 14},
	{ 2,  6,  5, 10,  4,  0, 15,  8},
	{ 9,  0,  5,  7,  2,  4, 10, 15},
	{14,  1, 11, 12,  6,  8,  3, 13},
	{ 2, 12,  6, 10,  0, 11,  8,  3},
	{ 4, 13,  7,  5, 15, 14,  1,  9},
	{12,  5,  1, 15, 14, 13,  4, 10},
	{ 0,  7,  6,  3,  9,  2,  8, 11},
	{13, 11,  7, 14, 12,  1,  3,  9},
	{ 5,  0, 15,  4,  8,  6,  2, 10},
	{ 6, 15, 14,  9, 11,  3,  0,  8},
	{12,  2, 13,  7,  1,  4, 10,  5},
	{10,  2,  8,  4,  7,  6,  1,  5},
	{15, 11,  9, 14,  3, 12, 13 , 0}
};

static void multfly_blake2s_qround(uint32_t v[16], const uint32_t m[16], int a, int b, int c, int d, int r, int qr) {
	v[a] = v[a] + v[b] + m[multfly_blake2s_sigma[r][2 * qr]];
	v[d] = multfly_rotl(v[d] ^ v[a], 32 - 16);
	v[c] = v[c] + v[d];
	v[b] = multfly_rotl(v[b] ^ v[c], 32 - 12);
	v[a] = v[a] + v[b] + m[multfly_blake2s_sigma[r][2 * qr + 1]];
	v[d] = multfly_rotl(v[d] ^ v[a], 32 - 8);
	v[c] = v[c] + v[d];
	v[b] = multfly_rotl(v[b] ^ v[c], 32 - 7);
}

static void multfly_blake2s_dround(uint32_t v[16], const uint32_t m[16], int r) {
	multfly_blake2s_qround(v, m, 0, 4,  8, 12, r, 0);
	multfly_blake2s_qround(v, m, 1, 5,  9, 13, r, 1);
	multfly_blake2s_qround(v, m, 2, 6, 10, 14, r, 2);
	multfly_blake2s_qround(v, m, 3, 7, 11, 15, r, 3);
	multfly_blake2s_qround(v, m, 0, 5, 10, 15, r + 1, 0);
	multfly_blake2s_qround(v, m, 1, 6, 11, 12, r + 1, 1);
	multfly_blake2s_qround(v, m, 2, 7,  8, 13, r + 1, 2);
	multfly_blake2s_qround(v, m, 3, 4,  9, 14, r + 1, 3);
}

void multfly_blake2s_compress(uint32_t h[8], const uint32_t m[16], uint64_t cnt, uint32_t fin0) {
	uint32_t v[16];

	for (int i = 0; i < 8; i++) {
		v[i] = h[i];
		v[i + 8] = multfly_blake2s_iv[i];
	}
	v[12] ^= (uint32_t)cnt;
	v[13] ^= (uint32_t)(cnt >> 32);
	v[14] ^= fin0;

	for (int r = 0; r < 20; r += 2) {
		multfly_blake2s_dround(v, m, r);
	}

	for (int i = 0; i < 8; i++) {
		h[i] ^= v[i] ^ v[i + 8];
	}
}

void multfly_init(uint32_t key[4][2], const char *seed, size_t len) {
	// Blake2b
	uint32_t *h = (uint32_t *)key;
	uint32_t m[16];

	for (int i = 0; i < 8; i++) {
		h[i] = multfly_blake2s_iv[i];
	}
	h[0] ^= 0x01010020;

	uint64_t cnt = 0;
	while (len - cnt > 64) {
		for (int i = 0; i < 16; i++) {
			m[i] =
				(uint32_t)seed[cnt + i * 4] |
				(uint32_t)seed[cnt + i * 4 + 1] << 8 |
				(uint32_t)seed[cnt + i * 4 + 2] << 16 |
				(uint32_t)seed[cnt + i * 4 + 3] << 24;
		}
		cnt += 64;
		multfly_blake2s_compress(h, m, cnt, 0);
	}

	for (int i = 0; i < 16; i++) {
		m[i] = 0;
	}
	for (int j = 0; cnt < len; cnt++, j++) {
		m[j / 4] |= (uint32_t)seed[cnt] << (j % 4 * 8);
	}
	multfly_blake2s_compress(h, m, cnt, -1);
}

void multfly_gen(const uint32_t key[4][2], uint64_t ctr, uint32_t result[4]) {
	uint32_t mulu = UINT32_C(2718281829);
	uint32_t mulv = UINT32_C(3141592653);
	uint32_t incr = UINT32_C(0x33123456);

	uint32_t us[4];
	uint32_t vs[4][4];

	for (int i = 0; i < 4; i++) {
		us[i] = key[i][0] ^ (uint32_t)ctr;
		vs[0][i] = key[i][1] ^ (uint32_t)(ctr >> 32);
		++ctr;
	}

	for (int r = 0; r < 3; r++) {
		for (int i = 0; i < 4; i++) {
			uint32_t u = us[i];
			uint32_t v = vs[r][i ^ r];
			u += multfly_rotl(incr, i);
			v += u;
			v ^= multfly_rotl(u, 8);
			v *= mulv;
			u ^= multfly_rotl(v, 9);
			u += v;
			u *= mulu;
			v ^= multfly_rotl(u, 16);
			v += u;
			us[i] = u;
			vs[r + 1][i] = v;
		}
	}

	for (int i = 0; i < 4; i++) {
		result[i] = vs[3][i];
	}
}

void multfly_gen_to_state(const uint32_t key[4][2], uint64_t ctr, uint32_t state[4][4]) {
	uint32_t result[4];
	multfly_gen(key, ctr, result);

	uint32_t incr = UINT32_C(0x33123456);
	for (int i = 0; i < 4; i++) {
		for (int j = 0; j < 4; j++) {
			state[i][j] = result[(j + i) % 4];
		}
		state[i][1] += multfly_rotl(incr, i);
	}
}

void multfly_gen_by_state(uint32_t state[4][4], uint32_t result[4]) {
	// Xoshiro128**
	for (int i = 0; i < 4; i++) {
		result[i] = multfly_rotl(state[i][0] * 5, 7) * 9;
		uint32_t t = state[i][1] << 9;
		state[i][2] ^= state[i][0];
		state[i][3] ^= state[i][1];
		state[i][1] ^= state[i][2];
		state[i][0] ^= state[i][3];
		state[i][2] ^= t;
		state[i][3] = multfly_rotl(state[i][3], 11);
	}
}

const static char seed[] = "Hello, World!";
const static uint64_t initial_ctr = UINT64_C(0x123456789ABCDEF0);

uint32_t multfly() {
	static uint32_t key[4][2];
	static uint64_t ctr;
	static uint32_t result[4];

	static int running = 0;
	if (!running) {
		running = 1;
		ctr = initial_ctr / 4 * 4;
		multfly_init(key, seed, strlen(seed));
	}

	if (ctr % 4 == 0) {
		multfly_gen(key, ctr, result);
	}
	return result[(ctr++) % 4];
}