   
   estimateSize(a, size, exp, nr_rowsA, nr_colsA, ownLimbSize);
 
-  if (size + 1 < maxNoMatrices) {
-    matrixMultSymmBasecase_cuBlas(handle, c, a, d_aS, tmpTransferLongToGMP, 
-				  tmp, sizeMatrix, realExpMatrix, signMatrix,  
-				  d_prodRes, d_res, nr_rowsA, nr_colsA, prec); 
+  if (size + 1 < maxNoMatrices) {    
+    // We need to allocate the matrix in memory here. This should be made more efficiently
+    generateLongMatrixFromGMPMatrix_GPU(a, d_aS, tmpTransferLongToGMP, sizeMatrix, realExpMatrix, signMatrix, 
+					size, nr_rowsA, nr_colsA, exp, ownLimbSize);
+    matrixMultSymmBasecaseArbPrec_cuBlas(handle, c, a, d_aS, tmpTransferLongToGMP, 
+					 tmp, sizeMatrix, realExpMatrix, signMatrix,  
+					 d_prodRes, d_res, nr_rowsA, nr_colsA, size, exp, ownLimbSize, prec); 
   } else {
-    toom2(handle, c, a, b, d_aS, d_bS, tmpTransferLongToGMP, 
-	  tmp, sizeMatrix, realExpMatrix, signMatrix,  
-	  d_prodRes, d_res, nr_rowsA, nr_colsA, nr_colsB);
+    toom2Symm(handle, c, a, d_aS, tmpTransferLongToGMP, 
+	      tmp, sizeMatrix, realExpMatrix, signMatrix,  
+	      d_prodRes, d_res, nr_rowsA, nr_colsA, nr_colsB, prec);
   }
 }
 
 // Implements Karatsuba algorithm for matrix multiplication c = a.b
 // We split each matrix element in two parts which are equal in length up to the size of one GMP limb. 
 void toom2(const cublasHandle_t handle, mpf_class *c, mpf_class *a, mpf_class *b, 
-	   double *d_aS, double *d_bS, double *tmpTransferLongToGMP, 
-	   long long *tmp, int *sizeMatrix, int *realExpMatrix, int *signMatrix,  
-	   double *d_prodRes, long long *d_res, const int nr_rowsA, const int nr_colsA, const int nr_colsB) {
+	   mpf_class *tmpC1, mpf_class *tmpC2, double *d_aS, double *d_bS, double *tmpTransferLongToGMP, 
+	   long long *tmp, int *sizeMatrix, int *realExpMatrix, int *signMatrix, double *d_prodRes, 
+	   long long *d_res, const int nr_rowsA, const int nr_colsA, const int nr_colsB, 
+	   const int size_aS, const int size_bS, const int leadingExp, const int whatOrder) {
+  mpf_class *a1, *a2, *b1, *b2; 
+  
+  int sizeA1 = (size_aS / 2 + size_aS % 2), sizeA2 = size_aS / 2, sizeB1 = (size_bS / 2 + size_bS % 2), sizeB2 = size_bS / 2; 
+ 
+  for(int i = 0; i < nr_colsA; ++i){
+    for(int k = 0; k < nr_rowsA; ++k){
+      if (sizeA1 + sizeA2 != abs(a[i *  nr_rowsA + k].get_mpf_t()->_mp_size))
+	std::cout << "Error?" << std::endl;
+      int sign =  a[i *  nr_rowsA + k].get_mpf_t()->_mp_size /
+	abs(a[i *  nr_rowsA + k].get_mpf_t()->_mp_size);
+      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_size = sign * sizeA1;
+      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_size = sign * sizeA2;
+      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_exp = a[i *  nr_rowsA + k].get_mpf_t()->_mp_exp; 
+      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_exp = (a[i *  nr_rowsA + k].get_mpf_t()->_mp_exp) - sizeA1;
+      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_d = &(a[i *  nr_rowsA + k].get_mpf_t()->_mp_d[sizeA2]); 
+      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_d = &(a[i *  nr_rowsA + k].get_mpf_t()->_mp_d[0]);
+    }
+  }
+  
+  for(int i = 0; i < nr_colsB; ++i){
+    for(int k = 0; k < nr_colsA; ++k){
+      if (sizeB1 + sizeB2 != abs(b[i *  nr_rowsA + k].get_mpf_t()->_mp_size))
+	std::cout << "Error?" << std::endl;
+      int sign =  b[i *  nr_rowsA + k].get_mpf_t()->_mp_size /
+	abs(b[i *  nr_rowsA + k].get_mpf_t()->_mp_size);
+      b1[i *  nr_colsA + k].get_mpf_t()->_mp_size = sign * sizeB1; 
+      b2[i *  nr_colsA + k].get_mpf_t()->_mp_size = sign * sizeB2;
+      b1[i *  nr_colsA + k].get_mpf_t()->_mp_exo = b[i *  nr_rowsA + k].get_mpf_t()->_mp_exp; 
+      b2[i *  nr_colsA + k].get_mpf_t()->_mp_exp = b[i *  nr_rowsA + k].get_mpf_t()->_mp_exp;
+      b1[i *  nr_colsA + k].get_mpf_t()->_mp_d = &(b[i *  nr_rowsA + k].get_mpf_t()->_mp_d[sizeB2]); 
+      b2[i *  nr_colsA + k].get_mpf_t()->_mp_d = &(b[i *  nr_rowsA + k].get_mpf_t()->_mp_d[0]); 
+    }
+  }
+ 
+  matrixMult_cuBlas(handle, tmpC1, a1, b1, d_aS, d_bS, tmpTransferLongToGMP,
+		    tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
+		    nr_rowsA, nr_colsA, nr_colsB);
+  matrixMult_cuBlas(handle, tmpC2, a2, b2, d_aS, d_bS, tmpTransferLongToGMP,
+		    tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
+		    nr_rowsA, nr_colsA, nr_colsB);
+  for(int i = 0; i < nr_colsA; ++i){
+    for(int k = 0; k < nr_rowsA; ++k){
+      a1[i *  nr_rowsA + k] += a2[i *  nr_rowsA + k]; 
+    }
+  }
+  for(int i = 0; i < nr_colsB; ++i){
+    for(int k = 0; k < nr_colsA; ++k){
+      b1[i *  nr_colsA + k] += b2[i *  nr_rowsA + k]; 
+    }
+  }
+
+  matrixMult_cuBlas(handle, C, a1, b1, d_aS, d_bS, tmpTransferLongToGMP,
+		    tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
+		    nr_rowsA, nr_colsA, nr_colsB);
+  for (int i = 0; i < nr_colsB; ++i) {
+    for (int j = 0; j < nr_rowsA; ++j) {
+      C[i *  nr_rowsA + j] -= tmpC1[i *  nr_rowsA + j];
+      C[i *  nr_rowsA + j] -= tmpC2[i *  nr_rowsA + j];
+      C[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= sizeA1;
+      tmpC2[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= size;
+      C[i * nr_rowsA + j] += (tmpC1[i * nr_rowsA + j] + tmpC2[i * nr_rowsA + j]);
+    }
+  }  
 }
 
 // Implements Toom-3 algorithm for matrix multiplciation c = a.b

      // This case corresponds to a precision of ~20000 bits in GMP
    }
    gettimeofday(&t1, NULL);
    cudaMemcpy(tmp, d_res, nr_rowsA * nr_rowsA * 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);
    addToGMPMatrixSymm(c, tmp, nr_rowsA, 2 * expA - (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;
    #pragma omp parallel for schedule(dynamic)
    for(int l = 0; l < nr_rowsA; ++l){
      for(int k = 0; k < nr_rowsA; ++k){
	tmp[l *  nr_rowsA + k] = 0;
      }
    }
  }
 
  printf("Transfer GPU = %fms\n", etTransfer);
  printf("Multiplication GPU = %fms\n", etMult);
  printf("Addition GPU = %fms\n", etAdd);

  double etBackAlloc = 0;
  double etaddBack = 0;
  
  for (int i = 0; i < size; i++) {
    gettimeofday(&t1, NULL);
    cudaMemcpy(d_res, tmp, nr_rowsA * nr_rowsA * sizeof(long long), cudaMemcpyHostToDevice);
    gettimeofday(&t2, NULL);
