   }
   std::cout << "Average difference between the two is " << sum/(nr_rows * nr_rows) << std::endl;
   return(sum/(nr_rows * nr_rows));
+}
+
+
+void estimateSize(const mpf_class *a, int &sizeOfArray, int &maxExp, const int nr_rows, const int nr_cols, const int ownLimbSize) {
+  assert(ownLimbSize <= INT64L);
+
+  int minSize = 0;
+  maxExp = a[0].get_mpf_t()->_mp_exp;
+  for(int i = 0; i < nr_rows; ++i)
+    for(int j = 0; j < nr_cols; ++j) {
+      int toCmp =  a[j * nr_rows + i].get_mpf_t()->_mp_exp;
+      if (toCmp >= maxExp) {
+        maxExp = toCmp;
+        minSize = a[j * nr_rows + i].get_mpf_t()->_mp_size;
+      }
+    }
+
+  maxExp *= INT64L;
+  sizeOfArray = (INT64L * minSize) / ownLimbSize;
+  if ((INT64L * minSize) % ownLimbSize != 0) sizeOfArray += 1;
+}
+
+
+// All arrays need to already have allocated memory: 
+// *a : array of mpf_class with size nr_rows * nr_cols
+// *d_aS: array allocated in GPU with size sizeOfArray * nr_rows * nr_cols. Note that we flatten this 
+//        array of matrices in order to speed up the access in the GPU.
+//        Thus, to access the k-th limb from the (i, j) matrix entry one calls d_aS[k * nr_rows * nr_cols + j * nr_rows + i] 
+// *tmpA: temporary array for doubles used for transfers that needs to be allocated nr_rows * nr_cols entries
+// sizeOfArray: store number of own limbs that are saved for each entry
+// nr_rows: number of matrix rows
+// nr_cols: number of matrix columns 
+// maxExp:  maximum power of 2 for the leading most limb among all the entries of the array *a 
+void generateLongMatrixFromGMPMatrix_GPU(const mpf_class *a, double *d_aS, double *tmpA, int &sizeOfArray, 
+					    const int nr_rows, const int nr_cols, int &maxExp, 
+					    const int ownLimbSize) {
+  assert(ownLimbSize <= INT64L);
+  timeval t1, t2;
+  double etAlloc = 0;
+  for (int k = 0; k < sizeOfArray - 1; k++) {
+    for(int i = 0; i < nr_rows; ++i) {
+      for(int j = 0; j < nr_cols; ++j) {
+	int size = a[j * nr_rows + i].get_mpf_t()->_mp_size;
+	int realExp = a[j * nr_rows + i].get_mpf_t()->_mp_exp * INT64L;
+	int padding = (maxExp - realExp) / ownLimbSize;
+	int padBitOffset = (maxExp - realExp) % ownLimbSize;
+	int sign  = getSign(a[j * nr_rows + i]);
+	if(k > padding) {
+	  int leftGmpLimb  = size - 1 - (((k - padding) * ownLimbSize - padBitOffset) / INT64L);
+	  int leftBit      = ((k - padding) * ownLimbSize - padBitOffset) % INT64L;
+	  int rightGmpLimb = size - 1 - (((k - padding + 1) * ownLimbSize - padBitOffset) / INT64L);
+	  int rightBit     = ((k - padding + 1) * ownLimbSize - padBitOffset) % INT64L;
+	  if (leftGmpLimb == rightGmpLimb) {
+	    long long tmp  = a[j * nr_rows + i].get_mpf_t()->_mp_d[leftGmpLimb];
+	    tmpA[j * nr_rows + i] = sign * getBitsFromOneLong(tmp, leftBit, rightBit);
+	  } else {
+	    long long tmp1 = a[j * nr_rows + i].get_mpf_t()->_mp_d[leftGmpLimb];
+	    long long tmp2 = a[j * nr_rows + i].get_mpf_t()->_mp_d[rightGmpLimb];
+	    tmpA[j * nr_rows + i] = sign * getBitsFromTwoLong(tmp1, tmp2, leftBit, rightBit);
+	  }
+	} else if (k < padding){
+	  tmpA[j * nr_rows + i] = 0;
+	} else {
+	  long long tmp  = a[j * nr_rows + i].get_mpf_t()->_mp_d[size - 1];
+	  tmpA[j * nr_rows + i] = sign * getBitsFromOneLong(tmp, 0, ownLimbSize - padBitOffset);
+	}
+      } 
+    }
+    // Transfer matrix to memeory
+    gettimeofday(&t1, NULL);
+    cudaMemcpy(&d_aS[k * nr_rows * nr_cols], tmpA, nr_rows * nr_cols * sizeof(double), cudaMemcpyHostToDevice);
+    cudaThreadSynchronize();
+    gettimeofday(&t2, NULL);
+    etAlloc +=  (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
+  }
+
+  for(int i = 0; i < nr_rows; ++i) {
+    for(int j = 0; j < nr_cols; ++j) {
+      int realExp = a[j * nr_rows + i].get_mpf_t()->_mp_exp * INT64L;
+      int padding = (maxExp - realExp) / ownLimbSize;
+      int padBitOffset = (maxExp - realExp) % ownLimbSize;
+      int sign  = getSign(a[j * nr_rows + i]);
+      int leftBit = ((sizeOfArray - padding - 1) * ownLimbSize + padBitOffset) % INT64L;
+      int tmp = a[j * nr_rows + i].get_mpf_t()->_mp_d[0];
+      tmpA[j * nr_rows + i] = sign * getBitsFromOneLong(tmp, leftBit, INT64L);
+    }
+  }
+  cudaMemcpy(&d_aS[(sizeOfArray - 1) * nr_cols * nr_rows], tmpA, nr_rows * nr_cols * sizeof(double),cudaMemcpyHostToDevice);
+  printf("Actual allocation time to GPU = %fms\n", etAlloc);
+}
+
+
+
+void checkFit() {
+}
+
+// Returns matrix product c = a.b where each entry is of the type mpf_class
+//
+// All arrays need to already have allocated memory:                                                                  
+// *c : array of mpf_class with size nr_rowsA * nr_colsB                   
+// *a : array of mpf_class with size nr_rowsA * nr_colsA
+// *b : array of mpf_class with size nr_colsA * nr_colsB
+// *d_aS: array allocated in GPU with size sizeOfArray * nr_rowsA * nr_colsA. Note that we flatten this               
+//        array of matrices in order to speed up the access in the GPU.                                              
+//        Thus, to access the k-th limb from the (i, j) matrix entry one calls d_aS[k * nr_rowsA * nr_colsA + j * nr_rowsA + i]
+// *d_bS: array allocated in GPU with size sizeOfArray * nr_colsA * nr_colsB. Note that we flatten this                                                              
+//        array of matrices in order to speed up the access in the GPU.                                                                                              
+//        Thus, to access the k-th limb from the (i, j) matrix entry one calls d_bS[k * nr_colsA * nr_colsB + j * nr_colsA + i] 
+// *tmpTransferLongToGMP : temporary array of doubles used for transfers that needs max(nr_rowsA * nr_colsA, nr_colsA * nr_colsB) entries                            
+// *tmp: temporary array of 64-bit vars used for transfers. Needs to be allocated nr_rowsA * nr_colsB entries
+// *d_prodRes: temporary matrix allocated in GPU memory which handles individual multiplications of matrices
+// *d_res : temporary matrix allocated in GPU memory in which we sum up all individual multiplications of matrices
+// *d_rem : in case we encounter overflow in the matrix d_res we use the entries in the matrix d_rem to store the remainders. 
+//          Note that this matrix is only used when using the function vecAdd__wRem. When using vecAdd__wSign we assume that 
+//          overflow never occurs. This is the case if the number of own generated limbs is greater than 1024. This correpons 
+//          to a GMP precision ~20,000. For current bootstrap applications such precision is not needed.
+// nr_rows_A: number of rows for matrix A as well as for matrix C                                                                                
+// nr_cols_A: number of columns for matrix A as well as number of rows for matrix B
+// nr_cols_B: number of columns for matrix B as well as for matrix C                                                                               
+// maxExp:  maximum power of 2 for the leading most limb among all the entries of the array *a 
+void matrixMultiplicationBasecase_cuBlas(const cublasHandle_t handle, mpf_class *c, mpf_class *a, mpf_class *b, 
+					 double *d_aS, double *d_bS, double *tmpTransferLongToGMP, long long *tmp,   
+					 double *d_prodRes, long long *d_res,  //long long *d_rem, 
+					 const int nr_rowsA, const int nr_colsA, const int nr_colsB) {
+  
+  int size_aS = 0; 
+  int size_bS = 0;
+  int expA = 0;
+  int expB = 0;
+  timeval t1, t2;
+  
+  int ownLimbSize = DOUBLE_MANT/2 - ceil(log2((double) nr_colsA));
+  gettimeofday(&t1, NULL);
+  estimateSize(a, size_aS, expA, nr_rowsA, nr_colsA, ownLimbSize);
+  estimateSize(b, size_bS, expB, nr_colsA, nr_colsB, ownLimbSize);
+  
+  generateLongMatrixFromGMPMatrix_GPU(a, d_aS, tmpTransferLongToGMP, size_aS, nr_rowsA, nr_colsA, expA, ownLimbSize);
+  generateLongMatrixFromGMPMatrix_GPU(b, d_bS, tmpTransferLongToGMP, size_bS, nr_colsA, nr_colsB, expB, ownLimbSize);
+  gettimeofday(&t2, NULL);
+  double etTransfer = (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
+  
+  int size = min(size_aS, size_bS);
+  std::cout << size_aS << " " << size_bS << " " << expA << " " << expB << std::endl;
+
+  for(int i = 0; i < nr_rowsA; ++i){
+    for(int j = 0; j < nr_colsB; ++j){
+      c[j * nr_rowsA +  i] = mpf_class("0.0");
+      tmp[j * nr_rowsA +  i] = 0;
+    }
+  }
+
+  // Number of threads in each thread block                                                                          
+  int blockSize = 256; // Make sure this number is good? I don't know how
+  // Number of thread blocks in grid                                                                                 
+  int gridSize = (int)ceil((float)(nr_rowsA * nr_colsB)/blockSize);
+  
+  double etMult = 0;
+  double etAdd = 0;
+  double etAlloc = 0;
+  double etBackAlloc = 0;
+  double etaddBack = 0;
+  for (int i = 0; i < size; i++) {
+    gettimeofday(&t1, NULL);
+    cudaMemcpy(d_res, tmp, nr_rowsA * nr_colsB * sizeof(long long), cudaMemcpyHostToDevice);
+    gettimeofday(&t2, NULL);
+    etAlloc += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
+    //cudaMemcpy(d_rem, tmp, nr_rowsA * nr_colsB * sizeof(long long), cudaMemcpyHostToDevice);
+    for (int j = 0; j < i + 1; j++) {
+      gettimeofday(&t1, NULL);
+      gpu_blas_mmul(handle, &d_aS[j * nr_rowsA * nr_colsA], &d_bS[(i - j) * nr_colsA * nr_colsB], d_prodRes, nr_rowsA, nr_colsA, nr_colsB);
+      cudaThreadSynchronize();
+      gettimeofday(&t2, NULL);
+      etMult += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
+      
+      gettimeofday(&t1, NULL);
+      vecAdd__wSign<<<gridSize, blockSize>>>(d_prodRes, d_res,  nr_rowsA * nr_colsB);
+      cudaThreadSynchronize();
+      gettimeofday(&t2, NULL);
+      etAdd += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
+      // This is safe with overflow until there are 1024 of our own limbs that need to be summed up
+      // This case corresponds to a precision of ~20000 bits in GMP
+    }
+    gettimeofday(&t1, NULL);
+    cudaMemcpy(tmp, d_res, nr_rowsA * nr_colsB * sizeof(long long), cudaMemcpyDeviceToHost);
+    gettimeofday(&t2, NULL);
+    etBackAlloc += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
+    
+    gettimeofday(&t1, NULL);
+    addToGMPMatrix(c, tmp, nr_rowsA, nr_colsB, expA + expB - (i + 2) * ownLimbSize);
+    gettimeofday(&t2, NULL);
+    etaddBack += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
+    //cudaMemcpy(tmp, d_rem, nr_rowsA * nr_colsB * sizeof(long long),cudaMemcpyDeviceToHost);
+    //addToGMPMatrix(c, tmp, nr_rowsA, nr_colsB, 2 * exp - (i + 2) * ownLimbSize + (INT64L - 1));
+    for(int k = 0; k < nr_rowsA; ++k){
+      for(int l = 0; l < nr_colsB; ++l){
+	tmp[l *  nr_rowsA + k] = 0;
+      }
+    }
+  }
+  printf("Transfer GPU = %fms\n", etTransfer);
+  printf("Multiplication GPU = %fms\n", etMult);
+  printf("Addition GPU = %fms\n", etAdd);
+  printf("Alloc of zero to GPU = %fms\n", etAlloc);
+  printf("Transfer from GPU to host = %fms\n", etBackAlloc);
+  printf("Addition on CPU = %fms\n", etaddBack);
+}
+
+
+int estimateMaxGPUAllocation(double maxFrac, int nr_rows, int nr_cols) {
+  size_t free_byte ;
+  size_t total_byte ;
+  cudaError_t cuda_status = cudaMemGetInfo( &free_byte, &total_byte ) ;
+  if ( cudaSuccess != cuda_status ){
+    printf("Error: cudaMemGetInfo fails, %s \n", cudaGetErrorString(cuda_status) );
+    exit(1);
+  }
+  double free_db = ((double)free_byte) * maxFrac;
+  return (int) (8 * free_db/(nr_rows * nr_cols * INT64L));
 }
 
  

  rr.seed(time(NULL));
  for(int i = 0; i < nr_rows; i++){
    for(int j = 0; j < nr_cols; j++){
      std::cout << j << " " << i << std::endl;
      a[j * nr_rows +  i] = rr.get_f(500);
    }
  }

  }
} 

// TESTED                                                                                                                                                                                
void longToGMPMatrixDouble(mpf_class *a, double *toAdd, const int size_toAdd, const int nr_rows, const int nr_cols, const int ownLimbSize, const int whereWeStart) {
  for(int i = 0; i < nr_rows; ++i){
    for(int j = 0; j < nr_cols; ++j){
      a[j * nr_rows +  i] = mpf_class("0.0");
      for (int k = 0; k < size_toAdd; ++k) {
        a[j * nr_rows +  i]  = addToGMP(a[j * nr_rows +  i], (long long)toAdd[k * nr_rows * nr_cols + j * nr_rows +  i], whereWeStart - (k + 1) * ownLimbSize);
      }
    }
  }
}


// Now tested
void numberMultiplicationBasecase(mpf_class &c, const mpf_class a, const mpf_class b) {
  if(a.get_prec() != b.get_prec()) {
