83 		__m128 m1x,m1y,m1z;

    84 		__m128 m2x,m2y,m2z;

    85 		__m128 m;

    86 		

    87 		// m1(xyz) = u * samples[i].x

    88 		m1x = _mm_set1_ps(u.x);

    89 		m1y = _mm_set1_ps(u.y);

    90 		m1z = _mm_set1_ps(u.z);

    91 		m = _mm_set_ps(samples[0].x, samples[1].x, samples[2].x, samples[3].x);

    92 		m1x = _mm_mul_ps(m1x, m);

    93 		m1y = _mm_mul_ps(m1y, m);

    94 		m1z = _mm_mul_ps(m1z, m);

    95 		

    96 		// m2(xyz) = v * samples[i].y

    97 		m2x = _mm_set1_ps(v.x);

    98 		m2y = _mm_set1_ps(v.y);

    99 		m2z = _mm_set1_ps(v.z);

   100 		m = _mm_set_ps(samples[0].y, samples[1].y, samples[2].y, samples[3].y);

   101 		m2x = _mm_mul_ps(m2x, m);

   102 		m2y = _mm_mul_ps(m2y, m);

   103 		m2z = _mm_mul_ps(m2z, m);

   104 		

   105 		// m1(xyz) = (m1 + m2) = (u*samples[i].x + v*samples[i].y)

   106 		m1x = _mm_add_ps(m1x, m2x);

   107 		m1y = _mm_add_ps(m1y, m2y);

   108 		m1z = _mm_add_ps(m1z, m2z);

   109 		

   110 		// m1(xyz) = m1*F = (u*samples[i].x + v*samples[i].y)*F

   111 		m = _mm_set_ps(F,F,F,F);

   112 		m1x = _mm_mul_ps(m1x, m);

   113 		m1y = _mm_mul_ps(m1y, m);

   114 		m1z = _mm_mul_ps(m1z, m);

   115 		

   116 		// m1(xyz) = p - m1 = p - (u*samples[i].x + v*samples[i].y)*F = dir

   117 		m2x = _mm_set1_ps(p.x);

   118 		m2y = _mm_set1_ps(p.y);

   119 		m2z = _mm_set1_ps(p.z);

   120 		m2x = _mm_sub_ps(m2x, m1x);

   121 		m2y = _mm_sub_ps(m2y, m1y);

   122 		m2z = _mm_sub_ps(m2z, m1z);

   123 		

   124 		// normalize dir

   125 		m1x = _mm_mul_ps(m2x, m2x); // x*x

   126 		m1y = _mm_mul_ps(m2y, m2y); // y*y

   127 		m1z = _mm_mul_ps(m2z, m2z); // z*z

   128 		m = _mm_add_ps(m1x, m1y);   // x*x + y*y

   129 		m = _mm_add_ps(m, m1z);     // m = x*x + y*y + z*z

   130 		m = _mm_sqrt_ps(m);         // m = sqrt(m)

   131 		m2x = _mm_div_ps(m2x, m);   // dir(xyz) /= m

   132 		m2y = _mm_div_ps(m2y, m);

   133 		m2z = _mm_div_ps(m2z, m);

   134