 
 }
 
+// CUDA kernel. Each thread takes care of one element of c
+__global__ void vecAdd(unsigned long long *a, unsigned long long *b, unsigned long long *c, int size)
+{
+  // Get our global thread ID
+  int id = blockIdx.x*blockDim.x+threadIdx.x;
+  
+  // Make sure we do not go out of bounds
+  if (id < size)
+    c[id] = a[id] + b[id];
+}
+
+
+// CUDA kernel. Each thread takes care of one element of c                                                           
+__global__ void vecAdd__wSign(double *a, long long *res, int size)
+{
+  // Get our global thread ID                                                                                       
+  int id = blockIdx.x*blockDim.x+threadIdx.x;
+  // Make sure we do not go out of bounds                                                                           
+  if (id < size)
+    res[id] += ((long long)a[id]);
+}
+
+
+
+
+// CUDA kernel. Each thread takes care of one element of c                                                          
+__global__ void vecAdd_withRem(double *a, long long *res, long long *rem,  int size)
+{
+  // Get our global thread ID                                                                                       
+  int id = blockIdx.x*blockDim.x+threadIdx.x;
+
+  // Make sure we do not go out of bounds                                                                           
+  if (id < size) {
+    if (a[id] > 0 && res[id] > LLONG_MAX - a[id]) {
+      /* handle overflow */
+      rem[id] += 1;					       
+      res[id] += (((long long)a[id]) - LLONG_MAX - 1);
+    } else if (a[id] < 0 && res[id] < LLONG_MIN - a[id]) {
+      /* handle underflow */
+      rem[id] -= 1; 
+      res[id] += (((long long)a[id]) - LLONG_MIN + 1);
+    } else {
+      /* handle regular cases */
+      res[id] += a[id];
+    }
+  }
+}
 
 void matrixSquareIntoBlock(const double *A, double * B, const int nr_rows_A, const int nr_cols_A) {
   

  }
} 

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

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