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

// GENERAL /////////////////////////////////////////////////////////////
#define QU 81
//#define QU 27
//#define QU 9
#define  DEGREE 155
//#define  DEGREE 68
//#define  DEGREE 12
#define  Nm 	531441
//#define  Nm 	19683
//#define  Nm 	729
#define  ROWm	QU+5

typedef long i64;
typedef unsigned long ui64;
typedef unsigned int uint128_t __attribute__((mode(TI)));

#define align __attribute__ ((aligned (32)))

// r = 33098280119090191028775580055082175056428495623;
// r has a size of 155 bits

// FIELD GF(2^233) /////////////////////////////////////////////////////


typedef align ui64 elt[3];
//typedef align ui64 elt_v[256][3];
typedef struct elt_v{
   elt x;
} elt_v;

typedef align ui64 *elt_pt;

void types_randgen(elt ax, int deg);

elt_v V[2][Nm];
elt_v Vo[Nm];
	
#define types_print(X)													\
printf("%s : %8.16lX %8.16lX %8.16lX\n", #X, X[2], X[1], X[0]);
//printf("%s : %8.16li %8.16li %8.16li\n", #X, X[2], X[1], X[0]);

#define types_fprint(F, X)													\
fprintf(F, "0x%8.16lX%8.16lX%8.16lX", X[2], X[1], X[0]);


#define types_print_double(X)													\
printf("%s : %8.16lX %8.16lX %8.16lX %8.16lX %8.16lX %8.16lX\n", #X, X[5], X[4], X[3], X[2], X[1], X[0]);


#define types_comp(a, b)												\
	(memcmp(a, b, sizeof(ui64)*3)==0)

#define types_setzero(a)												\
	memset(a, 0, sizeof(ui64)*3);
	
#define types_setone(a)													\
	memset(a, 0, sizeof(ui64)*3);										\
	a[0] = 0x1;

#define types_isone(a)	(a[0]==1 && a[1]==0 && a[2]==0) 

#define types_iszero(a)	(!(a[0] || a[1] || a[2])) 

#define types_iseven(a) !(a[0]&0x1)

#define types_copy(Y,X)													\
	memcpy(Y, X, sizeof(ui64)*3);

#define types_singleshiftright(ax) {\
	int k, sshift = 64-1;\
	for (k=0; k<2; k++) {\
			ax[k] >>= 1;\
			ax[k]  ^= (ax[k+1]&1)<<sshift;\
	}\
	ax[k] >>= 1;\
}

// Definition of the arithmetic functions
void types_randgen(elt ax, int deg);
int types_less(elt a, elt b );
int in_addN(ui64 *out, const ui64 *x, const ui64 *y); 
inline  Fp_addC(elt out, ui64 *x, ui64 *y, elt p);
inline in_subN(ui64 *out, const ui64 *x, const ui64 *y, ui64 s);
static inline void Fp_reduce(elt out, ui64 *x, elt p);
void Fp_subC(elt out, ui64 *x, ui64 *y, elt p);
void m_E(ui64* out, ui64 a[], ui64 b[], ui64 s, ui64 t);
void r_E_small_2113(elt out, ui64 in[], elt p);
void r_E_small_2887(elt out, ui64 in[], elt p);
void r_E_165301473942399079669(elt out, ui64 in[], elt p);
void Fp_mulC_s(elt out, ui64 a[], ui64 b[], elt p);
void r_E(elt out, ui64 in[], elt p);
void Fp_mulC(elt out, ui64 a[], ui64 b[], elt p);
void Fp_invC(elt out, elt a, elt p);
int BerlekampMassey_s(char c[], char s[], int n);
int BerlekampMassey(elt_v C[], elt_v S[], elt p, int N);