static inline
void __calculate_normal( float *p1, float *p2, float *p3, float *normal, unsigned char n)
{ unsigned char i;
	float alpha;
	float *v1, *v2, *v3;
	float *u1, *u2, *u3;

	v1 = malloc(n * sizeof(float));
	v2 = malloc(n * sizeof(float));
	v3 = malloc(n * sizeof(float));
	u1 = malloc(n * sizeof(float));
	u2 = malloc(n * sizeof(float));
	u3 = malloc(n * sizeof(float));

	/*
		Calculate a normal vector of a plain using Gram-Schmidt process
	*/
	{
		for (i = 0; i < n; ++i)
		{
			v1[i] = p2[i] - p1[i];
			v2[i] = p3[i] - p1[i];
			v3[i] = p1[i];
		}

		for (i = 0; i < n; ++i)
		{
			u1[i] = v1[i];
		}

		{
			float proj[n];
			float dot_v2_u1 = 0.0f, dot_u1_u1 = 0.0f;
			for (i = 0; i < n; ++i)
			{
				dot_v2_u1 += v2[i] * u1[i];
				dot_u1_u1 += u1[i] * u1[i];
			}
			alpha = dot_v2_u1 / dot_u1_u1;

			for (i = 0; i < n; ++i)
			{
				proj[i] = u1[i] * alpha;
				u2[i] = v2[i] - proj[i];
			}
		}

		{
			float proj1[n], proj2[n];
			float dot_v3_u1 = 0.0f, dot_u1_u1 = 0.0f;
			float dot_v3_u2 = 0.0f, dot_u2_u2 = 0.0f;

			for (i = 0; i < n; ++i)
			{
				dot_v3_u1 += v3[i] * u1[i];
				dot_u1_u1 += u1[i] * u1[i];
			}
			for (i = 0; i < n; ++i)
			{
				proj1[i] = u1[i] * (dot_v3_u1 / dot_u1_u1);
			}

			for (i = 0; i < n; ++i)
			{
				dot_v3_u2 += v3[i] * u2[i];
				dot_u2_u2 += u2[i] * u2[i];
			}
			for (i = 0; i < n; ++i)
			{
				proj2[i] = u2[i] * (dot_v3_u2 / dot_u2_u2);
				u3[i] = v3[i] - proj1[i] - proj2[i];
			}
		}

		float magnitude = 0.0f;
		for (i = 0; i < n; ++i)
		{
			magnitude += u3[i] * u3[i];
		}
		magnitude = sqrtf(magnitude);

		for (i = 0; i < n; ++i)
		{
			normal[i] = u3[i] / magnitude;
		}

		free(v1);
		free(v2);
		free(v3);
		free(u1);
		free(u2);
		free(u3);
		return;
	}
}

void klein_normalize(struct klein * klein)
{
	unsigned long i;
	unsigned char j;
	float *norm_vec;

	klein->normals = malloc((klein->dim * klein->vertex_size) * sizeof(float));
	norm_vec = malloc(klein->dim * sizeof(float));

	for (i = 0; i < klein->vertex_size; i += 3 * klein->dim)
	{
		__calculate_normal(klein->vertex + i, klein->vertex + i + klein->dim, klein->vertex + i + 2 * klein->dim, norm_vec, klein->dim);
		for (j = 0; j < klein->dim; ++j )
		{
			(klein->normals + i)[j]=norm_vec[j];
			(klein->normals + i + klein->dim)[j]=norm_vec[j];
			(klein->normals + i + 2*klein->dim)[j]=norm_vec[j];
		}
	}

	free(norm_vec);
}