 	}
 }
 
-__global__ void bp_weight_f(float d_weight[10][36][4][4], float d_preact[10], float p_output[36][4][4])
-{
-	FOR4(10, 36, 4, 4, d_weight[i1][i2][i3][i4] = d_preact[i1] * p_output[i2][i3][i4]);
+__global__ void bp_weight_f(float d_weight[10][6][6][6], float d_preact[10], float p_output[6][6][6])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 10*6*6*6;
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 10);
+		const int i2 = ((idx /= 10	) % 6);
+		const int i3 = ((idx /= 6	) % 6);
+		const int i4 = ((idx /= 6	) % 6);
+
+		d_weight[i1][i2][i3][i4] = d_preact[i1] * p_output[i2][i3][i4];
+	}
 }
 
 __global__ void bp_bias_f(float bias[10], float d_preact[10])
 {
-	FOR1(10, bias[i1] += dt * d_preact[i1]);
-}
-
-__global__ void bp_output_s2(float d_output[36][4][4], float n_weight[10][36][4][4], float nd_preact[10])
-{
-	FOR4(10, 36, 4, 4, atomicAdd(&d_output[i2][i3][i4], n_weight[i1][i2][i3][i4] * nd_preact[i1]));
-}
-
-template<int I1, int I2, int I3, int I4>
-__global__ void bp_output_s(float d_output[I1][I2+I3-1][I2+I3-1], float n_weight[I4][I2][I2], float nd_preact[I1*I4][I3][I3])
-{
-	FOR6(I1, I2, I2, I3, I3, I4,
-		atomicAdd(&d_output[i1][i2+i4-1][i3+i5-1], n_weight[i6][i2][i3] * nd_preact[i6 * I1 + i1][i4][i5])
-	);
-}
-
-template<int I1, int I2>
-__global__ void bp_preact_s(float d_preact[I1][I2][I2], float d_output[I1][I2][I2], float preact[I1][I2][I2])
-{
-	FOR3(I1, I2, I2, 
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 10;
+
+	for (int idx = N * pos / size; idx < N * (pos+1) / size; ++idx) {
+		bias[idx] += dt * d_preact[idx];
+	}
+}
+
+__global__ void bp_output_s1(float d_output[6][6][6], float n_weight[10][6][6][6], float nd_preact[10])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 10*6*6*6;
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 10);
+		const int i2 = ((idx /= 10	) % 6);
+		const int i3 = ((idx /= 6	) % 6);
+		const int i4 = ((idx /= 6	) % 6);
+
+		atomicAdd(&d_output[i2][i3][i4], n_weight[i1][i2][i3][i4] * nd_preact[i1]);
+	}
+}
+
+__global__ void bp_preact_s1(float d_preact[6][6][6], float d_output[6][6][6], float preact[6][6][6])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 6*6*6;
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 6);
+		const int i2 = ((idx /= 6	) % 6);
+		const int i3 = ((idx /= 6	) % 6);
+
 		const float o = step_function(preact[i1][i2][i3]);
 
 		d_preact[i1][i2][i3] = d_output[i1][i2][i3] * o * (1 - o);
-	);
-}
-
-template<int I1, int I2, int I3>
-__global__ void bp_weight_s(float d_weight[I1][I1], float d_preact[I2][I3][I3], float p_output[I2][I1*I3][I1*I3])
-{
-	const float d = I2 * I3 * I3;
-
-	FOR5(I1, I1, I2, I3, I3,
-		atomicAdd(&d_weight[i1][i2], d_preact[i3][i4][i5] * p_output[i3][i4 * I1 + i1][i5 * I1 + i2]) / d);
-}
-
-template<int I1, int I2>
-__global__ void bp_bias_s(float bias[1], float d_preact[I1][I2][I2])
-{
-	const float d = I1 * I2 * I2;
-
-	FOR3(I1, I2, I2, atomicAdd(&bias[0], dt * d_preact[i1][i2][i3] / d));
-}
-
-template<int I1, int I2, int I3>
-__global__ void bp_output_c(float d_output[I2][I1*I3][I1*I3], float n_weight[I1][I1], float nd_preact[I2][I3][I3])
-{
-	FOR5(I1, I1, I2, I3, I3,
-		atomicAdd(&d_output[i3][i4 * I1 + i1][i5 * I1 + i2], n_weight[i1][i2] * nd_preact[i3][i4][i5]));
-}
-
-template<int I1, int I2>
-__global__ void bp_preact_c(float d_preact[I1][I2][I2], float d_output[I1][I2][I2], float preact[I1][I2][I2])
-{
-	FOR3(I1, I2, I2, 
+	}
+}
+
+__global__ void bp_weight_s1(float d_weight[1][4][4], float d_preact[6][6][6], float p_output[6][24][24])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 1*4*4*6*6*6;
+	const float d = pow(6.0f, 3.0f);
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 1);
+		const int i2 = ((idx /= 1	) % 4);
+		const int i3 = ((idx /= 4	) % 4);
+		const int i4 = ((idx /= 4	) % 6);
+		const int i5 = ((idx /= 6	) % 6);
+		const int i6 = ((idx /= 6	) % 6);
+
+		atomicAdd(&d_weight[i1][i2][i3], d_preact[i4][i5][i6] * p_output[i4][i5 * 4 + i2][i6 * 4 + i3]);
+	}
+}
+
+__global__ void bp_bias_s1(float bias[1], float d_preact[6][6][6])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 6*6*6;
+	const float d = pow(6.0f, 3.0f);
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 6);
+		const int i2 = ((idx /= 6	) % 6);
+		const int i3 = ((idx /= 6	) % 6);
+
+		atomicAdd(&bias[0], dt * d_preact[i1][i2][i3] / d);
+	}
+}
+
+__global__ void bp_output_c1(float d_output[6][24][24], float n_weight[1][4][4], float nd_preact[6][6][6])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 1*4*4*6*6*6;
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 1);
+		const int i2 = ((idx /= 1	) % 4);
+		const int i3 = ((idx /= 4	) % 4);
+		const int i4 = ((idx /= 4	) % 6);
+		const int i5 = ((idx /= 6	) % 6);
+		const int i6 = ((idx /= 6	) % 6);
+
+		atomicAdd(&d_output[i4][i5 * 4 + i2][i6 * 4 + i3], n_weight[i1][i2][i3] * nd_preact[i4][i5][i6]);
+	}
+}
+
+__global__ void bp_preact_c1(float d_preact[6][24][24], float d_output[6][24][24], float preact[6][24][24])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 6*24*24;
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 6);
+		const int i2 = ((idx /= 6	) % 24);
+		const int i3 = ((idx /= 24	) % 24);
+
 		const float o = step_function(preact[i1][i2][i3]);
 
 		d_preact[i1][i2][i3] = d_output[i1][i2][i3] * o * (1 - o);
-	);
-}
-
-template<int I1, int I2, int I3, int I4>
-__global__ void bp_weight_c(float d_weight[I1][I2][I2], float d_preact[I1*I4][I3][I3],
-	float p_output[I4][I3 + I2 - 1][I3 + I2 - 1])
-{
-	const float d = I3 * I3 * I4;
-
-	FOR6(I1, I2, I2, I3, I3, I4,
-		atomicAdd(&d_weight[i1][i2][i3], d_preact[i6 * I1 + i1][i4][i5] * p_output[i6][i4 + i2][i5 + i3] / d));
-}
-
-template<int I1, int I2, int I3>
-__global__ void bp_bias_c(float bias[I1], float d_preact[I1 * I3][I2][I2])
-{
-	const float d = I2 * I2 * I3;
-
-	FOR4(I1, I2, I2, I3, atomicAdd(&bias[i1], dt * d_preact[i1 * I3 + i4][i2][i3] / d));
+	}
+}
+
+__global__ void bp_weight_c1(float d_weight[6][5][5], float d_preact[6][24][24], float p_output[28][28])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 6*5*5*24*24;
+	const float d = pow(24.0f, 2.0f);
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 6);
+		const int i2 = ((idx /= 6	) % 5);
+		const int i3 = ((idx /= 5	) % 5);
+		const int i4 = ((idx /= 5	) % 24);
+		const int i5 = ((idx /= 24	) % 24);
+
+		atomicAdd(&d_weight[i1][i2][i3], d_preact[i1][i4][i5] * p_output[i4 + i2][i5 + i3] / d);
+	}
+}
+
+__global__ void bp_bias_c1(float bias[6], float d_preact[6][24][24])
+{
+	const int pos = blockIdx.x * blockDim.x + threadIdx.x;
+	const int size = blockDim.x * gridDim.x;
+
+	const int N = 6*24*24;
+	const float d = pow(24.0f, 2.0f);
+
+	for (int n = N * pos / size; n < N * (pos+1) / size; ++n) {
+		int idx = n;
+		const int i1 = ((idx /= 1	) % 6);
+		const int i2 = ((idx /= 6	) % 24);
+		const int i3 = ((idx /= 24	) % 24);
+
+		atomicAdd(&bias[i1], dt * d_preact[i1][i2][i3] / d);
+	}
 }
 
 static void propagate(double data[28][28])

	cublasCreate(&blas);

	float err;
	int iter = 3;

	while (iter < 0 || iter-- > 0) {
		err = 0.0f;

#include <cublas_v2.h>

const static float dt = 1.0E-01f;
const static float threshold = 1.0E-02f;

struct Layer {
	float *output;
