#include <stdio.h>
#include <stdlib.h>
#include <time.h>
typedef unsigned long long u64;
typedef signed short s16;

const int N = 64;

u64 A[N*N];

u64 RndUniform() {
	return (((u64)rand()) << 32ull) ^ rand();
}

s16 RndNormal() {
	u64 x = RndUniform(), y = 0;
	y += x & 16383; x >>= 14;
	y += x & 16383; x >>= 14;
	y += x & 16383; x >>= 14;
	y += x & 16383; x >>= 14;
	return (x ^ y) - 32768;
}

int main() {
	srand(time(0) ^ rand());
	for(int i=0;i<N*N;i++) A[i] = RndUniform();

	s16 p[N];
	for(int i=0;i<N;i++) {
		p[i] = RndNormal();
	}
	printf("p = [%d, %d, %d, ...] (%d bits)\n", p[0], p[1], p[2], 16*N);

	u64 k[N];
	for(int i=0;i<N;i++) {
		k[i] = RndNormal();
		for(int j=0;j<N;j++) {
			k[i] += A[j*N + i] * p[j];
		}
	}
	printf("k = %016llX... (%d bits)\n", k[0], 64*N);

	u64 m[] = {2147483648, 2147483648, 2147483648, 2147483648};
	u64 c[N*4];
	for(int d=0;d<4;d++) {
		u64 r[N];
		for(int i=0;i<N;i++) {
			r[i] = RndNormal();
		}
		for(int i=0;i<N;i++) {
			c[i + N*d] = RndNormal();
			for(int j=0;j<N;j++) {
				c[i + N*d] += A[i*N + j] * r[j];
			}
		}
		for(int j=0;j<N;j++) {
			m[d] -= k[j] * r[j];
		}
		m[d] >>= 32;
	}
	printf("G(r) = %016llX... (%d bits)\n", c[0], 64*N*4);
	printf("FG(r)   = %08llX%08llX%08llX%08llX\n", m[0], m[1], m[2], m[3]);

	u64 v[] = {2147483648, 2147483648, 2147483648, 2147483648};
	for(int d=0;d<4;d++) {
		for(int j=0;j<N;j++) {
			v[d] -= p[j] * c[j + N*d];
		}
		v[d] >>= 32;
	}
	printf("F(G(r)) = %08llX%08llX%08llX%08llX\n", v[0], v[1], v[2], v[3]);
}