 
     //Merge sub-superaccs into work-group partial-accumulator ( ATTENTION: PartialSuperacc is transposed!)
     uint pos = threadIdx.x;
-//if(pos < WARP_COUNT) {
-//        int imin = 0;
-//        int imax = 38;
-//    NormalizeT( l_workingBase, &imin, &imax);
-//}
-//    __syncthreads();
+    if(pos < WARP_COUNT) {
+            int imin = 0;
+            int imax = 38;
+        NormalizeT( l_workingBase, &imin, &imax);
+    }
+    __syncthreads();
 
     if (pos < BIN_COUNT) {
         long long int sum = 0;

    return carry_in < 0;
}

__device__
double Round( long long int *accumulator) {
    int imin = 0;
    int imax = 38;
    int negative = Normalize(accumulator, &imin, &imax);

    //Find leading word
    int i;
    //Skip zeroes
    for (i = imax; accumulator[i] == 0 && i >= imin; --i) {
    }
    if (negative) {
        //Skip ones
        for(; (accumulator[i] & ((1l << DIGITS) - 1)) == ((1l << DIGITS) - 1) && i >= imin; --i) {
        }
    }
    if (i < 0)
        return 0.0;

    long long int hiword = negative ? ((1l << DIGITS) - 1) - accumulator[i] : accumulator[i];
    double rounded = (double) hiword;
    double hi = ldexp(rounded, (i - F_WORDS) * DIGITS);
    if (i == 0)
        return negative ? -hi : hi;  // Correct rounding achieved
    hiword -= (long long int) rint(rounded);
    double mid = ldexp((double) hiword, (i - F_WORDS) * DIGITS);

    //Compute sticky
    long long int sticky = 0;
    for (int j = imin; j != i - 1; ++j)
        sticky |= negative ? (1l << DIGITS) - accumulator[j] : accumulator[j];

    long long int loword = negative ? (1l << DIGITS) - accumulator[i - 1] : accumulator[i - 1];
    loword |= !!sticky;
    double lo = ldexp((double) loword, (i - 1 - F_WORDS) * DIGITS);

    //Now add3(hi, mid, lo)
    //No overlap, we have already normalized
    if (mid != 0)
        lo = OddRoundSumNonnegative_gpu(mid, lo);

    //Final rounding
    hi = hi + lo;
    return negative ? -hi : hi;
}


////////////////////////////////////////////////////////////////////////////////
// Main computation pass: compute partial superaccs

    //Read data from global memory and scatter it to sub-superaccs
    double a[NBFPE] = {0.0};
    for(uint pos = blockIdx.x*blockDim.x+threadIdx.x; pos < NbElements; pos += gridDim.x*blockDim.x) {
        double r = 0.0;
        double x = TwoProductFMA(d_a[pos], d_b[pos], &r);

        #pragma unroll
        for(uint i = 0; i != NBFPE; ++i) {
            double s;
            a[i] = KnuthTwoSum(a[i], x, &s);
            x = s;
        }
        if (x != 0.0) {
            Accumulate(l_workingBase, x);
            // Flush FPEs to superaccs
            #pragma unroll
            for(uint i = 0; i != NBFPE; ++i) {
                Accumulate(l_workingBase, a[i]);
                a[i] = 0.0;
            }
        }

        //if (r != 0.0) {//add the rest r in the same manner
        //    #pragma unroll
