#include <typeinfo>
#include <iostream>
#include <fstream>
#include <cstdio>
#include <algorithm>
#include <stack>
#include <queue>
#include <vector>
#include <map>
#include <set>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <ctime>

#include <mie/gmp_util.hpp>
#include <mie/fp.hpp>
#include "bn.h"

using namespace bn;
using namespace std;

struct TagFr; // any name to distinguish Fr from other finite classes

typedef mie::FpT<mie::Gmp, TagFr> Fr;

typedef pair<int,int> pii;

#define DEBUG_MODE 0
#define LEVEL 11

#define NUM_GATE 2047
#define NUM_WIRE 2046
#define NUM_WITNESS 6141

#define NAND_GATE 0
#define FORK_GATE 1

/*
  MSB first binary method. Multi exponentiation
*/
Fp ***gTable;
Fp2 ***gTable2;
mie::Vuint *Vv;
int mSize[(NUM_GATE + NUM_WIRE) * 6];
template<class FF, class INT>
inline void ScalarMultiMult(FF *out, FF (*in)[3], INT *m,
	const int &SIZE) {
	typedef typename INT::value_type value_type;

	bool allZero = true;
	for (int i = 0; i < SIZE; ++i) {
		if (m[i] != 0) {
			allZero = false;
			break;
		}
	}
	if (allZero) {
    out[2] = 0;
    return;
  }

	assert(SIZE % 4 == 0);

	// precompute pair of bases (sliding window technique)
	FF ***table;
	//table = new FF*[SIZE>>1];
	if (typeid(table).name() == typeid(gTable).name()) {
		table = (FF***)gTable;
	} else {
		table = (FF***)gTable2;
	}
	for (int i = 0; i < SIZE; i += 4) {
		//table[i>>1] = new FF[3];
		int tmpidx = i>>2;
		ecop::copy(table[tmpidx][1], in[i]);
		ecop::copy(table[tmpidx][2], in[i+1]);
		ecop::copy(table[tmpidx][4], in[i+2]);
		ecop::copy(table[tmpidx][8], in[i+3]);
		ecop::ECAdd(table[tmpidx][3], in[i], in[i+1]);
		ecop::ECAdd(table[tmpidx][5], in[i], in[i+2]);
		ecop::ECAdd(table[tmpidx][6], in[i+1], in[i+2]);
		ecop::ECAdd(table[tmpidx][9], in[i], in[i+3]);
		ecop::ECAdd(table[tmpidx][10], in[i+1], in[i+3]);
		ecop::ECAdd(table[tmpidx][12], in[i+2], in[i+3]);
		ecop::ECAdd(table[tmpidx][7], in[i], table[tmpidx][6]);
		ecop::ECAdd(table[tmpidx][11], in[i+3], table[tmpidx][3]);
		ecop::ECAdd(table[tmpidx][13], in[i], table[tmpidx][12]);
		ecop::ECAdd(table[tmpidx][14], in[i+2], table[tmpidx][10]);
		ecop::ECAdd(table[tmpidx][15], in[i], table[tmpidx][14]);
	}

	//int *mSize = new int[SIZE];
	int mSizeMax = -1;
	for (int i = 0; i < SIZE; ++i) {
		mSize[i] = (int)m[i].size();
		mSizeMax = max(mSizeMax, mSize[i] - 1);
	}

  const int vSize = (int)sizeof(value_type) * 8;
  const value_type mask = value_type(1) << (vSize - 1);

  /*
    Extract and process for MSB of most significant word.
  */
	FF temp[3];
	FF *temp2;
	value_type *v = (value_type*)Vv;//new value_type[SIZE];
	for (int i = 0; i < SIZE; ++i) {
		v[i] = (mSize[i] - 1) == mSizeMax ? m[i][mSize[i]-1] : 0;
	}
	bool foundMSB = false;
	int j = 0;
	while (!foundMSB) {
		for (int i = 0; i < SIZE; i += 4) {
			int x = ((v[i+3] & mask) | ((v[i+2] & mask) >> 1) |
				((v[i+1] & mask) >> 2) |
				((v[i] & mask) >> 3)) >> (vSize-4);
			if (x) {
				//temp2 = x == 1 ? in[i] : (x == 2 ? in[i+1] : table[i>>1]);
				temp2 = table[i>>2][x];
				if (foundMSB) {
					ecop::copy(temp, out);
					ecop::ECAdd(out, temp, temp2);
				} else {
					ecop::copy(out, temp2);
				}
				foundMSB = true;
			}
			v[i] <<= 1;
			v[i+1] <<= 1;
			v[i+2] <<= 1;
			v[i+3] <<= 1;
		}
		++j;
	}

  /*
    Process for most significant word.
  */
  for (; j != vSize; ++j) {
		ecop::copy(temp, out);
		ecop::ECDouble(out, temp);
		for (int i = 0; i < SIZE; i += 4) {
			int x = ((v[i+3] & mask) | ((v[i+2] & mask) >> 1) |
				((v[i+1] & mask) >> 2) |
				((v[i] & mask) >> 3)) >> (vSize-4);
			if (x) {
				//temp2 = x == 1 ? in[i] : (x == 2 ? in[i+1] : table[i>>1]);
				temp2 = table[i>>2][x];
				ecop::copy(temp, out);
				ecop::ECAdd(out, temp, temp2);
			}
			v[i] <<= 1;
			v[i+1] <<= 1;
			v[i+2] <<= 1;
			v[i+3] <<= 1;
		}
  }

  /*
    Process for non most significant words.
  */
  for (int i = mSizeMax - 1; i >= 0; --i) {
		for (j = 0; j < SIZE; ++j) {
			v[j] = i < mSize[j] ? m[j][i] : 0;
		}
    for (j = 0; j != vSize; ++j) {
			ecop::copy(temp, out);
			ecop::ECDouble(out, temp);
			for (int k = 0; k < SIZE; k += 4) {
				int x = ((v[k+3] & mask) | ((v[k+2] & mask) >> 1) |
				((v[k+1] & mask) >> 2) |
				((v[k] & mask) >> 3)) >> (vSize-4);
				if (x) {
					//temp2 = x == 1 ? in[k] : (x == 2 ? in[k+1] : table[k>>1]);
					temp2 = table[k>>2][x];
					ecop::copy(temp, out);
					ecop::ECAdd(out, temp, temp2);
				}
				v[k] <<= 1;
				v[k+1] <<= 1;
				v[k+2] <<= 1;
				v[k+3] <<= 1;
			}
    }
  }
}

/**
 *  Gate specification: gate_type, input1, input2, output
 *  all are variable number, e.g. x_1, x_2, ..., x_n
 */
int gateSpec[NUM_GATE][4]; /* = {
	{NAND_GATE, 1, 2, 5, -1},
	{NAND_GATE, 3, 4, 6, -1},
	{NAND_GATE, 7, 8, 9, -1}
};*/

/**
 *  Wire specification: var1, var2, var3
 *  all are variable number, e.g. x_1, x_2, ..., x_n
 *  shows which of the variables are in the same wire
 *  -1 means fanout is < 2, -1 should only happen on the end
 */
int wireSpec[NUM_WIRE][3]; /* = {
	{6, 8, -1},
	{5, 7, -1}
};*/

/**
 *  Reed Solomon generator for fanout 1 and 2
 */
int RSGeneratorF1[2][3] = {
	{1, 0, -1},
	{0, 1, 2}
};

int RSGeneratorF2[3][5] = {
	{1, 0, 0, 1, 3},
	{0, 1, 0, -3, -8},
	{0, 0, 1, 3, 6}
};

/**
 * Witness vector should be given by the order of gate appearance
 */

int witnessA[NUM_WITNESS][2]; /* = {
	{1, 1}, {2, 1}, {-5, 1},
	{3, 1}, {4, 1}, {-6, 1},
	{-7, 2}, {-8, -1}, {9, 1}
};*/

int witnessB[NUM_WITNESS][2]; /* = {
	{1, 2}, {2, -1}, {-5, 1},
	{3, 3}, {4, -2}, {-6, 1},
	{-7, 1}, {-8, 1}, {9, 1}
};*/

int varNum = 1, gateNum = 0, witnessNum = 0, wireNum = 0;
inline void genCircuit(int level) {
	if (level == LEVEL) {
		assert(gateNum == NUM_GATE);
		assert(witnessNum == NUM_WITNESS);
		assert(wireNum == NUM_WIRE);
		return;
	}
	int par = level == 0 ? -1 : varNum - (1 << (level - 1)) * 3;
	for (int i = 0; i < (1 << level); i++, gateNum++) {
		gateSpec[gateNum][0] = NAND_GATE;
		for (int j = 1; j < 4; j++, varNum++, witnessNum++) {
			gateSpec[gateNum][j] = varNum;
			int val;
			if (level & 1) {
				val = j == 3 ? -(witnessNum+1) : witnessNum+1;
				witnessA[witnessNum][1] = 1;
				switch (j) {
					case 1:
						witnessB[witnessNum][1] = 2;
						break;
					case 2:
						witnessB[witnessNum][1] = -1;
						break;
					case 3:
						witnessB[witnessNum][1] = 1;
						break;
				}
			} else {
				val = j == 3 ? witnessNum+1 : -(witnessNum+1);
				switch (j) {
					case 1:
						witnessA[witnessNum][1] = 2;
						break;
					case 2:
						witnessA[witnessNum][1] = -1;
						break;
					case 3:
						witnessA[witnessNum][1] = 1;
						break;
				}
				witnessB[witnessNum][1] = 1;
			}
			witnessA[witnessNum][0] = val;
			witnessB[witnessNum][0] = val;
			if (j == 3 && level != 0) {
				assert(varNum % 3 == 0);
				wireSpec[wireNum][0] = varNum;
				wireSpec[wireNum][1] = par;
				wireSpec[wireNum][2] = -1;
				wireNum++;
				par++;
				if (par % 3 == 0) {
					par++;
				}
			}
		}
	}
	genCircuit(level+1);
}

vector<int> widthG, heightG;
int totalWidthG = 0, totalHeightG = 0;
vector<int> aggregatedGateCheckerTarget;
// map of labelling to vectors of column number where value is 1
vector<map<int, vector<int> > > aggregatedGateCheckerSP;

inline void aggregateGateChecker() {
	// for each gate, generate the span program of that gate
	// after that, add that program to the aggregate gate checker
	for (int i = 0; i < NUM_GATE; i++) {
		map<int, vector<int> > spanProgram;
		switch (gateSpec[i][0]) {
			case NAND_GATE:
				widthG.push_back(3);
				heightG.push_back(6);
				totalWidthG += 3;
				totalHeightG += 6;
				aggregatedGateCheckerTarget.push_back(1);
				aggregatedGateCheckerTarget.push_back(1);
				aggregatedGateCheckerTarget.push_back(1);
				// assign labelling, > 0 is positive literal, < 0 otherwise
				for (int j = 0; j < 3; j++) {
					spanProgram[gateSpec[i][j+1]] = vector<int>();
					spanProgram[-gateSpec[i][j+1]] = vector<int>();
				}
				// construct span program
				spanProgram[gateSpec[i][1]].push_back(0);
				spanProgram[gateSpec[i][2]].push_back(1);
				spanProgram[gateSpec[i][3]].push_back(0);
				spanProgram[gateSpec[i][3]].push_back(1);
				spanProgram[-gateSpec[i][1]].push_back(2);
				spanProgram[-gateSpec[i][2]].push_back(2);
				spanProgram[-gateSpec[i][3]].push_back(2);
				aggregatedGateCheckerSP.push_back(spanProgram);
			break;
			case FORK_GATE:
				widthG.push_back(3);
				heightG.push_back(6);
				totalWidthG += 3;
				totalHeightG += 6;
				aggregatedGateCheckerTarget.push_back(1);
				aggregatedGateCheckerTarget.push_back(1);
				aggregatedGateCheckerTarget.push_back(1);
				// assign labelling, > 0 is positive literal, < 0 otherwise
				for (int j = 0; j < 3; j++) {
					spanProgram[gateSpec[i][j+1]] = vector<int>();
					spanProgram[-gateSpec[i][j+1]] = vector<int>();
				}
				// construct span program
				spanProgram[gateSpec[i][1]].push_back(1);
				spanProgram[gateSpec[i][2]].push_back(2);
				spanProgram[gateSpec[i][3]].push_back(0);
				spanProgram[-gateSpec[i][1]].push_back(0);
				spanProgram[-gateSpec[i][2]].push_back(1);
				spanProgram[-gateSpec[i][3]].push_back(2);
				aggregatedGateCheckerSP.push_back(spanProgram);
			break;
		}
	}

	if (DEBUG_MODE) {
		cout << "Aggregated Gate checker" << endl;
		cout << "Target: ";
		for (int i = 0; i < aggregatedGateCheckerTarget.size(); i++) {
			cout << aggregatedGateCheckerTarget[i] << " ";
		}
		cout << endl;
		for (int i = 0; i < aggregatedGateCheckerSP.size(); i++) {
			map<int, vector<int> >::iterator iter = aggregatedGateCheckerSP[i].begin();
			cout << "Gate " << i << endl;
			for ( ; iter != aggregatedGateCheckerSP[i].end(); iter++) {
				cout << "Var: " << iter->first << " === ";
				vector<int> position = iter->second;
				for (int k = 0; k < position.size(); k++) {
					cout << position[k] << " ";
				}
				cout << endl;
			}
			cout << endl;
		}
	}
}

vector<int> cumWidthW, cumHeightW;
int totalWidthW = 0, totalHeightW = 0;
map<int, int> wireCheckerMap; //map of variable to wire checker
vector<map<int, vector<pii> > > aggregatedWireCheckerQSP;

inline void aggregateWireChecker() {
	// for each wire, construct the wire checker
	for (int i = 0; i < NUM_WIRE; i++) {
		if (wireSpec[i][2] == -1) { //fanout 1
			cumWidthW.push_back(totalWidthW);
			cumHeightW.push_back(totalHeightW);
			totalWidthW += 6;
			totalHeightW += 4;
			// row 1 & 2 is positive, row 3 & 4 is negative
			map<int, vector<pii> > spanProgram;
			// assign labelling
			for (int j = 0; j < 2; j++) {
				spanProgram[wireSpec[i][j]] = vector<pii>();
				spanProgram[-wireSpec[i][j]] = vector<pii>();
				wireCheckerMap[wireSpec[i][j]] = i;
				wireCheckerMap[-wireSpec[i][j]] = i;
			}
			// construct the QSP
			for (int j = 0; j < 2; j++) {
				for (int k = 0; k < 3; k++) {
					if (RSGeneratorF1[j][k]) {
						spanProgram[wireSpec[i][j]].push_back(
							make_pair(k, RSGeneratorF1[j][k]));
						spanProgram[-wireSpec[i][j]].push_back(
							make_pair(3 + k, RSGeneratorF1[j][k]));
					}
				}
			}
			aggregatedWireCheckerQSP.push_back(spanProgram);
		} else { //fanout 2
			cumWidthW.push_back(totalWidthW);
			cumHeightW.push_back(totalHeightW);
			totalWidthW += 10;
			totalHeightW += 6;
			// row 1 - 3 is positive, row 4 - 6 is negative
			map<int, vector<pii> > spanProgram;
			// assign labelling
			for (int j = 0; j < 3; j++) {
				spanProgram[wireSpec[i][j]] = vector<pii>();
				spanProgram[-wireSpec[i][j]] = vector<pii>();
				wireCheckerMap[wireSpec[i][j]] = i;
				wireCheckerMap[-wireSpec[i][j]] = i;
			}
			// construct the QSP
			for (int j = 0; j < 3; j++) {
				for (int k = 0; k < 5; k++) {
					if (RSGeneratorF2[j][k]) {
						spanProgram[wireSpec[i][j]].push_back(
							make_pair(k, RSGeneratorF2[j][k]));
						spanProgram[-wireSpec[i][j]].push_back(
							make_pair(5 + k, RSGeneratorF2[j][k]));
					}
				}
			}
			aggregatedWireCheckerQSP.push_back(spanProgram);
		}
	}

	if (DEBUG_MODE) {
		cout << "Aggregated Wire checker" << endl;
		for (int i = 0; i < aggregatedWireCheckerQSP.size(); i++) {
			map<int, vector<pii> >::iterator iter = aggregatedWireCheckerQSP[i].begin();
			cout << "Wire " << i << endl;
			for (; iter != aggregatedWireCheckerQSP[i].end(); iter++) {
				cout << "Var: " << iter->first << " === ";
				vector<pii> position = iter->second;
				for (int k = 0; k < position.size(); k++) {
					cout << "(" << position[k].first << "," << position[k].second << ")" << " ";
				}
				cout << endl;
			}
			cout << endl;
		}
	}
}

inline void verifyCircuitChecker() {
	// verify the gate checker V and W part
	int sumVg[totalWidthG];
	int sumWg[totalWidthG];
	memset(sumVg, 0, sizeof(sumVg));
	memset(sumWg, 0, sizeof(sumWg));
	int s = 0;
	int currentWidth = 0;
	map<int, vector<int> > gateChecker = aggregatedGateCheckerSP[s];
	for (int i = 0; i < NUM_WITNESS; i++) {
		if (gateChecker.find(witnessA[i][0]) ==
			gateChecker.end()) {
			currentWidth += widthG[s];
			gateChecker = aggregatedGateCheckerSP[++s];
		}
		if (witnessA[i][1]) {
			vector<int> position = gateChecker[witnessA[i][0]];
			assert(position.size() > 0);
			for (int j = 0; j < position.size(); j++) {
				sumVg[currentWidth + position[j]] += witnessA[i][1];
			}
		}
		if (witnessB[i][1]) { // Notice: witnessA[i][0] = witnessB[i][0]
			vector<int> position = gateChecker[-witnessB[i][0]];
			assert(position.size() > 0);
			for (int j = 0; j < position.size(); j++) {
				sumWg[currentWidth + position[j]] += witnessB[i][1];
			}
		}
	}
	for (int i = 0; i < totalWidthG; i++) {
		assert(sumVg[i] == aggregatedGateCheckerTarget[i]);
		assert(sumWg[i] == aggregatedGateCheckerTarget[i]);
	}
	cout << "Gate checker Vg, Wg verification successful!" << endl;

	// verify the wire checker V and W part
	int sumVw[totalWidthW];
	int sumWw[totalWidthW];
	memset(sumVw, 0, sizeof(sumVw));
	memset(sumWw, 0, sizeof(sumWw));
	for (int i = 0; i < NUM_WITNESS; i++) {
		if (witnessA[i][1] &&
			wireCheckerMap.find(witnessA[i][0]) != wireCheckerMap.end()) {
			int s = wireCheckerMap[witnessA[i][0]];
			map<int, vector<pii> > wireChecker =
				aggregatedWireCheckerQSP[s];
			vector<pii> position = wireChecker[witnessA[i][0]];
			assert(position.size() > 0);
			for (int j = 0; j < position.size(); j++) {
				sumVw[cumWidthW[s] + position[j].first] +=
					witnessA[i][1] * position[j].second;
			}
		}
		if (witnessB[i][1] &&
			wireCheckerMap.find(-witnessB[i][0]) != wireCheckerMap.end()) {
			int s = wireCheckerMap[-witnessB[i][0]];
			map<int, vector<pii> > wireChecker =
				aggregatedWireCheckerQSP[s];
			vector<pii> position = wireChecker[-witnessB[i][0]];
			assert(position.size() > 0);
			for (int j = 0; j < position.size(); j++) {
				sumWw[cumWidthW[s] + position[j].first] +=
					witnessB[i][1] * position[j].second;
			}
		}
	}
	for (int i = 0; i < totalWidthW; i++) {
		assert(sumVw[i] == 0 || sumWw[i] == 0);
	}
	cout << "Wire checker Vw, Ww verification successful!" << endl;
}

/**
 *  CRS part
 *  VLabelling, WLabelling: labelling of the variable (0 is target)
 *  circuitV: the V matrix of circuit checker
 *  circuitW: the W matrix
 *  first row of each matrix is the -target (for v0 and w0)
 */
vector<int> VLabelling, WLabelling;
vector<vector<int> > circuitV, circuitW;

inline int generateCircuitVW(char writeToFile) {
	circuitV.resize(1, vector<int>(2 * totalWidthG + totalWidthW, 0));
	circuitW.resize(1, vector<int>(2 * totalWidthG + totalWidthW, 0));
	VLabelling.push_back(0);
	WLabelling.push_back(0);
	// for outer program usage
	ofstream myfile;
	if (writeToFile == '0') {
		myfile.open("circuit.txt");
		myfile << Param::r << endl;
		myfile << 1 << endl;
		myfile << "[ ";
	}
	/////////////////////

	for (int i = 0; i < aggregatedGateCheckerTarget.size(); i++) {
		circuitV[0][i] = -aggregatedGateCheckerTarget[i];
		circuitW[0][totalWidthG + i] = -aggregatedGateCheckerTarget[i];
	}
	for (int i = 0; i < totalWidthG; i++) {
		circuitV[0][totalWidthG + i] = -1;
		circuitW[0][i] = -1;
	}

	/////////////////////
	if (writeToFile == '0') {
		for (int j = 0; j < circuitV[0].size(); j++) {
			myfile << circuitV[0][j] << " ";
		}
		myfile << "]" << endl;
		myfile << "[ ";
		for (int j = 0; j < circuitW[0].size(); j++) {
			myfile << circuitW[0][j] << " ";
		}
		myfile << "]" << endl;
	}
	/////////////////////
	int row = 1;
	for (int i = 0; i < aggregatedGateCheckerSP.size(); i++) {
		map<int, vector<int> >::iterator iter = aggregatedGateCheckerSP[i].begin();
		for (; iter != aggregatedGateCheckerSP[i].end(); iter++) {
			VLabelling.push_back(iter->first);
			WLabelling.push_back(-iter->first);
			row++;
		}
	}

	/*
	if (DEBUG_MODE) {
		puts("Matrix V from circuit checker:");
		for (int i = 0; i < circuitV.size(); i++) {
			printf("%3d", VLabelling[i]);
			for (int j = 0; j < circuitV[i].size(); j++) {
				printf("%3d", circuitV[i][j]);
			}
			puts("");
		}
		puts("");
		puts("Matrix W from circuit checker:");
		for (int i = 0; i < circuitV.size(); i++) {
			printf("%3d", WLabelling[i]);
			for (int j = 0; j < circuitV[i].size(); j++) {
				printf("%3d", circuitW[i][j]);
			}
			puts("");
		}
		puts("");
	}
	*/

	/////////////////////
	if (writeToFile == '0') {
		myfile.flush();
		myfile.close();
	}
	/////////////////////

	return row;
}

/**
 *  Prover part
 *  sumV: the sum a_i . v_i (based on the witnessA)
 *  sumW: the sum b_i . w_i (based on the witnessB)
 */

inline vector<vector<int> > calcSumVW() {
	vector<int> sumV;
	vector<int> sumW;
	// make the witness map first
	map<int, int> witnessMapA;
	map<int, int> witnessMapB;
	for (int i = 0; i < NUM_WITNESS; i++) {
		witnessMapA[witnessA[i][0]] = witnessA[i][1];
		witnessMapB[witnessB[i][0]] = witnessB[i][1];
	}
	sumV.resize(2 * totalWidthG + totalWidthW, 0);
	sumW.resize(2 * totalWidthG + totalWidthW, 0);
	int currentWidthG = 0;
	for (int i = 0; i < aggregatedGateCheckerSP.size(); i++) {
		map<int, vector<int> >::iterator iter = aggregatedGateCheckerSP[i].begin();
		for (; iter != aggregatedGateCheckerSP[i].end(); iter++) {
			if (witnessMapA.find(iter->first) != witnessMapA.end()) {
				vector<int> position = iter->second;
				for (int j = 0; j < position.size(); j++) {
					sumV[currentWidthG + position[j]] += witnessMapA[iter->first];
				}
				if (wireCheckerMap.find(iter->first) != wireCheckerMap.end()) {
					int s = wireCheckerMap[iter->first];
					map<int, vector<pii> > wireChecker = aggregatedWireCheckerQSP[s];
					vector<pii> position2 = wireChecker[iter->first];
					for (int k = 0; k < position2.size(); k++) {
						sumV[2 * totalWidthG + cumWidthW[s] + position2[k].first] +=
							witnessMapA[iter->first] * position2[k].second;
					}
				}
			}
			if (witnessMapB.find(-iter->first) != witnessMapB.end()) {
				vector<int> position = iter->second;
				for (int j = 0; j < position.size(); j++) {
					sumW[totalWidthG + currentWidthG + position[j]] += witnessMapB[-iter->first];
				}
				if (wireCheckerMap.find(iter->first) != wireCheckerMap.end()) {
					int s = wireCheckerMap[iter->first];
					map<int, vector<pii> > wireChecker = aggregatedWireCheckerQSP[s];
					vector<pii> position2 = wireChecker[iter->first];
					for (int k = 0; k < position2.size(); k++) {
						sumW[2 * totalWidthG + cumWidthW[s] + position2[k].first] +=
							witnessMapB[-iter->first] * position2[k].second;
					}
				}
			}
		}
		currentWidthG += widthG[i];
	}
	if (DEBUG_MODE) {
		cout << "sumV" << endl;
		for (int i = 0; i < sumV.size(); i++) {
			cout << sumV[i] << " ";
		}
		cout << endl;
		cout << "sumW" << endl;
		for (int i = 0; i < sumW.size(); i++) {
			cout << sumW[i] << " ";
		}
		cout << endl;
	}
	vector<vector<int> > sumVW;
	sumVW.push_back(sumV);
	sumVW.push_back(sumW);
	return sumVW;
}

/**
 * precompute all basis polynomial for Lagrange interpolation
 */

vector<vector<Fr> > basisPoly;
vector<Fr> numerator;
Fr firstDenom;

inline void precomputeInterpolation(int degree) {
	// pre-compute (x-x1) . (x-x2) ... (x-xn)
	vector<Fr> numerator;
	numerator.push_back(Fr(1));
	numerator.push_back(Fr(-1));
	for (int j = 1; j < degree; j++) {
		numerator.push_back(Fr(0)); // times x
		// add with time -value
		Fr prev = Fr(numerator[0]);
		for (int k = 1; k < numerator.size(); k++) {
			Fr temp = Fr(-(j+1)) * prev;
			prev = numerator[k];
			numerator[k] += temp;
		}
	}
	// construct l_j(x)
	for (int j = 0; j < degree; j++) {
		// polynomial long division
		vector<Fr> result;
		vector<Fr> tempNumerator(numerator);
		for (int k = 0; k + 1 < tempNumerator.size(); k++) {
			result.push_back(tempNumerator[k]);
			tempNumerator[k + 1] -= Fr(-(j+1)) * tempNumerator[k];
		}
		assert(tempNumerator[tempNumerator.size() - 1] == Fr(0));
		basisPoly.push_back(result);
	}

	// first denom will be from -1 to -degree+1
	firstDenom = Fr(1);
	for (int i = -1; i >= -degree+1; i--) {
		firstDenom *= Fr(i);
	}
}


/**
 *  r_j = j+1
 */

inline vector<vector<mie::Vuint> > interpolate(vector<vector<int> >& V) {
	vector<vector<mie::Vuint> > Vx;
	int n = V.size();
	int degree = V[0].size();
	assert(degree == basisPoly[0].size());
	for (int i = 0; i < n; i++) {
		vector<Fr> polynomial(degree, Fr(0));
		Fr denom(firstDenom);
		int highest = -1;
		int lowest = -degree+1;
		for (int j = 0; j < degree; j++) {
			Fr a(V[i][j]);
			Fr c = a / denom;
			assert(basisPoly[j].size() == polynomial.size());
			for (int k = 0; k < basisPoly[j].size(); k++) {
				polynomial[k] += c * basisPoly[j][k];
			}
			denom /= Fr(lowest);
			lowest++;
			highest++;
			if (lowest == 0) {
				lowest++;
			}
			if (highest == 0) {
				highest++;
			}
			denom *= Fr(highest);
		}
		// just for testing purpose, comment out for efficiency !!!!!
		if (DEBUG_MODE) {
			for (int j = 0; j < degree; j++) {
				Fr x(j+1);
				Fr calc(0);
				for (int k = 0; k < polynomial.size(); k++) {
					calc *= x;
					calc += polynomial[k];
				}
				assert((calc == Fr(V[i][j])));
			}
		}
		vector<mie::Vuint> dummy;
		for (int i = 0; i < polynomial.size(); i++) {
			dummy.push_back(mie::Vuint(polynomial[i].toStr()));
		}
		Vx.push_back(dummy);
	}
	return Vx;
}


inline mie::Vuint calcZSigma(mie::Vuint& sigma1, int d) {
	Fr result(1);
	Fr sigma(sigma1.toString());
	for (int i = 1; i <= d; i++) {
		result *= (sigma - Fr(i));
	}
	return mie::Vuint(result.toStr());
}

inline vector<mie::Vuint> calcVWSigma(mie::Vuint& sigma1, mie::Vuint& zSigma1, int d) {
	Fr test = -1;
	vector<Fr> L, Lprime, basisPoly;
	Fr sigma(sigma1.toString()), zSigma(zSigma1.toString());
	for (int i = 1; i <= d; i++) {
		L.push_back(zSigma / (sigma - Fr(i)));
	}
	// compute Lprime_0
	Fr denom = 1;
	for (int i = -1; i >= -d + 1; i--) {
		denom *= Fr(i);
	}
	Lprime.push_back(denom);
	int highest = -1;
	int lowest = -d + 1;
	// compute rest of Lprime
	for (int i = 1; i < d; i++) {
		denom /= Fr(lowest);
		lowest++;
		highest++;
		if (lowest == 0) {
			lowest++;
		}
		if (highest == 0) {
			highest++;
		}
		denom *= Fr(highest);
		Lprime.push_back(denom);
	}
	for (int i = 0; i < d; i++) {
		basisPoly.push_back(L[i] / Lprime[i]);
	}

	vector<mie::Vuint> result, result2;
	Fr temp = 0, temp2 = 0;
	for (int j = 0; j < circuitV[0].size(); j++) {
		temp += circuitV[0][j] * basisPoly[j];
	}
	result.push_back(mie::Vuint(temp.toStr()));
	for (int j = 0; j < circuitW[0].size(); j++) {
		temp2 += circuitW[0][j] * basisPoly[j];
	}
	result2.push_back(mie::Vuint(temp2.toStr()));

	int currentWidthG = 0;
	for (int i = 0; i < aggregatedGateCheckerSP.size(); i++) {
		map<int, vector<int> >::iterator iter = aggregatedGateCheckerSP[i].begin();
		for (; iter != aggregatedGateCheckerSP[i].end(); iter++) {
			temp = 0;
			temp2 = 0;
			vector<int> position = iter->second;
			for (int j = 0; j < position.size(); j++) {
				temp += basisPoly[currentWidthG + position[j]];
				temp2 += basisPoly[totalWidthG + currentWidthG + position[j]];
			}
			if (wireCheckerMap.find(iter->first) != wireCheckerMap.end()) {
				int s = wireCheckerMap[iter->first];
				map<int, vector<pii> > wireChecker = aggregatedWireCheckerQSP[s];
				vector<pii> position2 = wireChecker[iter->first];
				for (int k = 0; k < position2.size(); k++) {
					temp += Fr(position2[k].second) *
						basisPoly[2 * totalWidthG + cumWidthW[s] + position2[k].first];
					temp2 += Fr(position2[k].second) *
						basisPoly[2 * totalWidthG + cumWidthW[s] + position2[k].first];
				}
			}
			result.push_back(mie::Vuint(temp.toStr()));
			result2.push_back(mie::Vuint(temp2.toStr()));
		}
		currentWidthG += widthG[i];
	}
	for (int i = 0; i < result2.size(); i++) {
		result.push_back(result2[i]);
	}
	return result;
}

inline vector<mie::Vuint> calcZx(int degree) {
	// compute (x-x1) . (x-x2) ... (x-xn)
	vector<Fr> result;
	result.push_back(Fr(1));
	result.push_back(Fr(-1));
	for (int j = 1; j < degree; j++) {
		result.push_back(Fr(0)); // times x
		// add with time -value
		Fr prev(result[0]);
		for (int k = 1; k < result.size(); k++) {
			Fr temp = Fr(-(j+1)) * prev;
			prev = result[k];
			result[k] += temp;
		}
	}
	if (DEBUG_MODE) {
		for (int j = 1; j <= degree; j++) {
			Fr y(0);
			for (int i = 0; i < result.size(); i++) {
				y *= Fr(j);
				y += result[i];
			}
			assert(y == Fr(0));
		}
	}
	vector<mie::Vuint> result2(result.size());
	for (int i = 0; i < result.size(); i++) {
		result2[i] = mie::Vuint(result[i].toStr());
	}
	return result2;
}

inline vector<mie::Vuint> addPoly(vector<mie::Vuint>& polyA, vector<mie::Vuint>& polyB) {
	assert(polyA.size() == polyB.size());
	vector<mie::Vuint> resultPoly(polyA.size());
	for (int i = 0; i < polyA.size(); i++) {
		resultPoly[i] = polyA[i] + polyB[i];
		resultPoly[i] %= Param::r;
	}
	return resultPoly;
}

inline vector<mie::Vuint> calcHx(vector<mie::Vuint>& VDagger, vector<mie::Vuint>& WDagger, vector<mie::Vuint>& Zx) {
	// calculate VDagger . WDagger
	assert(Zx.size() == VDagger.size() + 1);
	assert(VDagger.size() == WDagger.size());
	vector<mie::Vuint> multResult(VDagger.size() + WDagger.size() - 1, mie::Vuint(0));
	for (int i = 0; i < WDagger.size(); i++) {
		for (int j = 0; j < VDagger.size(); j++) {
			multResult[i + j] += VDagger[j] * WDagger[i];
			multResult[i + j] %= Param::r;
		}
	}
	vector<mie::Vuint> Hx;
	assert(Zx[0] == 1);
	// calculate multResult / Zx
	for (int i = 0; i + Zx.size() - 1 < multResult.size(); i++) {
		mie::Vuint quot(multResult[i]);
		Hx.push_back(quot);
		multResult[i] = 0;
		for (int j = 1; j < Zx.size(); j++) {
			mie::Vsint temp(quot);
			temp *= Zx[j];
			temp *= -1;
			temp %= Param::r;
			if (temp.isNegative()) {
				temp += Param::r;
			}
			multResult[i + j] += temp.abs();
			multResult[i + j] %= Param::r;
		}
	}

	if (DEBUG_MODE) {
		for (int i = 0; i < multResult.size(); i++) {
			assert(multResult[i] == 0);
		}
	}

	return Hx;
}

vector<mie::Vuint> vwSigma;

int main(int argc, char** argv) {
	if (argc < 2) {
		cout << "Select a mode!" << endl;
		return 0;
	}

	genCircuit(0);
	srand(time(NULL));
	Param::init(-1);
	string rStr = Param::r.toString();
	Fr::setModulo(rStr);
	const Fp g[3] = {
		Fp("1674578968009266105367653690721407808692458796109485353026408377634195183292"),
		Fp("8299158460239932124995104248858950945965255982743525836869552923398581964065"),
		Fp("1")
	};
	const Fp2 g2[3] = {
		Fp2(
			Fp("12723517038133731887338407189719511622662176727675373276651903807414909099441"),
			Fp("4168783608814932154536427934509895782246573715297911553964171371032945126671")
		),
		Fp2(
			Fp("13891744915211034074451795021214165905772212241412891944830863846330766296736"),
			Fp("7937318970632701341203597196594272556916396164729705624521405069090520231616")
		),
		Fp2(
			Fp("1"),
			Fp("0")
		)
	};

	aggregateGateChecker();
	aggregateWireChecker();
	int m = generateCircuitVW(argv[1][0]);
	vector<vector<mie::Vuint> > Vx, Wx;

	clock_t begin = clock();

	vector<vector<int> > sumVW = calcSumVW();
	vector<vector<int> > proverMatrix(sumVW);
	vector<vector<mie::Vuint> > proverPoly;

	clock_t end = clock();
	double interpolateTime, CRSInterpolateTime, proverComputeTime;

	int degree = 2 * totalWidthG + totalWidthW;
	assert(degree == circuitV[0].size());
	//vector<mie::Vuint> Zx;
	vector<mie::Vuint> Hx;

	if (argv[1][0] == '0') {
		cout << "Gates     = " << NUM_GATE << endl;
		cout << "Wires     = " << NUM_WIRE << endl;
		cout << "Variables = " << NUM_WITNESS << endl;
		cout << "M         = " << m << endl;
		cout << "d         = " << circuitV[0].size() << endl;
		//verifyCircuitChecker();
		ofstream myfile;
		myfile.open("circuit.txt", ofstream::out | ofstream::app);
		for (int i = 0; i < proverMatrix.size(); i++) {
			myfile << "[ ";
			for (int j = 0; j < proverMatrix[i].size(); j++) {
				myfile << proverMatrix[i][j] << " ";
			}
			myfile << "]" << endl;
		}
		myfile.flush();
		myfile.close();
		return 0;
	} else {
		ifstream myfile;
		myfile.open("poly.txt");
		myfile >> degree;
		vector<mie::Vuint> polynom;
		for (int j = 0; j < degree; j++) {
			string coef;
			myfile >> coef;
			polynom.push_back(mie::Vuint(coef));
		}
		Vx.push_back(polynom);

		polynom.clear();
		for (int j = 0; j < degree; j++) {
			string coef;
			myfile >> coef;
			polynom.push_back(mie::Vuint(coef));
		}
		Wx.push_back(polynom);

		for (int i = 0; i < proverMatrix.size(); i++) {
			polynom.clear();
			for (int j = 0; j < degree; j++) {
				string coef;
				myfile >> coef;
				polynom.push_back(mie::Vuint(coef));
			}
			proverPoly.push_back(polynom);
		}

		for (int j = 0; j < degree - 1; j++) {
			string coef;
			myfile >> coef;
			Hx.push_back(mie::Vuint(coef));
		}
		myfile >> interpolateTime;
		myfile >> proverComputeTime;
		myfile >> CRSInterpolateTime;
		myfile.close();
	}

	gTable = new Fp**[degree>>2];
	gTable2 = new Fp2**[degree>>2];
	for (int i = 0; i < (degree>>2); i++) {
		gTable[i] = new Fp*[16];
		gTable2[i] = new Fp2*[16];
		for (int j = 0; j < 16; j++) {
			gTable[i][j] = new Fp[3];
			gTable2[i][j] = new Fp2[3];
		}
	}
	Vv = new mie::Vuint[degree];

	Fp (*sigmaToJEnc)[3] = new Fp[degree][3];
	Fp (*alphaSigmaToJEnc)[3] = new Fp[degree][3];;
	Fp2 (*sigmaToJEnc2)[3] = new Fp2[degree][3]; // second group

	Fp (*betaVViSigmaEnc)[3] = new Fp[m-1][3];
	Fp (*betaWWiSigmaEnc)[3] = new Fp[m-1][3];
	Fp betaVZSigmaEnc[3], betaWZSigmaEnc[3];
	Fp (*ViSigmaEnc)[3] = new Fp[m-1][3];
	Fp (*WiSigmaEnc)[3] = new Fp[m-1][3];

	mie::Vuint* vHat = new mie::Vuint[degree];
	mie::Vuint* wHat = new mie::Vuint[degree];
	// hHat should be degree - 1, but we pad it to degree for scalarmultimult
	mie::Vuint* hHat = new mie::Vuint[degree];

	int mult4 = 2 * (m - 1) + 2;
	while (mult4 % 4 != 0) mult4++;
	Fp (*phiYMult)[3] = new Fp[mult4][3];
	mie::Vuint* phiYExp = new mie::Vuint[mult4];

	clock_t crsComputeBegin = clock();
	//Zx.resize(2 * totalWidthG + totalWidthW, mie::Vuint(0));//= calcZx(2 * totalWidthG + totalWidthW);
	mie::Vuint alpha(rand()), sigma(rand()), betaV(rand()), betaW(rand()), gamma(rand());
	while(sigma <= degree) {
		// sigma has to be larger than degree, due to the fact that
		// in computing zSigma, we use division of (sigma - i)
		// where i is between 1 to degree
		sigma *= rand();
	}
	mie::Vuint v0Sigma(0), alphav0Sigma(0), w0Sigma(0), alphaw0Sigma(0), zSigma(0), alphazSigma(0);
	zSigma = calcZSigma(sigma, degree);
	vwSigma = calcVWSigma(sigma, zSigma, degree);
	circuitV.clear();
	circuitW.clear();
	v0Sigma = vwSigma[0];
	w0Sigma = vwSigma[m];
	alphazSigma = zSigma * alpha;
	alphazSigma %= Param::r;
	alphav0Sigma = v0Sigma * alpha;
	alphav0Sigma %= Param::r;
	alphaw0Sigma = w0Sigma * alpha;
	alphaw0Sigma %= Param::r;
	Fp V0[3], V0star[3], W0[3], W0star[3], Z[3], Zstar[3];
	Fp2 W02[3], Z2[3]; //second group
	ecop::ScalarMult(V0, g, v0Sigma);
	ecop::ScalarMult(V0star, g, alphav0Sigma);
	ecop::ScalarMult(W0, g, w0Sigma);
	ecop::ScalarMult(W0star, g, alphaw0Sigma);
	ecop::ScalarMult(Z, g, zSigma);
	ecop::ScalarMult(Zstar, g, alphazSigma);
	ecop::ScalarMult(W02, g2, w0Sigma);
	ecop::ScalarMult(Z2, g2, zSigma);

	mie::Vuint sigmaToJ(1);
	for (int i = 0; i < degree; i++) {
		ecop::ScalarMult(sigmaToJEnc[i], g, sigmaToJ);
		ecop::ScalarMult(sigmaToJEnc2[i], g2, sigmaToJ);
		mie::Vuint alphaSigmaToJ(alpha * sigmaToJ);
		alphaSigmaToJ %= Param::r;
		ecop::ScalarMult(alphaSigmaToJEnc[i], g, alphaSigmaToJ);
		sigmaToJ *= sigma;
		sigmaToJ %= Param::r;
		if (DEBUG_MODE) {
			Fp temp[3], temp2[3];
			ecop::ScalarMult(temp, sigmaToJEnc[i], alpha);
			ecop::NormalizeJac(temp, temp);
			ecop::NormalizeJac(temp2, alphaSigmaToJEnc[i]);
			for (int j = 0; j < 3; j++) {
				assert(temp[j] == temp2[j]);
			}
		}
	}

	for (int i = 1; i < m; i++) {
		mie::Vuint betaVViSigma(betaV * vwSigma[i]);
		betaVViSigma %= Param::r;
		mie::Vuint betaWWiSigma(betaW * vwSigma[m + i]);
		betaWWiSigma %= Param::r;
		ecop::ScalarMult(betaVViSigmaEnc[i-1], g, betaVViSigma);
		ecop::ScalarMult(betaWWiSigmaEnc[i-1], g, betaWWiSigma);
		ecop::ScalarMult(ViSigmaEnc[i-1], g, vwSigma[i]);
		ecop::ScalarMult(WiSigmaEnc[i-1], g, vwSigma[m + i]);
	}
	mie::Vuint betaVZSigma(betaV * zSigma);
	betaVZSigma %= Param::r;
	ecop::ScalarMult(betaVZSigmaEnc, g, betaVZSigma);
	mie::Vuint betaWZSigma(betaW * zSigma);
	betaWZSigma %= Param::r;
	ecop::ScalarMult(betaWZSigmaEnc, g, betaWZSigma);

	// for verifier's CRS (second group)
	Fp2 oneEnc[3], alphaEnc[3], gammaEnc[3], betaVGammaEnc[3], betaWGammaEnc[3];
	ecop::ScalarMult(oneEnc, g2, mie::Vuint(1));
	ecop::ScalarMult(alphaEnc, g2, alpha);
	ecop::ScalarMult(gammaEnc, g2, gamma);
	ecop::ScalarMult(betaVGammaEnc, g2, betaV * gamma);
	ecop::ScalarMult(betaWGammaEnc, g2, betaW * gamma);
	clock_t crsComputeEnd = clock();
	printf("CRS interpolation time: %.3f sec(s)\n", CRSInterpolateTime);
	printf("CRS computation time: %.3f sec(s)\n", double(crsComputeEnd - crsComputeBegin) / CLOCKS_PER_SEC);
	cout << "end of CRS computation" << endl;
	// Begin Prover's Computation
	printf("Prover's interpolation time: %.3f sec(s)\n", (double)(end - begin) / CLOCKS_PER_SEC + interpolateTime);
	fflush(stdout);
	begin = clock();

	assert(proverPoly.size() == 2);
	//vector<mie::Vuint> VDagger = addPoly(Vx[0], proverPoly[0]);
	//vector<mie::Vuint> WDagger = addPoly(Wx[0], proverPoly[1]);
	//Hx.resize(degree - 1, mie::Vuint(0));//= calcHx(VDagger, WDagger, Zx);

	Fp proverV[3], proverVStar[3], proverW[3], proverWStar[3], proverH[3], proverHStar[3];
	Fp2 proverW2[3]; // second group
	assert(proverPoly[0].size() == degree);
	assert(proverPoly[1].size() == degree);
	assert(Hx.size() == degree - 1);

	for (int i = 0; i < degree; i++) {
		vHat[i] = proverPoly[0][degree - i - 1];
		wHat[i] = proverPoly[1][degree - i - 1];
	}
	for (int i = 0; i < degree - 1; i++) {
		hHat[i] = Hx[degree - i - 2];
	}
	hHat[degree - 1] = mie::Vuint(0); // padding

	clock_t multiExpBegin = clock();
	ScalarMultiMult(proverV, sigmaToJEnc, vHat, degree);
	clock_t multiExpEnd = clock();
	printf("Single multi-exp time: %.3f sec(s)\n", (double)(multiExpEnd - multiExpBegin) / CLOCKS_PER_SEC);
	fflush(stdout);
	ScalarMultiMult(proverVStar, alphaSigmaToJEnc, vHat, degree);
	ScalarMultiMult(proverW, sigmaToJEnc, wHat, degree);
	ScalarMultiMult(proverWStar, alphaSigmaToJEnc, wHat, degree);
	ScalarMultiMult(proverW2, sigmaToJEnc2, wHat, degree);
	ScalarMultiMult(proverH, sigmaToJEnc, hHat, degree);
	ScalarMultiMult(proverHStar, alphaSigmaToJEnc, hHat, degree);

	mie::Vuint rV(rand()), rW(rand());
	Fp phiV[3], phiVStar[3], phiW[3], phiWStar[3], phiH[3], phiHStar[3];
	Fp temp[3];
	Fp2 temp2[3];
	Fp2 phiW2[3]; // second group
	clock_t expBegin = clock();
	ecop::ScalarMult(temp, Z, rV); // temp = Z^rv
	clock_t expEnd = clock();
	printf("Single exp time: %.3f sec(s)\n", (double)(expEnd - expBegin) / CLOCKS_PER_SEC);
	fflush(stdout);
	ecop::ECAdd(phiV, proverV, temp);
	ecop::ScalarMult(temp, Zstar, rV); // temp = Zstar^rv
	ecop::ECAdd(phiVStar, proverVStar, temp);
	ecop::ScalarMult(temp, Z, rW); // temp = Z^rw
	ecop::ECAdd(phiW, proverW, temp);
	ecop::ScalarMult(temp, Zstar, rW); // temp = Zstar^rw
	ecop::ECAdd(phiWStar, proverWStar, temp);
	Fp (*phiHMult)[3] = new Fp[4][3];
	Fp (*phiHStarMult)[3] = new Fp[4][3];
	mie::Vuint* phiHExp = new mie::Vuint[4];
	phiHExp[0] = mie::Vuint(1);
	phiHExp[1] = mie::Vuint(rV);
	phiHExp[2] = mie::Vuint(rW);
	phiHExp[3] = mie::Vuint(rV * rW);
	phiHExp[3] %= Param::r;
	ecop::copy(phiHMult[0], proverH);
	ecop::ECAdd(phiHMult[1], W0, proverW);
	ecop::ECAdd(phiHMult[2], V0, proverV);
	ecop::copy(phiHMult[3], Z);
	ScalarMultiMult(phiH, phiHMult, phiHExp, 4);
	ecop::copy(phiHStarMult[0], proverHStar);
	ecop::ECAdd(phiHStarMult[1], W0star, proverWStar);
	ecop::ECAdd(phiHStarMult[2], V0star, proverVStar);
	ecop::copy(phiHStarMult[3], Zstar);
	ScalarMultiMult(phiHStar, phiHStarMult, phiHExp, 4);
	ecop::ScalarMult(temp2, Z2, rW); // temp2 = Z2^rw
	ecop::ECAdd(phiW2, proverW2, temp2);

	Fp phiY[3];
	// make the witness map first
	map<int, int> witnessMapA;
	map<int, int> witnessMapB;
	for (int i = 0; i < NUM_WITNESS; i++) {
		witnessMapA[witnessA[i][0]] = witnessA[i][1];
		witnessMapB[witnessB[i][0]] = witnessB[i][1];
	}
	assert(NUM_WITNESS * 2 == m - 1); // positive + negative literal
	assert(VLabelling.size() == m);
	for (int i = 0; i < m - 1; i++) {
		ecop::copy(phiYMult[2 * i], betaVViSigmaEnc[i]);
		ecop::copy(phiYMult[2 * i + 1], betaWWiSigmaEnc[i]);
		if (witnessMapA.find(VLabelling[i+1]) != witnessMapA.end()) {
			mie::Vsint val(witnessMapA[VLabelling[i+1]]);
			if (val.isNegative()) {
				val += Param::r;
			}
			assert(val == val.abs());
			phiYExp[2 * i] = val.abs();
		} else {
			phiYExp[2 * i] = mie::Vuint(0);
		}
		if (witnessMapB.find(WLabelling[i+1]) != witnessMapB.end()) {
			mie::Vsint val(witnessMapB[WLabelling[i+1]]);
			if (val.isNegative()) {
				val += Param::r;
			}
			assert(val == val.abs());
			phiYExp[2 * i + 1] = val.abs();
		} else {
			phiYExp[2 * i + 1] = mie::Vuint(0);
		}
	}
	ecop::copy(phiYMult[2 * (m - 1)], betaVZSigmaEnc);
	ecop::copy(phiYMult[2 * (m - 1) + 1], betaWZSigmaEnc);
	phiYExp[2 * (m - 1)] = mie::Vuint(rV);
	phiYExp[2 * (m - 1) + 1] = mie::Vuint(rW);
	for (int i = 2 * (m - 1) + 2; i < mult4; i++) {
		phiYExp[i] = mie::Vuint(0);
	}
	ScalarMultiMult(phiY, phiYMult, phiYExp, mult4);
	end = clock();
	printf("Prover's exponentiation time: %.3f sec(s)\n", ((double)(end - begin) / CLOCKS_PER_SEC) + proverComputeTime);
	cout << "end of Prover's computation" << endl;

	// Begin Verifier's Computation
	begin = clock();

	Fp12 e11, e21, e31, e41, e51, e12, e22, e32, e42, e521, e522;
	ecop::ECAdd(temp, V0, phiV);
	ecop::ECAdd(temp2, W02, phiW2);
	ecop::NormalizeJac(temp, temp);
	ecop::NormalizeJac(temp2, temp2);
	ecop::NormalizeJac(phiH, phiH);
	ecop::NormalizeJac(Z2, Z2);
	opt_atePairing<Fp>(e11, temp2, temp);
	opt_atePairing<Fp>(e12, Z2, phiH);
	if (e11 == e12) {
		cout << "First pairing verified!" << endl;
	} else {
		cout << "Fail on first pairing!" << endl;
	}

	ecop::NormalizeJac(alphaEnc, alphaEnc);
	ecop::NormalizeJac(phiVStar, phiVStar);
	ecop::NormalizeJac(oneEnc, oneEnc);
	ecop::NormalizeJac(phiV, phiV);
	opt_atePairing<Fp>(e21, oneEnc, phiVStar);
	opt_atePairing<Fp>(e22, alphaEnc, phiV);
	if (e21 == e22) {
		cout << "Second pairing verified!" << endl;
	} else {
		cout << "Fail on second pairing!" << endl;
	}

	ecop::NormalizeJac(phiWStar, phiWStar);
	ecop::NormalizeJac(phiW, phiW);
	opt_atePairing<Fp>(e31, oneEnc, phiWStar);
	opt_atePairing<Fp>(e32, alphaEnc, phiW);
	if (e31 == e32) {
		cout << "Third pairing verified!" << endl;
	} else {
		cout << "Fail on third pairing!" << endl;
	}

	ecop::NormalizeJac(phiHStar, phiHStar);
	opt_atePairing<Fp>(e41, oneEnc, phiHStar);
	opt_atePairing<Fp>(e42, alphaEnc, phiH);
	if (e41 == e42) {
		cout << "Fourth pairing verified!" << endl;
	} else {
		cout << "Fail on fourth pairing!" << endl;
	}

	ecop::NormalizeJac(gammaEnc, gammaEnc);
	ecop::NormalizeJac(phiY, phiY);
	ecop::NormalizeJac(betaVGammaEnc, betaVGammaEnc);
	ecop::NormalizeJac(betaWGammaEnc, betaWGammaEnc);
	opt_atePairing<Fp>(e51, gammaEnc, phiY);
	opt_atePairing<Fp>(e521, betaVGammaEnc, phiV);
	opt_atePairing<Fp>(e522, betaWGammaEnc, phiW);
	if (e51 == (e521 * e522)) {
		cout << "Fifth pairing verified!" << endl;
	} else {
		cout << "Fail on fifth pairing!" << endl;
	}
	end = clock();
	printf("Verifier's computation time: %.3f sec(s)\n", (double)(end - begin) / CLOCKS_PER_SEC);
	cout << "end of Verifier's computation" << endl;
	fflush(stdout);
	return 0;
}