// This program is an attachment of the article 
// An infinite family of MUB-triplets in dimension 6
// written by P. Jaming, M. Matolcsi, P. M\'{o}ra,
// F. Sz\"{o}ll\H{o}si and M. Weiner.
//
// This is a C++ program. I compile it with gcc compiler with
//
//  g++ -O3 fab_uvb.cpp -o fab_uvb
//
// command, no other files are required. There
// should be no problem compile it with a 32/64 bit compiler.
// 
// This program is a part of a proof. In order to repeat this
// proof run this program without any parameters
// 
//  ./fab_uvb
//
// However this program assumes that ort_19.txt exists and
// created by the program ort.cpp. Please compile and run
// ort.cpp first.
//
// Under Windows it can be compiled with Dev-Cpp (free compiler,
// http://www.bloodshed.net/devcpp.html). After opening this file
// you can compile and run this program with button F9. This program
// requires no command line parameters (otherwise it returns with
// error but you might not be able to read it because it closes too
// fast). Please ensure that under menu Execute-> Parameters... 
// the line "Parameters to pass to your program" is empty.
//
//
// If this program runs successfully then some output messages
// and vectors are written to fab_uvb.txt. It takes about six hours
// to run on a computer with a 3,2 GHz CPU.
//
// The aim of this program is to prove that the two parameter
// Fourier family F(x,y) cannot be extended to a MUB-quartet.
// We assume that the matrix A is the identity matrix.
//
// In the paper vectors are considered as columns of matrices. However, 
// please note that here in the program we deal with them as rows. 
//
// The set of vectors found_ubvs in this program is denoted by 
// $FUB_{N, N'}^{\tilde a, \tilde b}$ in the paper. 
//
// In the following the vector ((0,q,w,e,r,t)) at level m_n2
// represents a set of vectors in dimension 6, which first
// coordinate is 0 and the 2nd, 3rd, 4th, 5th, 6th coordinates
// are in the intervals [q/m_n2,(q+1)/m_n2], [w/m_n2,(w+1)/m_n2],
// [e/m_n2,(e+1)/m_n2], [r/m_n2,(r+1)/m_n2],
// [t/m_n2,(t+1)/m_n2], respectively.
//
// First of all we fix the parameters 'a' and 'b'. They are integers
// and the two parameters of the Fourier family are in
// the intervals [a/m_n1,(a+1)/m_n1] and [b/m_n1,(b+1)/m_n1],
// respectively. We call this matrix B (some of the elements of
// this matrix are therefore intervals).
//
// Then we list out all possible vectors ((0,q,w,e,r,t)) at level m_n2
// which are unbiased to the matrix B.
// We say that the vector ((0,q,w,e,r,t)) at level m_n2
// is unbiased to a matrix if there exists a vector in the set
// defined by ((0,q,w,e,r,t)) and this vector is unbiased to every
// row of the matrix (after applying t->e^{i*2*Pi*t} in all coordinates).
//
// Further, we try to search rows for the matrix C from the vectors which are
// unbiased to the matrix B.
//
// After that we try to search rows for the matrix D, but we cannot build up
// any matrix D which would meet all conditions. This
// completes the proof.
// 

#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <vector>
#include <time.h>
#include <algorithm>

using namespace std;

// For every calculation we use double as floating type.
// The precision of a double type variable is about 16 decimal digits,
// therefore all round-off error are less than EPSILON:=10^{-9}.
#define EPSILON 0.000000001
#define PI2 6.28318530717958647692

// The array m_shift contains 5*32 values. 
// There is 2^5 ways to make a 0-1 list with length 5, m_shift consists of all of these.
// Namely:
// 0, 0, 0, 0, 0
// 0, 0, 0, 0, 1
// 0, 0, 0, 1, 0
// 0, 0, 0, 1, 1
// ...
// 1, 1, 1, 1, 1
int     m_shift [] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1};

// The array m_perm contains all possible permutation of the set {1,2,3,4,5}.
// Namely:
// 1, 2, 3, 4, 5
// 1, 2, 3, 5, 4
// ...
// 5, 4, 3, 2, 1
int     m_perm [] = {1, 2, 3, 4, 5, 1, 2, 3, 5, 4, 1, 2, 4, 3, 5, 1, 2, 4, 5, 3, 1, 2, 5, 3, 4, 1, 2, 5, 4, 3, 1, 3, 2, 4, 5, 1, 3, 2, 5, 4, 1, 3, 4, 2, 5, 1, 3, 4, 5, 2, 1, 3, 5, 2, 4, 1, 3, 5, 4, 2, 1, 4, 2, 3, 5, 1, 4, 2, 5, 3, 1, 4, 3, 2, 5, 1, 4, 3, 5, 2, 1, 4, 5, 2, 3, 1, 4, 5, 3, 2, 1, 5, 2, 3, 4, 1, 5, 2, 4, 3, 1, 5, 3, 2, 4, 1, 5, 3, 4, 2, 1, 5, 4, 2, 3, 1, 5, 4, 3, 2, 2, 1, 3, 4, 5, 2, 1, 3, 5, 4, 2, 1, 4, 3, 5, 2, 1, 4, 5, 3, 2, 1, 5, 3, 4, 2, 1, 5, 4, 3, 2, 3, 1, 4, 5, 2, 3, 1, 5, 4, 2, 3, 4, 1, 5, 2, 3, 4, 5, 1, 2, 3, 5, 1, 4, 2, 3, 5, 4, 1, 2, 4, 1, 3, 5, 2, 4, 1, 5, 3, 2, 4, 3, 1, 5, 2, 4, 3, 5, 1, 2, 4, 5, 1, 3, 2, 4, 5, 3, 1, 2, 5, 1, 3, 4, 2, 5, 1, 4, 3, 2, 5, 3, 1, 4, 2, 5, 3, 4, 1, 2, 5, 4, 1, 3, 2, 5, 4, 3, 1, 3, 1, 2, 4, 5, 3, 1, 2, 5, 4, 3, 1, 4, 2, 5, 3, 1, 4, 5, 2, 3, 1, 5, 2, 4, 3, 1, 5, 4, 2, 3, 2, 1, 4, 5, 3, 2, 1, 5, 4, 3, 2, 4, 1, 5, 3, 2, 4, 5, 1, 3, 2, 5, 1, 4, 3, 2, 5, 4, 1, 3, 4, 1, 2, 5, 3, 4, 1, 5, 2, 3, 4, 2, 1, 5, 3, 4, 2, 5, 1, 3, 4, 5, 1, 2, 3, 4, 5, 2, 1, 3, 5, 1, 2, 4, 3, 5, 1, 4, 2, 3, 5, 2, 1, 4, 3, 5, 2, 4, 1, 3, 5, 4, 1, 2, 3, 5, 4, 2, 1, 4, 1, 2, 3, 5, 4, 1, 2, 5, 3, 4, 1, 3, 2, 5, 4, 1, 3, 5, 2, 4, 1, 5, 2, 3, 4, 1, 5, 3, 2, 4, 2, 1, 3, 5, 4, 2, 1, 5, 3, 4, 2, 3, 1, 5, 4, 2, 3, 5, 1, 4, 2, 5, 1, 3, 4, 2, 5, 3, 1, 4, 3, 1, 2, 5, 4, 3, 1, 5, 2, 4, 3, 2, 1, 5, 4, 3, 2, 5, 1, 4, 3, 5, 1, 2, 4, 3, 5, 2, 1, 4, 5, 1, 2, 3, 4, 5, 1, 3, 2, 4, 5, 2, 1, 3, 4, 5, 2, 3, 1, 4, 5, 3, 1, 2, 4, 5, 3, 2, 1, 5, 1, 2, 3, 4, 5, 1, 2, 4, 3, 5, 1, 3, 2, 4, 5, 1, 3, 4, 2, 5, 1, 4, 2, 3, 5, 1, 4, 3, 2, 5, 2, 1, 3, 4, 5, 2, 1, 4, 3, 5, 2, 3, 1, 4, 5, 2, 3, 4, 1, 5, 2, 4, 1, 3, 5, 2, 4, 3, 1, 5, 3, 1, 2, 4, 5, 3, 1, 4, 2, 5, 3, 2, 1, 4, 5, 3, 2, 4, 1, 5, 3, 4, 1, 2, 5, 3, 4, 2, 1, 5, 4, 1, 2, 3, 5, 4, 1, 3, 2, 5, 4, 2, 1, 3, 5, 4, 2, 3, 1, 5, 4, 3, 1, 2, 5, 4, 3, 2, 1};


// The pointers x and y are supposed to be the addresses of arrays of integers with length 6.
// The function less_set returns true if x[0], x[1], x[2], x[3], x[4], x[5] is less than
// y[0], y[1], y[2], y[3], y[4], y[5] lexicography.
struct less_set  : public binary_function<int*, int*, bool> {
           bool operator()(int* x, int* y) {
                   if (x[0] < y[0])
                           return true;
                   if (x[0] > y[0])
                           return false;
                   if (x[1] < y[1])
                           return true;
                   if (x[1] > y[1])
                           return false;
                   if (x[2] < y[2])
                           return true;
                   if (x[2] > y[2])
                           return false;
                   if (x[3] < y[3])
                           return true;
                   if (x[3] > y[3])
                           return false;
                   if (x[4] < y[4])
                           return true;
                   if (x[4] > y[4])
                           return false;
                   if (x[5] < y[5])
                           return true;
                   else
                           return false;
}};

// The function itostr converts an integer to string (for the type string see Standard Template Library).
string   itostr(int i)
{
        char s[50];
        sprintf(s,"%d",i);
        return string(s);
};


// All the computations are done by the class Mub. 
// We will declare one instance of it.
class Mub{

	// We consider a family of Hadamard matrices with
	// two parameters. These parameters are chosen from
	// the interval [0,1]. We use a finite approximation:
	// the first parameter is in the interval [a/m_n1,(a+1)/m_n1],
	// the second parameter is in [b/m_n1,(b+1)/m_n1], where
	// the variables "a" and "b" will be used in function start.
	int		m_n1;
	
	// We search for the matrices C and D. The rows of these matrices
	// are chosen from the array found_ubvs (see function start) of vectors.
	// All of these vectors have the form: 
	// (0,x_1,x_2,x_3,x_4,x_5) at level m_n2, where x_1, x_2, ..., x_5
	// are integers. It means that the (i+1)-th coordinate of the
	// vector is in the interval [x_i/m_n2,(x_i+1)/m_n2] and the first
	// coordinate is 0.
	int		m_n2;

	// m_max_n2 = m_n2 * 2 ^ deep, where ^ means the power function
	// and deep is the parameter of the constructor Mub.
	int		m_max_n2;

	// The array m_hadamard contains the second, third, fourth, fifth
	// and sixth rows of the matrix B. Some values of it are the correct
	// values, for example x/m_n2, others means intervals [x/m_n1,(x+1)/m_n1].
	int		m_hadamard[30];

	// We search for all those vectors which are unbiased to the matrix B.
	// To do so, we use the array m_row. The variable m_row[0] is always 0.
	int		m_row[6];

	// The function scalar returns true if the current value of m_row at level n 
	// is unbiased to the matrix B. m_row[0]=0, which means that the
	// first coordinate is 0. For i=1,2,3,4,5 the corresponding second, third,
	// fourth, fifth, sixth coordinates are in [m_row[1]/n,(m_row[1]+1)/n],
	// ..., [m_row[5]/n,(m_row[5]+1)/n], where n is the parameter of this function.
	bool		scalar(int n);


	// For error estimate and scalar product we need to calculate a lot of values
	// many times. To speed up the program we store these values in advance.

	// For i=1,2,...,m_max_n2 we have
	//  m_error_with_sure[i] =  (double)5 * PI2 / (double) (2*i);
	//  m_error_with_unsure[i] =  ((double)5 * PI2 / (double) (2*i) + (double)4 * PI2 / (double) (2*m_n1));
	double*		m_error_with_unsure;		
	double*		m_error_with_sure;

	// For i=0,1,...,4*m_max_n2 we have
	//  ssin[i] = sin(((double)i) * PI2 / (double)(2*m_max_n2));
	//  ccos[i] = cos(((double)i) * PI2 / (double)(2*m_max_n2));
	double*		ssin;
	double*		ccos;

	// The elements of the array m_hadamard represents either a point or an interval
	// of [0,1] (for i=1,2,4,5,13,14,16,17,25,26,28,29 it is an interval, otherwise it is
	// a point, for more details see function start). We consider the point as an
	// interval with zero length. For i=0,1,...,29 hcos[3*i+0], hsin[3*i+0] equal to
	// the real and the imaginary part of e^{2*Pi*i*x}, where x is the left end of
	// the interval defined by m_hadamard[i]. In all cases x=m_hadamard[i]/m_n1.
	//
	// hcos[3*i+1], hsin[3*i+1] are the real and imaginary part of e^{2*Pi*i*y}, where
	// y is the middle of the interval defined by m_hadamard[i].
	//
	// hcos[3*i+2], hsin[3*i+2] are the real and imaginary part of e^{2*Pi*i*z}, where
	// z is the right end of the interval defined by m_hadamard[i].
	double		hsin[90];
	double		hcos[90];

	// Sqrt(6)
	double		m_sqrt6;

	// For the correct proof we do not need the array deepest and the corresponding code. 
	// It is only interesting for checking the correctness of this program.
	int		m_deepest[6];

	// The matrix B represents matrices from a two parameter family of a Hadamard matrices.
	// Some rows are fixed, and others defines vectors with some restriction for theirs
	// coordinates. The array m_fix says which rows are fixed. In our example:
	// m_fix[0] = m_fix[2] = m_fix[4] = false;
	// m_fix[1] = m_fix[3] = true;
	bool		m_fix[5];

	// The array m_ort_eps contains int[6] type arrays. Let us suppose that x and y are int[6].
	// x[0] = y[0] = 0. They represents two sets of vectors:
	// (0,[x[1]/m_n2,(x[1]+1)/m_n2], ..., [x[5]/m_n2,(x[5]+1)/m_n2])
	// and (0,[y[1]/m_n2,(y[1]+1)/m_n2], ..., [y[5]/m_n2,(y[5]+1)/m_n2]).
	// If there exists a vector from the first set and a vector from the second set in that way
	// they are orthogonal after applying function t->e^{2*Pi*i*t} in every coordinates,
	// then (0,x[1]-y[1] mod m_n2,...,x[5]-y[5] mod m_n2) is in m_ort_eps
	// and the function scalar_with_eps(x,y) returns true. (It may also return true otherwise.)
	// m_ort_eps_database[(((q*m_n2+w)*m_n2+e)*m_n2+r)*m_n2+t] == true if and only if
	// [0,q,w,e,r,t] is in m_ort_eps;
        vector<int*>    m_ort_eps;
        int*            m_ort_eps_database;
	bool		scalar_with_eps(int* v1, int* v2);

	// To speed up we use our modulo for values in [-m_n2,m_n2-1] instead of the built-in
	// function %.
        int             modulo(int);

	// For testing.
	void		print(int* v);

	// After fixing matrices B and C we start to search for vectors from the arrays found_ubvs
	// (see function start) which are unbiased to the matrix C. The parameters:
	// m = C, and the currently tested vector is (0,[v[1]/n,(v[1]+1)/n],...,[v[5]/n,(v[5]+1)/n]).
	bool		ubv_with_matrix(int* m, int* v, int n);

	// The function start calls iterate while we are searching for vectors, which are
	// unbiased to the matrix B. The currently tested vector is 
	// (0,[m_row[1]/n,(m_row[1]+1)/n],...,[m_row[5]/n,(m_row[5]+1)/n]).
	bool		iterate(int n);
	
public:
	Mub(int n, int deep);
	~Mub();

	// We consider a family of Hadamard matrices with
	// two parameters. These parameters are chosen from
	// the interval [0,1]. We use a finite approximation:
	// the first parameter is in the interval [a/m_n1,(a+1)/m_n1],
	// the second parameter is in [b/m_n1,(b+1)/m_n1].
	void		start(int n1, int a, int b, FILE* out);
};

void	Mub::print(int* v)
{
	printf("-  %d %d %d %d %d %d\n",v[0],v[1],v[2],v[3],v[4],v[5]);
}

int     Mub::modulo(int a)
{
        if (a < 0)
                return a+ m_n2;
        return a;
}


bool    Mub::scalar_with_eps(int* v1, int* v2)
{
       	return m_ort_eps_database[(((modulo(v1[1]-v2[1])*m_n2+modulo(v1[2]-v2[2]))*m_n2+modulo(v1[3]-v2[3]))*m_n2+modulo(v1[4]-v2[4]))*m_n2+modulo(v1[5]-v2[5])];
}

Mub::Mub(int n, int deep)
{
	m_n2 = n;
	m_max_n2 = m_n2;
	for (int i = 0; i < deep; i++)
	{
		m_max_n2 *= 2;
	}

        string filename;
        filename = "ort_";
        filename = filename + itostr(m_n2) + ".txt"; // We read from the file "ort_19.txt"
						// if n=19.

        FILE*   ort;
        ort = fopen(filename.c_str(),"r");
	// This file should exist.
        if (!ort)
        {
                printf("The file %s must exist!!! Run the program ort first!\n",filename.c_str());
                exit(1);
        }

	// Calculating often used data in advance.
	ssin = new double[4*m_max_n2];
	ccos = new double[4*m_max_n2];
	for (int i = 0; i < 4*m_max_n2; i++)
	{
		ssin[i] = sin(((double)i) * PI2 / (double)(2*m_max_n2));
		ccos[i] = cos(((double)i) * PI2 / (double)(2*m_max_n2));
	}
	m_error_with_unsure = NULL;
	m_error_with_sure = NULL;
	m_sqrt6 = sqrt(6.0);


        vector<int*> temp;
        vector<int*> temp2;

	

        int d0, d1, d2, d3, d4, d5;
        char st[100];
        // Load lines from ort_19.txt to the array temp (if n=19).
	// These arrays with 6 length contains all of those x[6] "vectors", which
	// defines vectors in (0,[x[1]/m_n2,(x[1]+1)/m_n2],...,[x[5]/m_n2,(x[5]+1)/m_n2]).
	// These vectors might be orthogonal to (0,0,0,0,0,0) after applying
	// t->e^{2*Pi*i*t} in all coordinates.
	//
	// In the following we want to construct the array m_ort_eps.
	// Namely, we want all those vectors v, which have the form
	// v=x+y, where x is orthogonal to (0,0,0,0,0,0) in the same sense
	// as above and y is (0,q,w,e,r,t) where q,w,e,r,t are in the set {0,1}.
        while (fgets(st, 100, ort))
        {
                sscanf(st,"%d, %d, %d, %d, %d, %d\n",&d0,&d1,&d2,&d3,&d4,&d5);
                int* to_load = new int[6];
                to_load[0] = d0; to_load[1]=d1; to_load[2]=d2; to_load[3]=d3; to_load[4]=d4; to_load[5]=d5;
                temp.push_back(to_load);
        }

        //Shift all the 32 ways all vector in temp, and put it in temp2.
        for (int i=0; i < temp.size(); i++)
        {
                for (int j = 0; j < 160; j+=5)
                {
                        int* to_load = new int[6];
                        to_load[0] = temp[i][0];
                        to_load[1] = (temp[i][1]+m_shift[j]) % m_n2;
                        to_load[2] = (temp[i][2]+m_shift[j+1]) % m_n2;
                        to_load[3] = (temp[i][3]+m_shift[j+2]) % m_n2;
                        to_load[4] = (temp[i][4]+m_shift[j+3]) % m_n2;
                        to_load[5] = (temp[i][5]+m_shift[j+4]) % m_n2;
                        temp2.push_back(to_load);
                }
        }
        //sort the array temp2
        printf("Start sorting ..."); cout.flush();
        sort(temp2.begin(),temp2.end(),less_set());
        printf("OK\n"); cout.flush();


        //there might be multiple vectors in temp2
        //we make one copy from all multiple vectors to m_ort_eps
        d0 = -1;
        for (int i =0; i< temp2.size(); i++)
        {
                if (d0 != temp2[i][0] || d1 != temp2[i][1] || d2 != temp2[i][2] || d3 != temp2[i][3] || d4 != temp2[i][4] || d5 != temp2[i][5])
                {
                        m_ort_eps.push_back(temp2[i]);
                        d0 = temp2[i][0];
                        d1 = temp2[i][1];
                        d2 = temp2[i][2];
                        d3 = temp2[i][3];
                        d4 = temp2[i][4];
                        d5 = temp2[i][5];
                }
                else
                        delete[] temp2[i];
        }
        temp2.clear();
        printf("Size of ORTeps: %d\n",m_ort_eps.size());
        cout.flush();


        m_ort_eps_database = new int[m_n2*m_n2*m_n2*m_n2*m_n2];

        printf("generating database for ORTeps... ");
        cout.flush();
	for (int i = 0; i < m_n2*m_n2*m_n2*m_n2*m_n2; i++)
	{
		m_ort_eps_database[i] = false;
	}
        for (int i = 0; i < m_ort_eps.size(); i++)
        {
                m_ort_eps_database[(((m_ort_eps[i][1]*m_n2+m_ort_eps[i][2])*m_n2+m_ort_eps[i][3])*m_n2+m_ort_eps[i][4])*m_n2+m_ort_eps[i][5]] = true;
		delete[] m_ort_eps[i];
        }
	m_ort_eps.clear();
        printf("OK\n");
        cout.flush();

};

Mub::~Mub()
{
	delete[] ssin;
	delete[] ccos;
	delete[] m_ort_eps_database;
};


void	Mub::start(int n1, int a, int b,FILE* out)
{
	m_n1 = n1;


	if (m_n1 % 6 != 0)
	{
		printf("ERROR, m_n1 should be multiply of 6.\n");
		exit(1);
	}


	// The first row of the B matrix is 0,0,0,0,0,0.
	// m_fix[i] = false if and only if the (i+2). row of B
	// contains parameters.
	// The B matrix:

	m_fix[0] = false;
	m_hadamard[0] = 0;
	m_hadamard[1] = 1*m_n1/6+a; 
	m_hadamard[2] = 2*m_n1/6+b; 
	m_hadamard[3] = 3*m_n1/6;
	m_hadamard[4] = 4*m_n1/6+a; 
	m_hadamard[5] = 5*m_n1/6+b; 

	m_fix[1] = true;
	m_hadamard[6] = 0;
	m_hadamard[7] = 2*m_n1/6;
	m_hadamard[8] = 4*m_n1/6;
	m_hadamard[9] = 0*m_n1/6;
	m_hadamard[10] = 2*m_n1/6;
	m_hadamard[11] = 4*m_n1/6;

	m_fix[2] = false;
	m_hadamard[12] = 0;
	m_hadamard[13] = 3*m_n1/6+a; 
	m_hadamard[14] = 0*m_n1/6+b; 
	m_hadamard[15] = 3*m_n1/6;
	m_hadamard[16] = 0*m_n1/6+a; 
	m_hadamard[17] = 3*m_n1/6+b; 

	m_fix[3] = true;
	m_hadamard[18] = 0;
	m_hadamard[19] = 4*m_n1/6;
	m_hadamard[20] = 2*m_n1/6;
	m_hadamard[21] = 0*m_n1/6;
	m_hadamard[22] = 4*m_n1/6;
	m_hadamard[23] = 2*m_n1/6;

	m_fix[4] = false;
	m_hadamard[24] = 0;
	m_hadamard[25] = 5*m_n1/6+a; 
	m_hadamard[26] = 4*m_n1/6+b; 
	m_hadamard[27] = 3*m_n1/6;
	m_hadamard[28] = 2*m_n1/6+a; 
	m_hadamard[29] = 1*m_n1/6+b; 


	// We will need some sin and cos values related to this Hadamard matrix.
	// We store these values in advance.

	for (int i = 0; i < 30; i++)
	{
		hsin[3*i+0] = sin(((double)2*m_hadamard[i]) * PI2 / (double)(2*m_n1));
		hcos[3*i+0] = cos(((double)2*m_hadamard[i]) * PI2 / (double)(2*m_n1));
		if (m_fix[i/6] || i%6 == 0 || i%6==3)
		{
			hsin[3*i+1] = hsin[3*i+2] = hsin[3*i+0];
			hcos[3*i+1] = hcos[3*i+2] = hcos[3*i+0];
		}
		else
		{
			hsin[3*i+1] = sin(((double)2*m_hadamard[i]+1) * PI2 / (double)(2*m_n1));
			hcos[3*i+1] = cos(((double)2*m_hadamard[i]+1) * PI2 / (double)(2*m_n1));
			hsin[3*i+2] = sin(((double)2*m_hadamard[i]+2) * PI2 / (double)(2*m_n1));
			hcos[3*i+2] = cos(((double)2*m_hadamard[i]+2) * PI2 / (double)(2*m_n1));
		}
	}

	// In function scalar we will need the following values.
	m_error_with_sure = new double[m_max_n2+1];
	for (int i=1; i <= m_max_n2; i++)
		m_error_with_sure[i] =  (double)5 * PI2 / (double) (2*i);
	m_error_with_unsure = new double[m_max_n2+1];
	for (int i=1; i <= m_max_n2; i++)
		m_error_with_unsure[i] =  ((double)5 * PI2 / (double) (2*i) + (double)4 * PI2 / (double) (2*m_n1));


	// Searching for vectors which are unbiased to B.

	m_row[0] = 0;
	int d = m_n2;
	int number_of_found_ubvs = 0;
	fprintf(out,"%d, deepest level: %d\n",m_n2,m_max_n2);
	fprintf(out,"Vectors unbiased to the matrix B:\n");
	vector<int*>	found_ubvs;
	for (m_row[1] = 0; m_row[1] < m_n2; m_row[1]++)
	{
		for (m_row[2] = 0; m_row[2] < m_n2; m_row[2]++)
			for (m_row[3] = 0; m_row[3] < m_n2; m_row[3]++)
				for (m_row[4] = 0; m_row[4] < m_n2; m_row[4]++)
					for (m_row[5] = 0; m_row[5] < m_n2; m_row[5]++)
					{
						if(iterate(m_n2))
						{
							// We found an unbiased vector.
							number_of_found_ubvs++;
							fprintf(out,"%d, %d, %d, %d, %d, %d",m_row[0],m_row[1],m_row[2],m_row[3],m_row[4],m_row[5]);
							fprintf(out,", at the deepest level: %d, %d, %d, %d, %d, %d\n",m_deepest[0],m_deepest[1],m_deepest[2],m_deepest[3],m_deepest[4],m_deepest[5]);
							// Saving it for later.
							int* to_load = new int[6];
							for (int i =0; i < 6; i++)
								to_load[i] = m_row[i];
							found_ubvs.push_back(to_load);
						}

					}
	}
	printf("number of found unbiased vectors to matrix B: %d\n",number_of_found_ubvs);
	fprintf(out,"number of found unbiased vectors to matrix B: %d\n",number_of_found_ubvs);

	int number_of_cs=0;
	bool* good_for_d = new bool[found_ubvs.size()];
	// We use the following variable only for statistics.
	int cases_when_at_least_six = 0;
	// searching for 6 vectors from found_ubvs for matrix C
	// All of these vectors are from array found_ubvs and they are
	// orthogonal to each other.
	for (int c0 =  0; c0 < found_ubvs.size(); c0++)
	{
		printf("c0 = %d         \r",c0);
		cout.flush();
		for (int c1 = c0; c1 < found_ubvs.size(); c1++)
		{
			if (scalar_with_eps(found_ubvs[c0],found_ubvs[c1]) == false)
				continue;
			for (int c2 = c1; c2 < found_ubvs.size(); c2++)
			{
				if (scalar_with_eps(found_ubvs[c0],found_ubvs[c2]) == false || scalar_with_eps(found_ubvs[c1],found_ubvs[c2]) == false )
					continue;
				for (int c3 = c2; c3 < found_ubvs.size(); c3++)
				{
					if (scalar_with_eps(found_ubvs[c0],found_ubvs[c3]) == false || scalar_with_eps(found_ubvs[c1],found_ubvs[c3]) == false || scalar_with_eps(found_ubvs[c2],found_ubvs[c3]) == false )
						continue;
					for (int c4 = c3; c4 < found_ubvs.size(); c4++)
					{
						if (scalar_with_eps(found_ubvs[c0],found_ubvs[c4]) == false || scalar_with_eps(found_ubvs[c1],found_ubvs[c4]) == false|| scalar_with_eps(found_ubvs[c2],found_ubvs[c4]) == false || scalar_with_eps(found_ubvs[c3],found_ubvs[c4]) == false )
							continue;
						for (int c5 = c4; c5 < found_ubvs.size(); c5++)
						{
							if (scalar_with_eps(found_ubvs[c0],found_ubvs[c5]) == false || scalar_with_eps(found_ubvs[c1],found_ubvs[c5]) == false|| scalar_with_eps(found_ubvs[c2],found_ubvs[c5]) == false || scalar_with_eps(found_ubvs[c3],found_ubvs[c5]) == false|| scalar_with_eps(found_ubvs[c4],found_ubvs[c5]) == false )
								continue;
							number_of_cs++; //We found a suitable matrix for C
							int k = 0;
							for (int i = 0; i < found_ubvs.size(); i++)
							{
								// We start to find those vectors from found_ubvs which are also
								// unbiased to every row of the matrix C 
								// (we are searching for rows for the matrix D).
								// We copy the matrix C to the array m.
								int m[36];
								for (int j = 0; j < 6; j++)
								{
									m[j] = found_ubvs[c0][j];
									m[j+6] = found_ubvs[c1][j];
									m[j+12] = found_ubvs[c2][j];
									m[j+18] = found_ubvs[c3][j];
									m[j+24] = found_ubvs[c4][j];
									m[j+30] = found_ubvs[c5][j];
								}
								if( ubv_with_matrix(&m[0],found_ubvs[i],m_n2))
								{
									// It would be good for D.
									good_for_d[i] = true;
									k++;
								}
								else
									good_for_d[i] = false;
							}
							// This part is for testing
						//	printf("size of good_for_d = %d, ",k);
						//	if (k > 1)
						//	{
						//		printf("Maybe interesting, the C matrix is:\n");
						//		print(found_ubvs[c0]);
						//		print(found_ubvs[c1]);
						//		print(found_ubvs[c2]);
						//		print(found_ubvs[c3]);
						//		print(found_ubvs[c4]);
						//		print(found_ubvs[c5]);
						//		printf("----- and vectors for D are: \n");
						//		for (int i =0; i < found_ubvs.size(); i++)
						//			if (good_for_d[i])
						//				print(found_ubvs[i]);
						//	}
							// We want to search for a matrix D. If k <=5 then we can not
							// choose 6 different vectors from it.
							if (k > 5)
							{
								cases_when_at_least_six++;
								for (int d0 = 0; d0 < found_ubvs.size(); d0++)
								{
									if (good_for_d[d0] == false)
										continue;
									for (int d1 = d0; d1 < found_ubvs.size(); d1++)
									{
										if (good_for_d[d1] == false || scalar_with_eps(found_ubvs[d0],found_ubvs[d1]) == false)
											continue;
										for (int d2 = d1; d2 < found_ubvs.size(); d2++)
										{
											if (good_for_d[d2] == false || scalar_with_eps(found_ubvs[d0],found_ubvs[d2]) == false || scalar_with_eps(found_ubvs[d1],found_ubvs[d2]) == false)
												continue;
											for (int d3 = d2; d3 < found_ubvs.size(); d3++)
											{
												if (good_for_d[d3] == false || scalar_with_eps(found_ubvs[d0],found_ubvs[d3]) == false || scalar_with_eps(found_ubvs[d1],found_ubvs[d3]) == false || scalar_with_eps(found_ubvs[d2],found_ubvs[d3]) == false)
													continue;
												for (int d4 = d3; d4 < found_ubvs.size(); d4++)
												{
													if (good_for_d[d4] == false || scalar_with_eps(found_ubvs[d0],found_ubvs[d4]) == false || scalar_with_eps(found_ubvs[d1],found_ubvs[d4]) == false || scalar_with_eps(found_ubvs[d2],found_ubvs[d4]) == false || scalar_with_eps(found_ubvs[d3],found_ubvs[d4]) == false)
														continue;
													for (int d5 = d4; d5 < found_ubvs.size(); d5++)
													{
														if (good_for_d[d5] == false || scalar_with_eps(found_ubvs[d0],found_ubvs[d5]) == false|| scalar_with_eps(found_ubvs[d1],found_ubvs[d5]) == false || scalar_with_eps(found_ubvs[d2],found_ubvs[d5]) == false || scalar_with_eps(found_ubvs[d3],found_ubvs[d5]) == false || scalar_with_eps(found_ubvs[d4],found_ubvs[d5]) == false)
															continue;
														printf("ERROR, THIS SHOULD NOT HAPPEN. \n --- \n");
														print(found_ubvs[c0]);
														print(found_ubvs[c1]);
														print(found_ubvs[c2]);
														print(found_ubvs[c3]);
														print(found_ubvs[c4]);
														print(found_ubvs[c5]);
														printf(" and \n");
														print(found_ubvs[d0]);
														print(found_ubvs[d1]);
														print(found_ubvs[d2]);
														print(found_ubvs[d3]);
														print(found_ubvs[d4]);
														print(found_ubvs[d5]);
														printf(" --- \n");
														fprintf(out,"ERROR, THIS SHOULD NOT HAPPEN. \n --- \n");
														fclose(out);
														exit(1);	//We stop the program and exit.
													}
												}
											}
										}
									}
								}
							}
						}
					}
				}
			}
		}
	}
	delete[] good_for_d;
	printf("number of C matrices: %d\n",number_of_cs);
	fprintf(out,"number of C matrices: %d\n",number_of_cs);
	printf("number of cases when we had to search for matrix D: %d\n",cases_when_at_least_six);
	fprintf(out,"number of cases when we had to search for matrix D: %d\n",cases_when_at_least_six);

	if (m_error_with_unsure)
		delete[] m_error_with_unsure;
	if (m_error_with_sure)
		delete[] m_error_with_sure;
	for (int i = 0; i < found_ubvs.size(); i++)
		if (found_ubvs[i])
			delete[] found_ubvs[i];
	found_ubvs.clear();
	
}
bool	Mub::iterate(int n)
{
	// If scalar(n) returns false then
	// the vectors defined by m_row can not be
	// unbiased to the matrix B.
	if (!scalar(n))
		return false;

	// Otherwise if n < m_max_n2 then
	// we double the value of n, and this
	// way we half the intervals of the coordinates.
	if (n < m_max_n2)
	{
		int r1=m_row[1]*2;
		int r2=m_row[2]*2;
		int r3=m_row[3]*2;
		int r4=m_row[4]*2;
		int r5=m_row[5]*2;
	

		bool ret = false;

		for (int d1 = 0; d1 < 2; d1++)
		for (int d2 = 0; d2 < 2; d2++)
		for (int d3 = 0; d3 < 2; d3++)
		for (int d4 = 0; d4 < 2; d4++)
		for (int d5 = 0; d5 < 2; d5++)
		{
			m_row[1]=r1+d1;
			m_row[2]=r2+d2;
			m_row[3]=r3+d3;
			m_row[4]=r4+d4;
			m_row[5]=r5+d5;
			
			if (iterate(2*n))
			{
				ret = true;
				d1 = 2; d2 = 2; d3 = 2; d4 = 2; d5 = 2;
				break;
			}
		}

		m_row[1]=r1/2;
		m_row[2]=r2/2;
		m_row[3]=r3/2;
		m_row[4]=r4/2;
		m_row[5]=r5/2;
		return ret;
	
	}
	else
	{
		for (int i= 0; i<6;i++)
			m_deepest[i] = m_row[i];	//for testing we save the vector from the deepest level
		return true;
	}
}

bool	Mub::ubv_with_matrix(int* m, int* v, int n)
{
	if (n > m_max_n2)
		return true;

	int mult = m_max_n2/n;
	int multn2 = m_max_n2/m_n2;
	double err = (double)5 * PI2 / (double) (2*m_n2) + (double)5 * PI2 / (double) (2*n);

	double x,y,r;
	// This is a basic error estimate. We choose the middle point of 
	// all intervals.
	for (int j = 0; j < 6; j++)
	{
		x = 1;
		y = 0;
		for ( int i = 1; i < 6; i++)
		{
			x += ccos[ (2*v[i]+1)*(mult) ]*ccos[ (2*m[6*j+i]+1)*multn2] + ssin[ (2*v[i]+1)*(mult) ]*ssin[ (2*m[6*j+i]+1)*multn2 ];
			x += ssin[ (2*v[i]+1)*(mult) ]*ccos[ (2*m[6*j+i]+1)*multn2] - ccos[ (2*v[i]+1)*(mult) ]*ssin[ (2*m[6*j+i]+1)*multn2 ];
		}

		if ( (r=sqrt(x*x+y*y)) > m_sqrt6 + err + EPSILON || r < m_sqrt6 - err - EPSILON )
			return false;
	}

	int r1=v[1]*2;
	int r2=v[2]*2;
	int r3=v[3]*2;
	int r4=v[4]*2;
	int r5=v[5]*2;


	bool ret = false;

	for (int d1 = 0; d1 < 2; d1++)
		for (int d2 = 0; d2 < 2; d2++)
			for (int d3 = 0; d3 < 2; d3++)
				for (int d4 = 0; d4 < 2; d4++)
					for (int d5 = 0; d5 < 2; d5++)
					{
						v[1]=r1+d1;
						v[2]=r2+d2;
						v[3]=r3+d3;
						v[4]=r4+d4;
						v[5]=r5+d5;

						if (ubv_with_matrix(m,v,2*n))
						{
							ret = true;
							d1 = 2; d2 = 2; d3 = 2; d4 = 2; d5 = 2;
							break;
						}
					}
	v[1]=r1/2;
	v[2]=r2/2;
	v[3]=r3/2;
	v[4]=r4/2;
	v[5]=r5/2;
	return ret;
}

bool	Mub::scalar(int n)
{
	double x = 1;
	double y = 0;
	int mult = m_max_n2/n;
	double r;
	// This is a basic error estimate. We choose the middle point of 
	// all intervals.
	for ( int i = 1; i < 6; i++)
	{
		x += ccos[ (2*m_row[i]+1)*(mult) ];
		y += ssin[ (2*m_row[i]+1)*(mult) ];
	}
	if ( (r=sqrt(x*x+y*y)) > m_sqrt6 + m_error_with_sure[n] + EPSILON || r < m_sqrt6 - m_error_with_sure[n] - EPSILON )
	{
		return false;
	}


	bool possible = false;
	double xx,yy;
	
	// With our parameters the error estimate is always less than 1. 
	// We make a more sophisticated error estimate, for more details please read
	// our article.
	// (Briefly: Consider one of the ends of the intervals (in all possible ways)
	// and the corresponding scalar product. If it is further to x,y than r,
	// then this vector might be unbiased to (0,0,0,0,0,0).)
	if (r-m_sqrt6 < 0)
		r = m_sqrt6 -r;
	else
		r = r-m_sqrt6;

	for (int i = 0; i < 160; i+=5)
	{
		xx=1.0+ccos[ 2*(m_row[1]+m_shift[i])*mult ]+
			ccos[ 2*(m_row[2]+m_shift[i+1])*mult ]+
			ccos[ 2*(m_row[3]+m_shift[i+2])*mult ]+
			ccos[ 2*(m_row[4]+m_shift[i+3])*mult ]+
			ccos[ 2*(m_row[5]+m_shift[i+4])*mult ] - x;
		yy=ssin[ 2*(m_row[1]+m_shift[i])*mult ]+
			ssin[ 2*(m_row[2]+m_shift[i+1])*mult ]+
			ssin[ 2*(m_row[3]+m_shift[i+2])*mult ]+
			ssin[ 2*(m_row[4]+m_shift[i+3])*mult ]+
			ssin[ 2*(m_row[5]+m_shift[i+4])*mult ] - y;
		if (sqrt(xx*xx+yy*yy) > r-EPSILON)
		{
			possible = true;
			break;
		}
	}
	if(!possible)
		return false;
	// Check whether m_row is unbiased the rows of the matrix B.
	double qqx[2*5];
	double qqy[2*5];
	for (int j = 0; j < 5; j++)
	{
		x = 1;
		y = 0;
		for ( int i = 1; i < 6; i++)
		{
			x += ccos[ (2*m_row[i]+1)*(mult) ]*hcos[3*(6*j+i)+1] + ssin[ (2*m_row[i]+1)*(mult) ]*hsin[3*(6*j+i)+1];
			y += ssin[ (2*m_row[i]+1)*(mult) ]*hcos[3*(6*j+i)+1] - ccos[ (2*m_row[i]+1)*(mult) ]*hsin[3*(6*j+i)+1];
		}
		if (m_fix[j])
		{
			if ( (r=sqrt(x*x+y*y)) > m_sqrt6 + m_error_with_sure[n] + EPSILON || r < m_sqrt6 - m_error_with_sure[n] - EPSILON )
				return false;
		}
		else
		{
			if ( (r=sqrt(x*x+y*y)) > m_sqrt6 + m_error_with_unsure[n] + EPSILON || r < m_sqrt6 - m_error_with_unsure[n] - EPSILON )
				return false;
		}
		//We make the sophisticated error estimate for this row of B as well.
		if (r-m_sqrt6 < 0)
			r = m_sqrt6 -r;
		else
			r = r-m_sqrt6;

		for (int q = 1; q < 6; q++)
		{
			qqx[(q-1)*2+0] = ccos[ (2*m_row[q]+0)*mult ]*hcos[ 18*j+3*q + 2] + ssin[ (2*m_row[q]+0)*mult ]*hsin[ 18*(j)+3*q+2] ;
			qqy[(q-1)*2+0] = ssin[ (2*m_row[q]+0)*mult ]*hcos[ 18*j+3*q + 2] - ccos[ (2*m_row[q]+0)*mult ]*hsin[ 18*(j)+3*q+2] ;
			qqx[(q-1)*2+1] = ccos[ (2*m_row[q]+2)*mult ]*hcos[ 18*j+3*q + 0] + ssin[ (2*m_row[q]+2)*mult ]*hsin[ 18*(j)+3*q+0] ;
			qqy[(q-1)*2+1] = ssin[ (2*m_row[q]+2)*mult ]*hcos[ 18*j+3*q + 0] - ccos[ (2*m_row[q]+2)*mult ]*hsin[ 18*(j)+3*q+0] ;
		}
		possible = false;
		for (int i = 0; i < 160; i+=5)
		{ 
			xx=1.0+qqx[ 0+m_shift[i] ]+
				qqx[ 2+m_shift[i+1] ]+
				qqx[ 4+m_shift[i+2] ]+
				qqx[ 6+m_shift[i+3] ]+
				qqx[ 8+m_shift[i+4] ] - x;
			yy=qqy[ 0+m_shift[i] ]+
				qqy[ 2+m_shift[i+1] ]+
				qqy[ 4+m_shift[i+2] ]+
				qqy[ 6+m_shift[i+3] ]+
				qqy[ 8+m_shift[i+4] ] - y;
			if (sqrt(xx*xx+yy*yy) > r-EPSILON)
			{
				possible = true;
				break;
			}
		}
		if (!possible)
			return false;
	}

	return true;
}



int	main(int argc, const char* argv[])
{
	if (argc != 1) // This program should be called without any parameters.
	{
		printf("Error, usage: fab_ubv\nWithout any parameters.");
		return 1;
	}

	Mub mub(19,7); // n2=19, deep=7

	FILE*	out;
	out = fopen("fab_ubv.txt","w");
	
	printf("START\n");

	int n1 = 180; 

	// It is enough to check a=0,1,2,...,29
	// and b=0,1,2,...,a/2+1 cases
	for (int a=0; a < 30; a++)
		for (int b=0; b <= a/2+1; b++)
		{
			printf("a=%d, b=%d\n",a,b);
			fprintf(out,"a=%d, b=%d\n",a,b);
			mub.start(n1,a,b,out);
			fprintf(out,"\n");
		}
	fclose(out);

	return 0;
}