 	
 	__shared__ uint s_permutations[160];
 
-    for(int i = 0; i < 16; i++)
-    {
+	for(int i = 0; i < 16; i++)
+	{
 		int pidx = idx + NUM_THREADS * i;
 		if (pidx >= 992) break;
-        
-        ushort start, end;
-        uint permutation = permutations[pidx];
-        if (pidx < 160) s_permutations[pidx] = permutation;
-                
-        float error = evalPermutation4(colors, colorSum, permutation, &start, &end);
-        
-        if (error < bestError)
-        {
-            bestError = error;
-            bestPermutation = permutation;
-            bestStart = start;
-            bestEnd = end;
-        }
-    }
-
-    if (bestStart < bestEnd)
-    {
-        swap(bestEnd, bestStart);
-        bestPermutation ^= 0x55555555;    // Flip indices.
-    }
-
-    for(int i = 0; i < 3; i++)
-    {
+		
+		ushort start, end;
+		uint permutation = permutations[pidx];
+		if (pidx < 160) s_permutations[pidx] = permutation;
+				
+		float error = evalPermutation4(colors, colorSum, permutation, &start, &end);
+		
+		if (error < bestError)
+		{
+			bestError = error;
+			bestPermutation = permutation;
+			bestStart = start;
+			bestEnd = end;
+		}
+	}
+
+	if (bestStart < bestEnd)
+	{
+		swap(bestEnd, bestStart);
+		bestPermutation ^= 0x55555555;	// Flip indices.
+	}
+
+	for(int i = 0; i < 3; i++)
+	{
 		int pidx = idx + NUM_THREADS * i;
-        if (pidx >= 160) break;
-        
-        ushort start, end;
-        uint permutation = s_permutations[pidx];
-        float error = evalPermutation3(colors, colorSum, permutation, &start, &end);
-        
-        if (error < bestError)
-        {
-            bestError = error;
-            bestPermutation = permutation;
-            bestStart = start;
-            bestEnd = end;
-            
-            if (bestStart > bestEnd)
-            {
-                swap(bestEnd, bestStart);
-                bestPermutation ^= (~bestPermutation >> 1) & 0x55555555;    // Flip indices.
-            }
-        }
-    }
+		if (pidx >= 160) break;
+		
+		ushort start, end;
+		uint permutation = s_permutations[pidx];
+		float error = evalPermutation3(colors, colorSum, permutation, &start, &end);
+		
+		if (error < bestError)
+		{
+			bestError = error;
+			bestPermutation = permutation;
+			bestStart = start;
+			bestEnd = end;
+			
+			if (bestStart > bestEnd)
+			{
+				swap(bestEnd, bestStart);
+				bestPermutation ^= (~bestPermutation >> 1) & 0x55555555;	// Flip indices.
+			}
+		}
+	}
 
 	errors[idx] = bestError;
 }

////////////////////////////////////////////////////////////////////////////////
inline __device__ float3 roundAndExpand(float3 v, ushort * w)
{
    v.x = rintf(__saturatef(v.x) * 31.0f);
    v.y = rintf(__saturatef(v.y) * 63.0f);
    v.z = rintf(__saturatef(v.z) * 31.0f);
    *w = ((ushort)v.x << 11) | ((ushort)v.y << 5) | (ushort)v.z;
    v.x *= 0.03227752766457f; // approximate integer bit expansion.
    v.y *= 0.01583151765563f;
    v.z *= 0.03227752766457f;
    return v;
}




__device__ float evalPermutation4(const float3 * colors, float3 color_sum, uint permutation, ushort * start, ushort * end)
{
	// Compute endpoints using least squares.
	float3 alphax_sum = make_float3(0.0f, 0.0f, 0.0f);
	uint akku = 0;

	// Compute alpha & beta for this permutation.
	#pragma unroll
	for (int i = 0; i < 16; i++)
	{
		const uint bits = permutation >> (2*i);
		
		alphax_sum += alphaTable4[bits & 3] * colors[i];
		akku += prods4[bits & 3];
	}

	float alpha2_sum = float(akku >> 16);
	float beta2_sum = float((akku >> 8) & 0xff);
	float alphabeta_sum = float(akku & 0xff);
	float3 betax_sum = 9.0f * color_sum - alphax_sum;

	const float factor = 1.0f / (alpha2_sum * beta2_sum - alphabeta_sum * alphabeta_sum);

	float3 a = (alphax_sum * beta2_sum - betax_sum * alphabeta_sum) * factor;
	float3 b = (betax_sum * alpha2_sum - alphax_sum * alphabeta_sum) * factor;
	
	// Round a, b to the closest 5-6-5 color and expand...
	a = roundAndExpand(a, start);
	b = roundAndExpand(b, end);

	// compute the error
	float3 e = a * a * alpha2_sum + b * b * beta2_sum + 2.0f * (a * b * alphabeta_sum - a * alphax_sum - b * betax_sum);

	return (1.0f / 9.0f) * dot(e, kColorMetricSqr);
}

__device__ float evalPermutation3(const float3 * colors, float3 color_sum, uint permutation, ushort * start, ushort * end)
{
	// Compute endpoints using least squares.
	float3 alphax_sum = make_float3(0.0f, 0.0f, 0.0f);
	uint akku = 0;

	// Compute alpha & beta for this permutation.
	#pragma unroll
	for (int i = 0; i < 16; i++)
	{
		const uint bits = permutation >> (2*i);

		alphax_sum += alphaTable3[bits & 3] * colors[i];
		akku += prods3[bits & 3];
	}

	float alpha2_sum = float(akku >> 16);
	float beta2_sum = float((akku >> 8) & 0xff);
	float alphabeta_sum = float(akku & 0xff);
	float3 betax_sum = 4.0f * color_sum - alphax_sum;

	const float factor = 1.0f / (alpha2_sum * beta2_sum - alphabeta_sum * alphabeta_sum);

	float3 a = (alphax_sum * beta2_sum - betax_sum * alphabeta_sum) * factor;
	float3 b = (betax_sum * alpha2_sum - alphax_sum * alphabeta_sum) * factor;
	
	// Round a, b to the closest 5-6-5 color and expand...
	a = roundAndExpand(a, start);
	b = roundAndExpand(b, end);

	// compute the error
	float3 e = a * a * alpha2_sum + b * b * beta2_sum + 2.0f * (a * b * alphabeta_sum - a * alphax_sum - b * betax_sum);

	return (1.0f / 4.0f) * dot(e, kColorMetricSqr);
}

__device__ float evalPermutation4(const float3 * colors, const float * weights, float3 color_sum, uint permutation, ushort * start, ushort * end)
{
	// Compute endpoints using least squares.
	float alpha2_sum = 0.0f;
	float beta2_sum = 0.0f;
	float alphabeta_sum = 0.0f;
	float3 alphax_sum = make_float3(0.0f, 0.0f, 0.0f);

	// Compute alpha & beta for this permutation.
	for (int i = 0; i < 16; i++)
	{
		const uint bits = permutation >> (2*i);
		
		float beta = (bits & 1);
		if (bits & 2) beta = (1 + beta) / 3.0f;
		float alpha = 1.0f - beta;
		
		alpha2_sum += alpha * alpha * weights[i];
		beta2_sum += beta * beta * weights[i];
		alphabeta_sum += alpha * beta * weights[i];
		alphax_sum += alpha * colors[i];
	}

	float3 betax_sum = color_sum - alphax_sum;

	const float factor = 1.0f / (alpha2_sum * beta2_sum - alphabeta_sum * alphabeta_sum);

	float3 a = (alphax_sum * beta2_sum - betax_sum * alphabeta_sum) * factor;
	float3 b = (betax_sum * alpha2_sum - alphax_sum * alphabeta_sum) * factor;
	
	// Round a, b to the closest 5-6-5 color and expand...
	a = roundAndExpand(a, start);
	b = roundAndExpand(b, end);

	// compute the error
	float3 e = a * a * alpha2_sum + b * b * beta2_sum + 2.0f * (a * b * alphabeta_sum - a * alphax_sum - b * betax_sum);

	return dot(e, kColorMetricSqr);
}

/*
