#include <cstdlib>
# include<iostream>
# include<string>
# include<random>
# include<cmath>
# include<array>
#include <fstream>
using namespace std;

typedef long long int Lint; // 64 bits 
typedef double Ldouble;
struct security {
  int num_shares;
  int num_required;
};

struct shareStruct {
  int x;
  Lint y;
};

bool isPrime(Lint n) {
  int flag = 0;
  for (int i = 2; i <= n / i; ++i) {
    if (n % i == 0) {
      flag = 1;
      break;
    }
  }
  if (flag == 0) return true;
  else return false;
}

Lint genRandInt(int n) {
    // Returns a random number 
    // between 2**(n-1)+1 and 2**n-1
    //long max = (long)powl(2, n) - 1;
    //long min = (long)powl(2, n - 1) + 1;
    long max = (long)pow(2, n) - 1;
    long min = (long)pow(2, n - 1) + 1;
    Lint result = min + (rand() % ( max - min + 1 ) );
    return result;
}

Lint genPrime() {
  Lint prime = 10;

  while (isPrime(prime) == false) {
    int complexity = 50;
    prime = genRandInt(complexity);
  }
  return prime;
}

int* encodeSecret(int* poly, const int secret, const int num_required) {
  poly[num_required-1] = secret;
  return poly;
}

Lint getPolyY(const int* poly, int poly_len, int poly_x, const Lint prime) {
  Lint total = 0;
  Lint poly_y = 0;

  for (int i=0; i<poly_len+1; i++) {
    int power = poly_len - i;
    int coefficient = poly[i];
    poly_y = coefficient * pow(poly_x, power);
    total = total + poly_y;
  }

  return total;
}

shareStruct* genShares(int num_shares, int num_required, const int* poly, const Lint prime){
  shareStruct* shares = new shareStruct[num_shares];
  for (int i=1; i<=num_shares; i++) {
    shareStruct share;
    share.x = i;
    share.y = getPolyY(poly, num_required-1, share.x, prime);
    shares[i-1] = share;
  }
  return shares;
}

int* genPoly(int degree, const Lint prime, const Lint secret) {
  int* poly = new int[degree];

  for (int i = 0; i < degree; i++) {
    int random_num = genRandInt(10);
    poly[i] = prime % random_num;
  }
  return poly;
}

// solving polynomials
struct inputStruct {
  int required;
  shareStruct* shares;
};

struct polyTerm {
  Lint coefficient;
  int power;
};

struct linearEquation {
  shareStruct point;
  polyTerm* terms;
};

linearEquation* constructEquations(const int required, shareStruct shares[]) {
  linearEquation* equations = new linearEquation[required];
  shareStruct share;
  polyTerm term;

  for (int i = 0; i < required; i++) {
    share = shares[i];
    linearEquation equation;
    polyTerm* terms = new polyTerm[required];
    
    for (int j = 0; j < required; j++) {
      term.power = required - 1 - j;
      terms[j] = term;
    }

    equation.terms = terms;
    equation.point.x = share.x;
    equation.point.y = share.y;
    
    equations[i] = equation;
    // dont delete terms from memory as its referanced in equations
  }
  return equations;
}

struct matrix{
  Lint** matrix;
  int dimension_x;
  int dimension_y;
};
struct matrix_system {
  matrix A;
  matrix B;
  matrix X;
};

matrix_system formMatrix(const linearEquation* equations, int required) {
  Lint** matrixA = new Lint*[required];
  Lint** matrixB = new Lint*[required];

  for (int i=0; i < required; i++) {
    linearEquation equation = equations[i];
    Lint* lineA = new Lint[required];
    for (int j=0; j < required; j++) {
      lineA[j] = pow(equation.point.x, equation.terms[j].power);
    }
    matrixA[i] = lineA;

    Lint* lineB = new Lint[1];
    lineB[0] = equation.point.y;
    matrixB[i] = lineB;
  }

  matrix matrixA_data; matrix matrixB_data;
  matrixA_data.matrix = matrixA; matrixB_data.matrix = matrixB;

  matrixA_data.dimension_x = required; matrixB_data.dimension_x = 1;
  matrixA_data.dimension_y = required; matrixB_data.dimension_y = required;

  matrix_system matricies;
  matricies.A = matrixA_data; matricies.B = matrixB_data;

  return matricies;
}

Lint** findMinor(Lint** matrixA, const int dimension, const int pos_x, const int pos_y) {
  Lint** matrixB = new Lint*[dimension-1];
  int matrixB_pos_x = 0; int matrixB_pos_y = 0;

  for (int i=0; i<dimension; i++) {
    Lint* line = new Lint[dimension-1];
    for (int j=0; j<dimension; j++) {
      if (i != pos_y and j != pos_x) {
        line[matrixB_pos_x] = matrixA[i][j];
        matrixB_pos_x++;
      }
    }
    if (matrixB_pos_x != 0) {
      matrixB[matrixB_pos_y] = line;
      matrixB_pos_y++;
    }
    else {
      delete[] line;
    }
    matrixB_pos_x = 0;
  }

  return matrixB;
}

Lint findDet(Lint** matrixA, const int dimension) {
  Lint det = 0;
  if (dimension == 0) {
    det = 1;
  }
  else if (dimension == 1) {
    det = matrixA[0][0];
  }
  else if (dimension == 2) {
    det = matrixA[0][0] * matrixA[1][1] - matrixA[0][1] * matrixA[1][0];
  }
  else {
    for (int i=0; i<dimension; i++) {
      // reuse form matrix? pottentially split it up into formMatrixA and formMatrixB?
      Lint** matrixB = findMinor(matrixA, dimension, i, 0);
      Lint matrixB_det = findDet(matrixB, dimension-1);
      Lint term = matrixA[0][i] * matrixB_det;

      if ((i+1)%2 == 0) {
        term = 0-term;
      }
      det = det + term;
    }
  }

  return det;
}

matrix formMatrixCofactors(Lint** matrixA, const int dimension) {
  Lint** matrixB = new Lint*[dimension];

  for (int i=0; i<dimension; i++) {
    Lint* line = new Lint[dimension];

    int sign = 1;
    if ((i+1)%2 == 0) {
      sign = -1;
    }
    for (int j=0; j<dimension; j++) {
      Lint** minor = findMinor(matrixA, dimension, j, i);
      Lint cofactor = findDet(minor, dimension-1) * sign;
      sign = -sign;
      line[j] = cofactor;
    }
    matrixB[i] = line;
  }

  matrix matrix_data; matrix_data.matrix = matrixB;
  matrix_data.dimension_x = dimension; matrix_data.dimension_y = dimension;
  return matrix_data;
}

matrix transposeMatrix(Lint** cofactors, const int dimension) {
  Lint** matrixB = new Lint*[dimension];

  for (int i=0; i<dimension; i++) {
    Lint* line = new Lint[dimension];
    for (int j=0; j<dimension; j++) {
      line[j] = cofactors[j][i];
    }
    matrixB[i] = line;
  }

  matrix matrixB_data; matrixB_data.matrix = matrixB;
  matrixB_data.dimension_x = dimension; matrixB_data.dimension_y = dimension;
  return matrixB_data;
}

struct float_matrix{
  Ldouble** matrix;
  int dimension_x;
  int dimension_y;
};
struct float_matrix_system {
  matrix A;
  matrix B;
  matrix X;
};

float_matrix multiplyConstant(matrix matrixA_data, const int dimension, const Lint det) {
  Ldouble** matrixB = new Ldouble*[dimension];
  Lint** matrixA = matrixA_data.matrix;

  for (int i=0; i<dimension; i++) {
    Ldouble* line = new Ldouble[dimension];
    for (int j=0; j<dimension; j++) {
      line[j] = (1.0/det) * matrixA[i][j];
    }
    matrixB[i] = line;
  }
  float_matrix matrixB_data; matrixB_data.matrix = matrixB;
  matrixB_data.dimension_x = matrixA_data.dimension_x; matrixB_data.dimension_y = matrixA_data.dimension_y;
  
  return matrixB_data;
}

float_matrix multiplyMatricies(float_matrix inverseA_data, matrix matrixB_data) {
  int dimension_x = inverseA_data.dimension_x;
  int dimension_y = inverseA_data.dimension_y;

  Ldouble** matrixC = new Ldouble*[matrixB_data.dimension_y];
  Ldouble** inverseA = inverseA_data.matrix;
  Lint** matrixB = matrixB_data.matrix;

  for (int i=0; i<dimension_y; i++) {
    Ldouble* line = new Ldouble[0];
    Ldouble result = 0;
    for (int j=0; j<dimension_x; j++) {
      result = result + inverseA[i][j] * matrixB[j][0];
    }
    line[0] = result;
    matrixC[i] = line;
  }
  float_matrix matrixC_data; matrixC_data.matrix = matrixC;
  matrixC_data.dimension_x = matrixB_data.dimension_x; matrixC_data.dimension_y = matrixB_data.dimension_y;
  
  return matrixC_data;
}

Lint** StructToArray(shareStruct* struct_array, int len_array) {
  Lint** array = new Lint*[len_array];
  for (int i=0; i<len_array; i++) {
    array[i] = new Lint[2];
    array[i][0] = struct_array[i].x;
    array[i][1] = struct_array[i].y;

  }
  return array;
}

shareStruct* ArrayToStruct(Lint** array, int len_array) {
  shareStruct* share_array = new shareStruct[len_array];
  for (int i=0; i<len_array; i++) {
    shareStruct share;
    share.x = array[i][0];
    share.y = array[i][1];
    share_array[i] = share;
  }
  return share_array;
}

void writeShares(shareStruct* shares, const int num_shares, const int num_required, string root_path) {
  cout << root_path << endl;
  for (int i=0; i<num_shares; i++) {
    shareStruct share = shares[i];
    string share_path = root_path + "share-" + to_string(share.x) + ".txt";
    ofstream share_file(share_path);
    share_file << "Share number: " << share.x << endl;
    share_file << "Share secret: " << share.y << endl;
    share_file << "Minimum share required: " << to_string(num_required) << endl << endl;
    share_file << "IMPORTANT: Please remind your admin that its there job to distribute and delete shares from the server";
  }
}

extern "C" Lint solveInternal(shareStruct* shares, int required) {
  inputStruct inputs;
  inputs.shares = shares;
  inputs.required = required;

  linearEquation* equations = new linearEquation[inputs.required];
  equations = constructEquations(inputs.required, inputs.shares);

  matrix_system matricies = formMatrix(equations, inputs.required);
  delete[] equations;
  Lint det = findDet(matricies.A.matrix, matricies.A.dimension_x);

  matrix cofactors = formMatrixCofactors(matricies.A.matrix, matricies.A.dimension_x);
  matrix transposition = transposeMatrix(cofactors.matrix, cofactors.dimension_x);

  float_matrix inverseA = multiplyConstant(transposition, transposition.dimension_x, det);
  float_matrix matrixC = multiplyMatricies(inverseA, matricies.B);

  Lint secret = matrixC.matrix[matrixC.dimension_y-1][0];
  return secret;
}

extern "C" void newSecretInternal(const Lint secret, const int num_shares, const int num_required, char* root_path) {
  string str(root_path);
  const Lint prime = genPrime();
  int* poly = genPoly(num_required-1, prime, secret);

  poly = encodeSecret(poly, secret, num_required);
  shareStruct* shares = genShares(num_shares, num_required, poly, prime);

  writeShares(shares, num_shares, num_required, root_path);
}


int main() {
}