namespace VRUT
{

	inline static __m256 __forceinline RCPNR(__m256 f)
	{
		return _mm256_div_ps(_mm256_set1_ps(1.0f), f);
		//return _mm256_rcp_ps(f);
	}

	template <int OCTANT>
	static void Swizzle(__m256 & d, const __m256 & a,  const __m256 &  b)
	{
		d = _mm256_blend_ps(a, b, (OCTANT | OCTANT << 4 ));
	}

	inline static void Swizzle(__m256 & d, const __m256 & a,  const __m256 &  b, unsigned int octant)
	{
		typedef void (*SwizzleOctant)(__m256 & d, const __m256 & a,  const __m256 &  b);
		static SwizzleOctant SwizzleDispatch[] = {
			&Swizzle<0>,
			&Swizzle<1>,
			&Swizzle<2>,
			&Swizzle<3>,
			&Swizzle<4>,
			&Swizzle<5>,
			&Swizzle<6>,
			&Swizzle<7>,
		};

		SwizzleDispatch[octant](d, a, b);
	}

	inline void RayBundle::InitTraversalCache()
	{
		RayPacket & packet = *activePackets[0];
		traversalCache.ODx._p = _mm256_mul_ps(packet.DInvx._p, packet.Ox._p);
		traversalCache.ODy._p = _mm256_mul_ps(packet.DInvy._p, packet.Oy._p);
		traversalCache.ODz._p = _mm256_mul_ps(packet.DInvz._p, packet.Oz._p);
	}

	inline void RayBundle::ReorderRays(unsigned int numberOfPackets)
	{
		static const __m128i SHUFFLE_TABLE[16] = {
			_mm_setr_epi8(12,13,14,15, 8, 9,10,11, 4, 5, 6, 7, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3,12,13,14,15, 8, 9,10,11, 4, 5, 6, 7),
			_mm_setr_epi8( 4, 5, 6, 7,12,13,14,15, 8, 9,10,11, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3, 4, 5, 6, 7,12,13,14,15, 8, 9,10,11),
			_mm_setr_epi8( 8, 9,10,11,12,13,14,15, 4, 5, 6, 7, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3, 8, 9,10,11,12,13,14,15, 4, 5, 6, 7),
			_mm_setr_epi8( 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15),

			_mm_setr_epi8(12,13,14,15, 8, 9,10,11, 4, 5, 6, 7, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3,12,13,14,15, 8, 9,10,11, 4, 5, 6, 7),
			_mm_setr_epi8( 4, 5, 6, 7,12,13,14,15, 8, 9,10,11, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3, 4, 5, 6, 7,12,13,14,15, 8, 9,10,11),
			_mm_setr_epi8( 8, 9,10,11,12,13,14,15, 4, 5, 6, 7, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3, 8, 9,10,11,12,13,14,15, 4, 5, 6, 7),
			_mm_setr_epi8( 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 0, 1, 2, 3),
			_mm_setr_epi8( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15),
		};

		unsigned int rayOffsetsBytes[MAX_RAYS_IN_BUNDLE];
		unsigned int rayHits = 0;
		unsigned int rayMissed = MAX_RAYS_IN_PACKET * numberOfPackets;
		for (unsigned int packetIndex = 0; packetIndex < numberOfPackets; ++packetIndex)
		{
			unsigned int * packetOffsetBytes = activePackets[packetIndex]->rayOffsets;
			unsigned long long hits = intersections[packetIndex];
			unsigned long long hit_lo = (hits & 0x0f);
			unsigned long long hit_hi = (hits & 0xf0) >> 4;
			unsigned long long pop_lo = _mm_popcnt_u64(hit_lo);
			unsigned long long pop_hi = _mm_popcnt_u64(hit_hi);
			// load lo/hi ID pairs
			__m128i id_lo = _mm_load_si128((__m128i*)packetOffsetBytes);
			__m128i id_hi = _mm_load_si128((__m128i*)(packetOffsetBytes + 4));

			// store hit/miss iDs
			__m128i shuf_lo = _mm_shuffle_epi8(id_lo, SHUFFLE_TABLE[hit_lo]);
			__m128i shuf_hi = _mm_shuffle_epi8(id_hi, SHUFFLE_TABLE[hit_hi]);
			_mm_storeu_si128((__m128i*)&rayOffsetsBytes[rayHits], shuf_lo);
			rayHits += pop_lo;
			_mm_storeu_si128((__m128i*)&rayOffsetsBytes[rayHits], shuf_hi);
			rayHits += pop_hi;

			// NOTE: Hits MUST be written before misses, or a full-hit packet can corrupt the miss area
			_mm_storeu_si128((__m128i*)&rayOffsetsBytes[rayMissed - 4], shuf_lo);
			rayMissed -= 4 - pop_lo;
			_mm_storeu_si128((__m128i*)&rayOffsetsBytes[rayMissed - 4], shuf_hi);
			rayMissed -= 4 - pop_hi;
		}
		// read rays
		for (unsigned int packetIndex = 0; packetIndex < numberOfPackets; ++packetIndex)
			ReadRay(*activePackets[packetIndex], &rayOffsetsBytes[MAX_RAYS_IN_PACKET * packetIndex]);
	}

	inline void RayBundle::ReadDir(RayPacket & packet, Point & Dx, Point & Dy, Point & Dz)
	{
#if USE_GATHERS && defined(__AVX2__)
		__m256i idx = _mm256_load_si256((__m256i*)packet.rayOffsets);
		Dx._p = _mm256_i32gather_ps(&rays->direction.x, idx, 1);
		Dy._p = _mm256_i32gather_ps(&rays->direction.y, idx, 1);
		Dz._p = _mm256_i32gather_ps(&rays->direction.z, idx, 1);
#else
#ifndef __WXMSW__
#define _mm256_set_m128(va, vb) \
	_mm256_insertf128_ps(_mm256_castps128_ps256(vb), va, 1)
#endif

		char * bytes = ((char*)rays) + 16;
#define LOADPS(x) _mm_load_ps((float*)(x))
		__m256 v0 = _mm256_set_m128(LOADPS(bytes + packet.rayOffsets[4]), LOADPS(bytes + packet.rayOffsets[0]));  // 0000 4444
		__m256 v1 = _mm256_set_m128(LOADPS(bytes + packet.rayOffsets[6]), LOADPS(bytes + packet.rayOffsets[2]));  // 2222 6666
		__m256 v2 = _mm256_set_m128(LOADPS(bytes + packet.rayOffsets[5]), LOADPS(bytes + packet.rayOffsets[1]));  // 1111 5555
		__m256 v3 = _mm256_set_m128(LOADPS(bytes + packet.rayOffsets[7]), LOADPS(bytes + packet.rayOffsets[3]));  // 3333 7777
#undef LOADPS
		__m256 t0 = _mm256_unpacklo_ps(v0,v1); // 02 02 46 46
		__m256 t1 = _mm256_unpackhi_ps(v0,v1); // 02 02 46 46
		__m256 t2 = _mm256_unpacklo_ps(v2,v3); // 13 13 57 57
		__m256 t3 = _mm256_unpackhi_ps(v2,v3); // 13 13 57 57
		Dx._p  = _mm256_unpacklo_ps(t0,t2);  // 01 23 45 67
		Dy._p  = _mm256_unpackhi_ps(t0,t2);
		Dz._p  = _mm256_unpacklo_ps(t1,t3);
#endif
	}

	inline void RayBundle::ReadRay(RayPacket & packet, const unsigned int * packetOffsetBytes)
	{
#if USE_GATHERS && defined(__AVX2__)
		__m256i idx = _mm256_loadu_si256((__m256i*)packetOffsetBytes);
		__m256 Ox = _mm256_i32gather_ps(&rays->origin.x, idx, 1);
		__m256 Oy = _mm256_i32gather_ps(&rays->origin.y, idx, 1);
		__m256 Oz = _mm256_i32gather_ps(&rays->origin.z, idx, 1);
		__m256 tmax = _mm256_i32gather_ps(&rays->t, idx, 1 );
		__m256 Dx = _mm256_i32gather_ps(&rays->direction.x, idx, 1);
		__m256 Dy = _mm256_i32gather_ps(&rays->direction.y, idx, 1);
		__m256 Dz = _mm256_i32gather_ps(&rays->direction.z, idx, 1);
		__m256 RID = _mm256_i32gather_ps((float*)&(rays->offset), idx,1);
		packet.Ox._p = Ox;
		packet.Oy._p = Oy;
		packet.Oz._p = Oz;
		packet.DInvx._p = RCPNR(Dx);
		packet.DInvy._p = RCPNR(Dy);
		packet.DInvz._p = RCPNR(Dz);
		packet.t._p = tmax;
		_mm256_store_ps((float*)(packet.rayOffsets), RID);
#else
		// unpacklo(x,y) -->   x0 y0 x1 y1 x4 y4 x5 y5
		// 0 1 2 3   0 1 2 3
		//  
		//   0 1 2 3   0 1 2 3
		//   0 1 2 3   0 1 2 3
		//   0 1 2 3   0 1 2 3
		//   0 1 2 3   0 1 2 3

		char * bytes = ((char*)rays);
#define LOADPS(x) _mm_load_ps((float*)(x))
		__m256 l0 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[1]), LOADPS(bytes + packetOffsetBytes[0]));
		__m256 l1 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[3]), LOADPS(bytes + packetOffsetBytes[2]));
		__m256 l2 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[5]), LOADPS(bytes + packetOffsetBytes[4]));
		__m256 l3 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[7]), LOADPS(bytes + packetOffsetBytes[6]));
		__m256 l4 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[1] + 16), LOADPS(bytes + packetOffsetBytes[0] + 16));
		__m256 l5 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[3] + 16), LOADPS(bytes + packetOffsetBytes[2] + 16));
		__m256 l6 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[5] + 16), LOADPS(bytes + packetOffsetBytes[4] + 16));
		__m256 l7 = _mm256_set_m128(LOADPS(bytes + packetOffsetBytes[7] + 16), LOADPS(bytes + packetOffsetBytes[6] + 16));
#undef LOADPS

		__m256 t0 = _mm256_unpacklo_ps(l0,l1); // 00 11 00 11
		__m256 t1 = _mm256_unpacklo_ps(l2,l3); // 00 11 00 11
		__m256 t2 = _mm256_unpackhi_ps(l0,l1); // 22 33 22 33
		__m256 t3 = _mm256_unpackhi_ps(l2,l3); // 22 33 22 33
		__m256 t4 = _mm256_unpacklo_ps(l4,l5); // 44 55 44 55
		__m256 t5 = _mm256_unpacklo_ps(l6,l7); // 44 55 44 55
		__m256 t6 = _mm256_unpackhi_ps(l4,l5); // 66 77 66 77 
		__m256 t7 = _mm256_unpackhi_ps(l6,l7); // 66 77 66 77

		__m256 Ox = _mm256_unpacklo_ps(t0,t1); // 00 00 00 00
		__m256 Oy = _mm256_unpackhi_ps(t0,t1); // 11 11 11 11
		__m256 Oz = _mm256_unpacklo_ps(t2,t3); // 22 22 22 22
		__m256 RID = _mm256_unpackhi_ps(t2,t3); // 33 33 33 33
		__m256 Dx = _mm256_unpacklo_ps(t4,t5); 
		__m256 Dy = _mm256_unpackhi_ps(t4,t5); 
		__m256 Dz = _mm256_unpacklo_ps(t6,t7);
		__m256 tmax = _mm256_unpackhi_ps(t6,t7);

		_mm256_store_ps((float*)&packet.Ox, Ox);
		_mm256_store_ps((float*)&packet.DInvx, RCPNR(Dx));
		_mm256_store_ps((float*)&packet.Oy, Oy);    
		_mm256_store_ps((float*)&packet.DInvy, RCPNR(Dy));
		_mm256_store_ps((float*)&packet.Oz, Oz);    
		_mm256_store_ps((float*)&packet.DInvz, RCPNR(Dz));
		_mm256_store_ps((float*)&packet.t, tmax );
		_mm256_store_ps((float*)packet.rayOffsets, RID);

		// unpacklo(x,y) -->   x0 y0 x1 y1 x4 y4 x5 y5
		// 0 1 2 3   0 1 2 3
		//  
		//   0 1 2 3   0 1 2 3
		//   0 1 2 3   0 1 2 3
		//   0 1 2 3   0 1 2 3
		//   0 1 2 3   0 1 2 3
#endif
	}

	inline void RayBundle::TransformRayUnit(unsigned int packetIndex, float t[], Point & tOx, Point & tOy, Point& tOz, Point & tDx, Point & tDy, Point& tDz)
	{
		RayPacket & p = *activePackets[packetIndex];
		memcpy(t, p.t.v, sizeof(float) * MAX_RAYS_IN_PACKET);
		tOx = p.Ox;
		tOy = p.Oy;
		tOz = p.Oz;
		ReadDir(p, tDx, tDy, tDz);
	}

	inline void RayBundle::TransformRay(const Matrix & matrix, unsigned int packetIndex, float t[], Point & tOx, Point & tOy, Point& tOz, Point & tDx, Point & tDy, Point& tDz)
	{
		RayPacket & p = *activePackets[packetIndex];
		memcpy(t, p.t.v, sizeof(float) * MAX_RAYS_IN_PACKET);
		RayBundle::Point Dx, Dy, Dz;
		ReadDir(p, Dx, Dy, Dz);
#ifdef __AVX2__
		__m256 tmp = _mm256_broadcast_ss(&matrix._m11);
		tOx._p = _mm256_mul_ps(p.Ox._p, tmp);
		tDx._p = _mm256_mul_ps(Dx._p, tmp);
		tmp = _mm256_broadcast_ss(&matrix._m21);
		tOx._p = _mm256_fmadd_ps(p.Oy._p, tmp, tOx._p);
		tDx._p = _mm256_fmadd_ps(Dy._p, tmp, tDx._p);
		tmp = _mm256_broadcast_ss(&matrix._m31);
		tOx._p = _mm256_fmadd_ps(p.Oz._p, tmp, tOx._p);
		tDx._p = _mm256_fmadd_ps(Dz._p, tmp, tDx._p);
		tmp = _mm256_broadcast_ss(&matrix._m41);
		tOx._p = _mm256_add_ps(tOx._p, tmp);

		tmp = _mm256_broadcast_ss(&matrix._m12);
		tOy._p = _mm256_mul_ps(p.Ox._p, tmp);
		tDy._p = _mm256_mul_ps(Dx._p, tmp);
		tmp = _mm256_broadcast_ss(&matrix._m22);
		tOy._p = _mm256_fmadd_ps(p.Oy._p, tmp, tOy._p);
		tDy._p = _mm256_fmadd_ps(Dy._p, tmp, tDy._p);
		tmp = _mm256_broadcast_ss(&matrix._m32);
		tOy._p = _mm256_fmadd_ps(p.Oz._p, tmp, tOy._p);
		tDy._p = _mm256_fmadd_ps(Dz._p, tmp, tDy._p);
		tmp = _mm256_broadcast_ss(&matrix._m42);
		tOy._p = _mm256_add_ps(tOy._p, tmp);

		tmp = _mm256_broadcast_ss(&matrix._m13);
		tOz._p = _mm256_mul_ps(p.Ox._p, tmp);
		tDz._p = _mm256_mul_ps(Dx._p, tmp);
		tmp = _mm256_broadcast_ss(&matrix._m23);
		tOz._p = _mm256_fmadd_ps(p.Oy._p, tmp, tOz._p);
		tDz._p = _mm256_fmadd_ps(Dy._p, tmp, tDz._p);
		tmp = _mm256_broadcast_ss(&matrix._m33);
		tOz._p = _mm256_fmadd_ps(p.Oz._p, tmp, tOz._p);
		tDz._p = _mm256_fmadd_ps(Dz._p, tmp, tDz._p);
		tmp = _mm256_broadcast_ss(&matrix._m43);
		tOz._p = _mm256_add_ps(tOz._p, tmp);
#else
		__m256 tmp = _mm256_broadcast_ss(&matrix._m11);
		tOx._p = _mm256_mul_ps(p.Ox._p, tmp);
		tDx._p = _mm256_mul_ps(Dx._p, tmp);
		tmp = _mm256_broadcast_ss(&matrix._m21);
		tOx._p = _mm256_add_ps(tOx._p, _mm256_mul_ps(p.Oy._p, tmp));
		tDx._p = _mm256_add_ps(tDx._p, _mm256_mul_ps(Dy._p, tmp));
		tmp = _mm256_broadcast_ss(&matrix._m31);
		tOx._p = _mm256_add_ps(tOx._p, _mm256_mul_ps(p.Oz._p, tmp));
		tDx._p = _mm256_add_ps(tDx._p, _mm256_mul_ps(Dz._p, tmp));
		tmp = _mm256_broadcast_ss(&matrix._m41);
		tOx._p = _mm256_add_ps(tOx._p, tmp);

		tmp = _mm256_broadcast_ss(&matrix._m12);
		tOy._p = _mm256_mul_ps(p.Ox._p, tmp);
		tDy._p = _mm256_mul_ps(Dx._p, tmp);
		tmp = _mm256_broadcast_ss(&matrix._m22);
		tOy._p = _mm256_add_ps(tOy._p, _mm256_mul_ps(p.Oy._p, tmp));
		tDy._p = _mm256_add_ps(tDy._p, _mm256_mul_ps(Dy._p, tmp));
		tmp = _mm256_broadcast_ss(&matrix._m32);
		tOy._p = _mm256_add_ps(tOy._p, _mm256_mul_ps(p.Oz._p, tmp));
		tDy._p = _mm256_add_ps(tDy._p, _mm256_mul_ps(Dz._p, tmp));
		tmp = _mm256_broadcast_ss(&matrix._m42);
		tOy._p = _mm256_add_ps(tOy._p, tmp);

		tmp = _mm256_broadcast_ss(&matrix._m13);
		tOz._p = _mm256_mul_ps(p.Ox._p, tmp);
		tDz._p = _mm256_mul_ps(Dx._p, tmp);
		tmp = _mm256_broadcast_ss(&matrix._m23);
		tOz._p = _mm256_add_ps(tOz._p, _mm256_mul_ps(p.Oy._p, tmp));
		tDz._p = _mm256_add_ps(tDz._p, _mm256_mul_ps(Dy._p, tmp));
		tmp = _mm256_broadcast_ss(&matrix._m33);
		tOz._p = _mm256_add_ps(tOz._p, _mm256_mul_ps(p.Oz._p, tmp));
		tDz._p = _mm256_add_ps(tDz._p, _mm256_mul_ps(Dz._p, tmp));
		tmp = _mm256_broadcast_ss(&matrix._m43);
		tOz._p = _mm256_add_ps(tOz._p, tmp);
#endif
	}

	inline unsigned int RayBundle::IntersectAABBCull(const float * fbox, unsigned int numberOfPackets)
	{
		int frontPacketIndex, backPacketIndex;
		unsigned int hitRays = 0;

		for (frontPacketIndex = 0; frontPacketIndex < int(numberOfPackets); ++frontPacketIndex)
		{
			RayPacket & packet = *activePackets[frontPacketIndex];
			__m256 Bmin0  = _mm256_broadcast_ss(fbox + 0);
			__m256 Bmax0  = _mm256_broadcast_ss(fbox + 3);
			__m256 tmin0 = _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Ox._p), packet.DInvx._p);
			__m256 tmax0 = _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Ox._p), packet.DInvx._p);

			Bmin0 = _mm256_broadcast_ss(fbox + 1);
			Bmax0 = _mm256_broadcast_ss(fbox + 4);
			tmin0 = _mm256_max_ps(tmin0, _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Oy._p), packet.DInvy._p));
			tmax0 = _mm256_min_ps(tmax0, _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Oy._p), packet.DInvy._p));

			Bmin0 = _mm256_broadcast_ss(fbox + 2);
			Bmax0 = _mm256_broadcast_ss(fbox + 5);
			tmin0 = _mm256_max_ps(tmin0, _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Oz._p), packet.DInvz._p));
			tmax0 = _mm256_min_ps(tmax0, _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Oz._p), packet.DInvz._p));

			__m256 l0  = _mm256_min_ps(tmax0, packet.t._p);
			__m256 hit = _mm256_andnot_ps(tmax0, _mm256_cmp_ps(tmin0, l0, _CMP_LE_OQ));

			unsigned int mask = (unsigned int)_mm256_movemask_ps(hit);
			hitRays += _mm_popcnt_u64(mask);
			intersections[frontPacketIndex] = (unsigned char)mask;

			if (mask)
				break;
		}

		for (backPacketIndex = numberOfPackets - 1; backPacketIndex > frontPacketIndex; --backPacketIndex)
		{
			RayPacket & packet = *activePackets[backPacketIndex];
			__m256 Bmin0  = _mm256_broadcast_ss(fbox + 0);
			__m256 Bmax0  = _mm256_broadcast_ss(fbox + 3);
			__m256 tmin0 = _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Ox._p), packet.DInvx._p);
			__m256 tmax0 = _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Ox._p), packet.DInvx._p);

			Bmin0 = _mm256_broadcast_ss(fbox + 1);
			Bmax0 = _mm256_broadcast_ss(fbox + 4);
			tmin0 = _mm256_max_ps(tmin0, _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Oy._p), packet.DInvy._p));
			tmax0 = _mm256_min_ps(tmax0, _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Oy._p), packet.DInvy._p));

			Bmin0 = _mm256_broadcast_ss(fbox + 2);
			Bmax0 = _mm256_broadcast_ss(fbox + 5);
			tmin0 = _mm256_max_ps(tmin0, _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Oz._p), packet.DInvz._p));
			tmax0 = _mm256_min_ps(tmax0, _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Oz._p), packet.DInvz._p));

			__m256 l0  = _mm256_min_ps(tmax0, packet.t._p);
			__m256 hit = _mm256_andnot_ps(tmax0, _mm256_cmp_ps(tmin0, l0, _CMP_LE_OQ));

			unsigned int mask = (unsigned int)_mm256_movemask_ps(hit);
			hitRays += _mm_popcnt_u64(mask);
			intersections[backPacketIndex] = (unsigned char)mask;

			if (mask)
				break;
		}

		int remaining = backPacketIndex - frontPacketIndex - 1;
		if (remaining > 0)
		{
			memset(intersections + frontPacketIndex + 1, 0xff, remaining);
			hitRays += MAX_RAYS_IN_PACKET * remaining;
		}

		return hitRays;
	}

	inline unsigned int RayBundle::IntersectAABB(const float * fbox, unsigned int numberOfPackets)
	{
		unsigned int packetIndex;
		unsigned int hitRays = 0;
		unsigned int numberOfPackets2 = (numberOfPackets & ~1); 
		for (packetIndex = 0; packetIndex < numberOfPackets2; packetIndex += 2)
		{
			RayPacket & packet0 = *activePackets[packetIndex + 0];
			RayPacket & packet1 = *activePackets[packetIndex + 1];
#ifdef __AVX2__
			__m256 Bmin = _mm256_broadcast_ss(fbox + 0);
			__m256 Bmax = _mm256_broadcast_ss(fbox + 3);
			__m256 O0 = _mm256_load_ps((float*)&packet0.Ox);
			__m256 O1 = _mm256_load_ps((float*)&packet1.Ox);
			__m256 D0 = _mm256_load_ps((float*)&packet0.DInvx);
			__m256 D1 = _mm256_load_ps((float*)&packet1.DInvx);
			__m256 r0 = _mm256_mul_ps(O0, D0);
			__m256 r1 = _mm256_mul_ps(O1, D1);
			__m256 tmin0 = _mm256_fmsub_ps(Bmin, D0, r0);
			__m256 tmax0 = _mm256_fmsub_ps(Bmax, D0, r0);
			__m256 tmin1 = _mm256_fmsub_ps(Bmin, D1, r1);
			__m256 tmax1 = _mm256_fmsub_ps(Bmax, D1, r1);

			Bmin = _mm256_broadcast_ss(fbox + 1);
			Bmax = _mm256_broadcast_ss(fbox + 4);
			O0 = _mm256_load_ps((float*)&packet0.Oy);
			O1 = _mm256_load_ps((float*)&packet1.Oy);
			D0 = _mm256_load_ps((float*)&packet0.DInvy);
			D1 = _mm256_load_ps((float*)&packet1.DInvy);
			r0 = _mm256_mul_ps(O0, D0);
			r1 = _mm256_mul_ps(O1, D1);
			__m256 t0 = _mm256_fmsub_ps(Bmin, D0, r0);
			__m256 t1 = _mm256_fmsub_ps(Bmax, D0, r0);
			__m256 t2 = _mm256_fmsub_ps(Bmin, D1, r1);
			__m256 t3 = _mm256_fmsub_ps(Bmax, D1, r1);
			tmin0 = _mm256_max_ps(tmin0, t0);
			tmax0 = _mm256_min_ps(tmax0, t1);
			tmin1 = _mm256_max_ps(tmin1, t2);
			tmax1 = _mm256_min_ps(tmax1, t3);

			Bmin = _mm256_broadcast_ss(fbox + 2);        
			Bmax = _mm256_broadcast_ss(fbox + 5); 
			O0 = _mm256_load_ps((float*)&packet0.Oz);
			O1 = _mm256_load_ps((float*)&packet1.Oz);
			D0 = _mm256_load_ps((float*)&packet0.DInvz);
			D1 = _mm256_load_ps((float*)&packet1.DInvz);
			r0 = _mm256_mul_ps(O0, D0);
			r1 = _mm256_mul_ps(O1, D1);
			t0 = _mm256_fmsub_ps(Bmin, D0, r0);
			t1 = _mm256_fmsub_ps(Bmax, D0, r0);
			t2 = _mm256_fmsub_ps(Bmin, D1, r1);
			t3 = _mm256_fmsub_ps(Bmax, D1, r1);
			tmin0 = _mm256_max_ps(tmin0, t0);
			tmax0 = _mm256_min_ps(tmax0, t1);
			tmin1 = _mm256_max_ps(tmin1, t2);
			tmax1 = _mm256_min_ps(tmax1, t3);
#else
			__m256 Bmin = _mm256_broadcast_ss(fbox + 0);
			__m256 Bmax = _mm256_broadcast_ss(fbox + 3);
			__m256 O0 = _mm256_load_ps((float*)&packet0.Ox);
			__m256 O1 = _mm256_load_ps((float*)&packet1.Ox);
			__m256 D0 = _mm256_load_ps((float*)&packet0.DInvx);
			__m256 D1 = _mm256_load_ps((float*)&packet1.DInvx);
			__m256 tmin0 = _mm256_mul_ps(_mm256_sub_ps(Bmin, O0), D0);
			__m256 tmax0 = _mm256_mul_ps(_mm256_sub_ps(Bmax, O0), D0);
			__m256 tmin1 = _mm256_mul_ps(_mm256_sub_ps(Bmin, O1), D1);
			__m256 tmax1 = _mm256_mul_ps(_mm256_sub_ps(Bmax, O1), D1);

			Bmin = _mm256_broadcast_ss(fbox + 1);
			Bmax = _mm256_broadcast_ss(fbox + 4);
			O0 = _mm256_load_ps((float*)&packet0.Oy);
			O1 = _mm256_load_ps((float*)&packet1.Oy);
			D0 = _mm256_load_ps((float*)&packet0.DInvy);
			D1 = _mm256_load_ps((float*)&packet1.DInvy);
			__m256 t0 = _mm256_mul_ps(_mm256_sub_ps(Bmin, O0), D0);
			__m256 t1 = _mm256_mul_ps(_mm256_sub_ps(Bmax, O0), D0);
			__m256 t2 = _mm256_mul_ps(_mm256_sub_ps(Bmin, O1), D1);
			__m256 t3 = _mm256_mul_ps(_mm256_sub_ps(Bmax, O1), D1);
			tmin0 = _mm256_max_ps(tmin0, t0);
			tmax0 = _mm256_min_ps(tmax0, t1);
			tmin1 = _mm256_max_ps(tmin1, t2);
			tmax1 = _mm256_min_ps(tmax1, t3);

			Bmin = _mm256_broadcast_ss(fbox + 2);        
			Bmax = _mm256_broadcast_ss(fbox + 5); 
			O0 = _mm256_load_ps((float*)&packet0.Oz);
			O1 = _mm256_load_ps((float*)&packet1.Oz);
			D0 = _mm256_load_ps((float*)&packet0.DInvz);
			D1 = _mm256_load_ps((float*)&packet1.DInvz);
			t0 = _mm256_mul_ps(_mm256_sub_ps(Bmin, O0), D0);
			t1 = _mm256_mul_ps(_mm256_sub_ps(Bmax, O0), D0);
			t2 = _mm256_mul_ps(_mm256_sub_ps(Bmin, O1), D1);
			t3 = _mm256_mul_ps(_mm256_sub_ps(Bmax, O1), D1);
			tmin0 = _mm256_max_ps(tmin0, t0);
			tmax0 = _mm256_min_ps(tmax0, t1);
			tmin1 = _mm256_max_ps(tmin1, t2);
			tmax1 = _mm256_min_ps(tmax1, t3);
#endif

			// andnot -> uses sign-bit trick for tmax>=0
			__m256 l0 = _mm256_min_ps(tmax0, packet0.t._p);
			__m256 l1 = _mm256_min_ps(tmax1, packet1.t._p);
			__m256 hit0 = _mm256_andnot_ps(tmax0, _mm256_cmp_ps(tmin0, l0, _CMP_LE_OQ));
			__m256 hit1 = _mm256_andnot_ps(tmax1, _mm256_cmp_ps(tmin1, l1, _CMP_LE_OQ));

			unsigned int mask0 = (unsigned int)_mm256_movemask_ps(hit0);
			unsigned int mask1 = (unsigned int)_mm256_movemask_ps(hit1);
			unsigned int nMergedMask = (mask1 << 8) | mask0;

			*((unsigned short int*)(intersections + packetIndex)) = (unsigned short int)nMergedMask;
			hitRays += _mm_popcnt_u64(nMergedMask);        
		}

		if (packetIndex < numberOfPackets)
		{
			RayPacket & packet = *activePackets[packetIndex];
			__m256 Bmin0  = _mm256_broadcast_ss(fbox + 0);
			__m256 Bmax0  = _mm256_broadcast_ss(fbox + 3);
			__m256 tmin0 = _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Ox._p), packet.DInvx._p);
			__m256 tmax0 = _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Ox._p), packet.DInvx._p);

			Bmin0 = _mm256_broadcast_ss(fbox + 1);
			Bmax0 = _mm256_broadcast_ss(fbox + 4);
			tmin0 = _mm256_max_ps(tmin0, _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Oy._p), packet.DInvy._p));
			tmax0 = _mm256_min_ps(tmax0, _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Oy._p), packet.DInvy._p));

			Bmin0 = _mm256_broadcast_ss(fbox + 2);
			Bmax0 = _mm256_broadcast_ss(fbox + 5);
			tmin0 = _mm256_max_ps(tmin0, _mm256_mul_ps(_mm256_sub_ps(Bmin0, packet.Oz._p), packet.DInvz._p));
			tmax0 = _mm256_min_ps(tmax0, _mm256_mul_ps(_mm256_sub_ps(Bmax0, packet.Oz._p), packet.DInvz._p));

			__m256 l0  = _mm256_min_ps(tmax0, packet.t._p);
			__m256 hit = _mm256_andnot_ps(tmax0, _mm256_cmp_ps(tmin0, l0, _CMP_LE_OQ));

			unsigned int mask = (unsigned int)_mm256_movemask_ps(hit);
			hitRays += _mm_popcnt_u64(mask);
			intersections[packetIndex] = (unsigned char)mask;
		}

		return hitRays;
	}

	inline unsigned int RayBundle::IntersectAABB(const AbstractBVH::BVHNode * nodes, unsigned int nodeIndex, unsigned int octant)
	{
		RayPacket & packet = *activePackets[0];

		const AbstractBVH::BVHNode * leftNode = nodes + nodeIndex;
		const AbstractBVH::BVHNode * rightNode = leftNode + 1;

		const float * fLeftBox  = (const float*)&leftNode->box;
		const float * fRightBox = (const float*)&rightNode->box;

		// fetch AABB.  each of the min/max vectors has the same data replicated into low and high lanes
		// bbmin(xyz), bbmax(xyz)
		__m256 min0 = _mm256_broadcast_ps((const __m128*)(fLeftBox));
		__m256 min1 = _mm256_broadcast_ps((const __m128*)(fRightBox));
		__m256 max0 = _mm256_broadcast_ps((const __m128*)(fLeftBox + 3));
		__m256 max1 = _mm256_broadcast_ps((const __m128*)(fRightBox + 3));

		// swap planes based on octant
		__m256 bbmin0, bbmax0, bbmin1, bbmax1;
		Swizzle(bbmin0, min0, max0, octant);
		Swizzle(bbmax0, max0, min0, octant);
		Swizzle(bbmin1, min1, max1, octant);
		Swizzle(bbmax1, max1, min1, octant);

		// prefetch left/right children since we'll probably need them
		_mm_prefetch((const char*)(nodes + leftNode->t), _MM_HINT_T0);
		_mm_prefetch((const char*)(nodes + rightNode->t), _MM_HINT_T0);

		// axis tests
#ifdef __AVX2__
		// TODO: test AVX2
		__m256 D = _mm256_load_ps((float*)&packet.DInvx);
		__m256 O = _mm256_load_ps((float*)&traversalCache.ODx);
		__m256 Bmin0 = _mm256_permute_ps(bbmin0, 0x00);
		__m256 Bmax0 = _mm256_permute_ps(bbmax0, 0x00);
		__m256 Bmin1 = _mm256_permute_ps(bbmin1, 0x00);
		__m256 Bmax1 = _mm256_permute_ps(bbmax1, 0x00);
		__m256 tmin0 = _mm256_fmsub_ps(Bmin0, D, O);
		__m256 tmax0 = _mm256_fmsub_ps(Bmax0, D, O);
		__m256 tmin1 = _mm256_fmsub_ps(Bmin1, D, O);
		__m256 tmax1 = _mm256_fmsub_ps(Bmax1, D, O);

		D = _mm256_load_ps((float*)&packet.DInvy);
		O = _mm256_load_ps((float*)&traversalCache.ODy);
		Bmin0 = _mm256_permute_ps(bbmin0, 0x55); // 0101
		Bmax0 = _mm256_permute_ps(bbmax0, 0x55);
		Bmin1 = _mm256_permute_ps(bbmin1, 0x55);
		Bmax1 = _mm256_permute_ps(bbmax1, 0x55);
		__m256 t0 = _mm256_fmsub_ps(Bmin0, D, O);
		__m256 t1 = _mm256_fmsub_ps(Bmax0, D, O);
		__m256 t2 = _mm256_fmsub_ps(Bmin1, D, O);
		__m256 t3 = _mm256_fmsub_ps(Bmax1, D, O);
		tmin0 = _mm256_max_ps(tmin0, t0);
		tmax0 = _mm256_min_ps(tmax0, t1);
		tmin1 = _mm256_max_ps(tmin1, t2);
		tmax1 = _mm256_min_ps(tmax1, t3);

		D = _mm256_load_ps((float*)&packet.DInvz);
		O = _mm256_load_ps((float*)&traversalCache.ODz);
		Bmin0 = _mm256_permute_ps(bbmin0, 0xAA); // 1010
		Bmax0 = _mm256_permute_ps(bbmax0, 0xAA);
		Bmin1 = _mm256_permute_ps(bbmin1, 0xAA);
		Bmax1 = _mm256_permute_ps(bbmax1, 0xAA);
		t0 = _mm256_fmsub_ps(Bmin0, D, O);
		t1 = _mm256_fmsub_ps(Bmax0, D, O);
		t2 = _mm256_fmsub_ps(Bmin1, D, O);
		t3 = _mm256_fmsub_ps(Bmax1, D, O);
		tmin0 = _mm256_max_ps(tmin0, t0);
		tmax0 = _mm256_min_ps(tmax0, t1);
		tmin1 = _mm256_max_ps(tmin1, t2);
		tmax1 = _mm256_min_ps(tmax1, t3);
#else
		__m256 D = _mm256_load_ps((float*)&packet.DInvx);
		__m256 O = _mm256_load_ps((float*)&packet.Ox);
		__m256 Bmin0 = _mm256_permute_ps(bbmin0, 0x00);
		__m256 Bmax0 = _mm256_permute_ps(bbmax0, 0x00);
		__m256 Bmin1 = _mm256_permute_ps(bbmin1, 0x00);
		__m256 Bmax1 = _mm256_permute_ps(bbmax1, 0x00);
		__m256 tmin0 = _mm256_mul_ps(_mm256_sub_ps(Bmin0, O), D);
		__m256 tmax0 = _mm256_mul_ps(_mm256_sub_ps(Bmax0, O), D);
		__m256 tmin1 = _mm256_mul_ps(_mm256_sub_ps(Bmin1, O), D);
		__m256 tmax1 = _mm256_mul_ps(_mm256_sub_ps(Bmax1, O), D);

		D = _mm256_load_ps((float*)&packet.DInvy);
		O = _mm256_load_ps((float*)&packet.Oy);
		Bmin0 = _mm256_permute_ps(bbmin0, 0x55); // 0101
		Bmax0 = _mm256_permute_ps(bbmax0, 0x55);
		Bmin1 = _mm256_permute_ps(bbmin1, 0x55);
		Bmax1 = _mm256_permute_ps(bbmax1, 0x55);
		__m256 t0 = _mm256_mul_ps(_mm256_sub_ps(Bmin0, O), D);
		__m256 t1 = _mm256_mul_ps(_mm256_sub_ps(Bmax0, O), D);
		__m256 t2 = _mm256_mul_ps(_mm256_sub_ps(Bmin1, O), D);
		__m256 t3 = _mm256_mul_ps(_mm256_sub_ps(Bmax1, O), D);
		tmin0 = _mm256_max_ps(tmin0, t0);
		tmax0 = _mm256_min_ps(tmax0, t1);
		tmin1 = _mm256_max_ps(tmin1, t2);
		tmax1 = _mm256_min_ps(tmax1, t3);

		D = _mm256_load_ps((float*)&packet.DInvz);
		O = _mm256_load_ps((float*)&packet.Oz);
		Bmin0 = _mm256_permute_ps(bbmin0, 0xAA); // 1010
		Bmax0 = _mm256_permute_ps(bbmax0, 0xAA);
		Bmin1 = _mm256_permute_ps(bbmin1, 0xAA);
		Bmax1 = _mm256_permute_ps(bbmax1, 0xAA);
		t0 = _mm256_mul_ps(_mm256_sub_ps(Bmin0, O), D);
		t1 = _mm256_mul_ps(_mm256_sub_ps(Bmax0, O), D);
		t2 = _mm256_mul_ps(_mm256_sub_ps(Bmin1, O), D);
		t3 = _mm256_mul_ps(_mm256_sub_ps(Bmax1, O), D);
		tmin0 = _mm256_max_ps(tmin0, t0);
		tmax0 = _mm256_min_ps(tmax0, t1);
		tmin1 = _mm256_max_ps(tmin1, t2);
		tmax1 = _mm256_min_ps(tmax1, t3);
#endif

		__m256 limL = _mm256_min_ps(tmax0, packet.t._p);
		__m256 limR = _mm256_min_ps(tmax1, packet.t._p);

		// using sign-bit trick for tmax >= 0
		t0 = _mm256_cmp_ps(tmin0, limL, _CMP_LE_OQ);
		t1 = _mm256_cmp_ps(tmin1, limR, _CMP_LE_OQ);
		__m256 hitL = _mm256_andnot_ps(tmax0, t0);
		__m256 hitR = _mm256_andnot_ps(tmax1, t1);

		unsigned int maskhitL = _mm256_movemask_ps(hitL);
		unsigned int maskhitR = _mm256_movemask_ps(hitR);
		unsigned int maskhitB = maskhitL & maskhitR;
		unsigned int orderFlags = 0x00;

		if (maskhitB)
		{
			__m256 LFirst = _mm256_cmp_ps(tmin0, tmin1, _CMP_LT_OQ);
			unsigned int lf = _mm256_movemask_ps(LFirst) & maskhitB;
			unsigned int rf = ~lf & maskhitB;
			if(_mm_popcnt_u64(lf) > _mm_popcnt_u64(rf))
				orderFlags = 0x03; // LR
			else
				orderFlags = 0x06; // RL
		}
		else
		{
			if(maskhitL)
				orderFlags |= 0x01; // L
			if(maskhitR)
				orderFlags |= 0x02; // R
		}

		return orderFlags;
	}

	inline void RayBundle::IntersectTriangles(const TriangleCol * triangleCols, int triangleBegin, int triangleEnd, unsigned int numberOfPackets)
	{
		for (unsigned int packetIndex = 0; packetIndex < numberOfPackets; ++packetIndex)
		{

			RayPacket & packet = *activePackets[packetIndex];
			__m256 ha = _mm256_setzero_ps();
			__m256 hb = _mm256_setzero_ps();
			__m256 h = _mm256_setzero_ps();
			__m256 id = _mm256_setzero_ps();
			__m256 tmax = packet.t._p;

			Point D[3];
			ReadDir(packet, D[0], D[1], D[2]);

			for (int triangleIndex = triangleBegin; triangleIndex < triangleEnd; ++triangleIndex)
			{

				const TriangleCol * triangle = &triangleCols[triangleIndex];
				const float * P0 = (const float*)&triangle->v1;
				const float * v10 = (const float*)&triangle->e21;
				const float * v02 = (const float*)&triangle->e13;
				const float * v10x02 = (const float*)&triangle->e21x13;

				_mm_prefetch((char*)(triangle + 1), _MM_HINT_T0);

				__m256 vDx = D[0]._p;
				__m256 vDy = D[1]._p;
				__m256 vDz = D[2]._p;

#ifdef __AVX2__
				__m256 v0A[3] = {
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[0]), packet.Ox._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[1]), packet.Oy._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[2]), packet.Oz._p),
				};
				__m256 v0AxD[] = {
					_mm256_fmsub_ps(v0A[1], vDz, _mm256_mul_ps(v0A[2], vDy)),
					_mm256_fmsub_ps(v0A[2], vDx, _mm256_mul_ps(v0A[0], vDz)),
					_mm256_fmsub_ps(v0A[0], vDy, _mm256_mul_ps(v0A[1], vDx))
				};
#else
				__m256 v0A[3] = {
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[0]), packet.Ox._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[1]), packet.Oy._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[2]), packet.Oz._p),
				};
				__m256 v0AxD[] = {
					_mm256_sub_ps(_mm256_mul_ps(v0A[1], vDz), _mm256_mul_ps(v0A[2], vDy)),
					_mm256_sub_ps(_mm256_mul_ps(v0A[2], vDx), _mm256_mul_ps(v0A[0], vDz)),
					_mm256_sub_ps(_mm256_mul_ps(v0A[0], vDy), _mm256_mul_ps(v0A[1], vDx))
				};
#endif

				//V = ((p1 - p0)x(p0 -p2)).d
				//Va = ((p1 - p0)x(p0 -p2)).(p0 -a)
				//V1 = ((p0 - a)×d).(p0 -p2)
				//V2 = ((p0 - a)×d).(p1 -p0)

				__m256 v10x02_x = _mm256_broadcast_ss(&v10x02[0]);
				__m256 v10x02_y = _mm256_broadcast_ss(&v10x02[1]);
				__m256 v10x02_z = _mm256_broadcast_ss(&v10x02[2]);

#ifdef __AVX2__
				__m256 t = _mm256_fmadd_ps(v10x02_z, v0A[2], _mm256_fmadd_ps(v10x02_y, v0A[1], _mm256_mul_ps(v10x02_x, v0A[0])));
				__m256 v = _mm256_fmadd_ps(v10x02_z, vDz, _mm256_fmadd_ps(v10x02_y, vDy, _mm256_mul_ps(v10x02_x, vDx)));
#else
				__m256 t = _mm256_add_ps(_mm256_mul_ps(v10x02_z, v0A[2]), _mm256_add_ps(_mm256_mul_ps(v10x02_y, v0A[1]), _mm256_mul_ps(v10x02_x, v0A[0])));
				__m256 v = _mm256_add_ps(_mm256_mul_ps(v10x02_z, vDz), _mm256_add_ps(_mm256_mul_ps(v10x02_y, vDy), _mm256_mul_ps(v10x02_x, vDx)));
#endif

				v = RCPNR(v);

				__m256 v02_x = _mm256_broadcast_ss(&v02[0]);
				__m256 v02_y = _mm256_broadcast_ss(&v02[1]);
				__m256 v02_z = _mm256_broadcast_ss(&v02[2]);

#ifdef __AVX2__
				__m256 a = _mm256_fmadd_ps(v0AxD[2], v02_z, _mm256_fmadd_ps(v0AxD[1], v02_y, _mm256_mul_ps(v0AxD[0],  v02_x)));
#else
				__m256 a = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v02_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1], v02_y), _mm256_mul_ps(v0AxD[0],  v02_x)));
#endif

				__m256 v10_x = _mm256_broadcast_ss(&v10[0]);
				__m256 v10_y = _mm256_broadcast_ss(&v10[1]);
				__m256 v10_z = _mm256_broadcast_ss(&v10[2]);

#ifdef __AVX2__
				__m256 b = _mm256_fmadd_ps(v0AxD[2], v10_z, _mm256_fmadd_ps(v0AxD[1], v10_y, _mm256_mul_ps(v0AxD[0], v10_x)));
#else
				__m256 b = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v10_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1],  v10_y), _mm256_mul_ps(v0AxD[0],  v10_x)));
#endif

				a = _mm256_mul_ps(a, v);
				b = _mm256_mul_ps(b, v);
				t = _mm256_mul_ps(t, v);

				__m256 front = _mm256_and_ps(_mm256_cmp_ps(t, _mm256_setzero_ps(), _CMP_GT_OQ), _mm256_cmp_ps(tmax, t, _CMP_GT_OQ));
				__m256 in = _mm256_and_ps(_mm256_cmp_ps(a, _mm256_setzero_ps(), _CMP_GE_OQ ), _mm256_cmp_ps(b, _mm256_setzero_ps(), _CMP_GE_OQ));

				in = _mm256_and_ps(in, _mm256_cmp_ps(_mm256_add_ps(a, b), _mm256_set1_ps(1.0f), _CMP_LE_OQ ));

				__m256 hit = _mm256_and_ps(in, front);

				ha = _mm256_blendv_ps(ha, a, hit);
				hb = _mm256_blendv_ps(hb, b, hit);
				tmax = _mm256_blendv_ps(tmax, t, hit);
				h = _mm256_or_ps(h, hit);
				id = _mm256_blendv_ps(id, _mm256_broadcast_ss((float*)&triangleIndex), hit);

			}

			// now write results back out
			unsigned int nHit = _mm256_movemask_ps(h);
			if(nHit)
			{
				// merge T values back into packet
				_mm256_store_ps((float*)&packet.t , tmax);

				// scatter out results
				__m256 u = ha;
				__m256 v = hb;

				VRUT_ALIGNED(32) float pu[MAX_RAYS_IN_PACKET];
				VRUT_ALIGNED(32) float pv[MAX_RAYS_IN_PACKET];
				VRUT_ALIGNED(32) float pt[MAX_RAYS_IN_PACKET];
				VRUT_ALIGNED(32) unsigned int pid[MAX_RAYS_IN_PACKET];

				_mm256_store_ps((float*)pu, u);
				_mm256_store_ps((float*)pv, v);
				_mm256_store_ps((float*)pt, tmax);
				_mm256_store_ps((float*)pid, id);

				for (int i = 0; i < MAX_RAYS_IN_PACKET; ++i)
				{
					if ((1 << i) & nHit)
					{
						unsigned int id = packet.rayOffsets[i] / sizeof(Ray);
						rays[id].t = pt[i];
						infos[id].dist = pt[i];
						infos[id].u = pu[i];
						infos[id].v = pv[i];
						infos[id].index = pid[i];
					}
				}
			}
		}
	}

	inline void RayBundle::IntersectTriangles(const TriangleCol * triangleCols, int triangleBegin, int triangleEnd)
	{
		RayPacket & packet = *activePackets[0];
		__m256 tmax = packet.t._p;

		for (int triangleIndex = triangleBegin; triangleIndex < triangleEnd; ++triangleIndex)
		{

			const TriangleCol * triangle = &triangleCols[triangleIndex];
			const float * P0 = (const float*)&triangle->v1;
			const float * v10 = (const float*)&triangle->e21;
			const float * v02 = (const float*)&triangle->e13;
			const float * v10x02 = (const float*)&triangle->e21x13;

			_mm_prefetch((char*)(triangle + 1), _MM_HINT_T0);

			__m256 vDx = _mm256_load_ps(intersectionCache.Dx.v);
			__m256 vDy = _mm256_load_ps(intersectionCache.Dy.v);
			__m256 vDz = _mm256_load_ps(intersectionCache.Dz.v);
			__m256 vOx = _mm256_load_ps(packet.Ox.v);
			__m256 vOy = _mm256_load_ps(packet.Oy.v);
			__m256 vOz = _mm256_load_ps(packet.Oz.v);

			__m256 v0A[3] = {
				_mm256_sub_ps(_mm256_broadcast_ss(&P0[0]), vOx),
				_mm256_sub_ps(_mm256_broadcast_ss(&P0[1]), vOy),
				_mm256_sub_ps(_mm256_broadcast_ss(&P0[2]), vOz),
			};

			//V = ((p1 - p0)x(p0 -p2)).d
			//Va = ((p1 - p0)x(p0 -p2)).(p0 -a)
			//V1 = ((p0 - a)×d).(p0 -p2)
			//V2 = ((p0 - a)×d).(p1 -p0)

			__m256 v10x02_x = _mm256_broadcast_ss(&v10x02[0]);
			__m256 v10x02_y = _mm256_broadcast_ss(&v10x02[1]);
			__m256 v10x02_z = _mm256_broadcast_ss(&v10x02[2]);

#ifdef __AVX2__
			__m256 t = _mm256_fmadd_ps(v10x02_z, v0A[2], _mm256_fmadd_ps(v10x02_y, v0A[1], _mm256_mul_ps(v10x02_x, v0A[0])));
			__m256 v = _mm256_fmadd_ps(v10x02_z, vDz, _mm256_fmadd_ps(v10x02_y, vDy, _mm256_mul_ps(v10x02_x, vDx)));
#else
			__m256 t = _mm256_add_ps(_mm256_mul_ps(v10x02_z, v0A[2]), _mm256_add_ps(_mm256_mul_ps(v10x02_y, v0A[1]), _mm256_mul_ps(v10x02_x, v0A[0])));
			__m256 v = _mm256_add_ps(_mm256_mul_ps(v10x02_z, vDz), _mm256_add_ps(_mm256_mul_ps(v10x02_y, vDy), _mm256_mul_ps(v10x02_x, vDx)));
#endif

			v = RCPNR(v);

			__m256 v02_x = _mm256_broadcast_ss(&v02[0]);
			__m256 v02_y = _mm256_broadcast_ss(&v02[1]);
			__m256 v02_z = _mm256_broadcast_ss(&v02[2]);

#ifdef __AVX2__
			__m256 v0AxD[] = {
				_mm256_fmsub_ps( v0A[1], vDz, _mm256_mul_ps( v0A[2], vDy ) ),
				_mm256_fmsub_ps( v0A[2], vDx, _mm256_mul_ps( v0A[0], vDz ) ),
				_mm256_fmsub_ps( v0A[0], vDy, _mm256_mul_ps( v0A[1], vDx ) )
			};
			__m256 a = _mm256_fmadd_ps(v0AxD[2], v02_z, _mm256_fmadd_ps(v0AxD[1],  v02_y, _mm256_mul_ps(v0AxD[0], v02_x)));
#else
			__m256 v0AxD[] = {
				_mm256_sub_ps(_mm256_mul_ps(v0A[1], vDz), _mm256_mul_ps(v0A[2], vDy)),
				_mm256_sub_ps(_mm256_mul_ps(v0A[2], vDx), _mm256_mul_ps(v0A[0], vDz)),
				_mm256_sub_ps(_mm256_mul_ps(v0A[0], vDy), _mm256_mul_ps(v0A[1], vDx))
			};
			__m256 a = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v02_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1], v02_y), _mm256_mul_ps(v0AxD[0], v02_x)));
#endif

			__m256 v10_x = _mm256_broadcast_ss(&v10[0]);
			__m256 v10_y = _mm256_broadcast_ss(&v10[1]);
			__m256 v10_z = _mm256_broadcast_ss(&v10[2]);

#ifdef __AVX2__
			__m256 b = _mm256_fmadd_ps(v0AxD[2], v10_z, _mm256_fmadd_ps( v0AxD[1],  v10_y, _mm256_mul_ps(v0AxD[0], v10_x)));
#else
			__m256 b = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v10_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1],  v10_y), _mm256_mul_ps(v0AxD[0], v10_x)));
#endif

			a = _mm256_mul_ps(a, v);
			b = _mm256_mul_ps(b, v);
			t = _mm256_mul_ps(t, v);

			__m256 front = _mm256_and_ps(_mm256_cmp_ps(t ,_mm256_setzero_ps(), _CMP_GT_OQ), _mm256_cmp_ps(tmax, t, _CMP_GT_OQ));
			__m256 in = _mm256_and_ps(_mm256_cmp_ps(a, _mm256_setzero_ps(),_CMP_GE_OQ), _mm256_cmp_ps(b, _mm256_setzero_ps(),_CMP_GE_OQ));

			in = _mm256_and_ps(in, _mm256_cmp_ps(_mm256_add_ps(a, b), _mm256_set1_ps(1.0f), _CMP_LE_OQ));

			__m256 hit = _mm256_and_ps(in, front);

			unsigned int nHit = _mm256_movemask_ps(hit);
			if (nHit)
			{
				tmax = _mm256_blendv_ps(tmax, t, hit);
				_mm256_store_ps((float*)&packet.t._p, tmax);
				__m256 stored_a = _mm256_load_ps(intersectionCache.ha.v);
				__m256 stored_b = _mm256_load_ps(intersectionCache.hb.v);
				__m256 stored_id = _mm256_load_ps((float*)intersectionCache.id);
				__m256 id = _mm256_broadcast_ss((float*)&triangleIndex);
				stored_a = _mm256_blendv_ps(stored_a, a, hit);
				stored_b = _mm256_blendv_ps(stored_b, b, hit);
				stored_id = _mm256_blendv_ps(stored_id, id, hit);
				_mm256_store_ps(intersectionCache.ha.v, stored_a);
				_mm256_store_ps(intersectionCache.hb.v, stored_b);
				_mm256_store_ps((float*)intersectionCache.id, stored_id);
				intersectionCache.mask |= nHit;
			}

		}

	}

	inline unsigned int RayBundle::IntersectTrianglesShadow(const TriangleCol * triangleCols, int triangleBegin, int triangleEnd, unsigned int numberOfPackets, unsigned int numberOfAllPackets)
	{
		unsigned int unoccludedPackets = 0;
		for (unsigned int packetIndex = 0; packetIndex < numberOfPackets; ++packetIndex)
		{
			RayPacket & packet = *activePackets[packetIndex];

			unsigned int occluded = 0;
			__m256 h = _mm256_setzero_ps();
			__m256 tmax = packet.t._p;

			Point D[3];
			ReadDir(packet, D[0], D[1], D[2]);

			for (int triangleIndex = triangleBegin; triangleIndex < triangleEnd; ++triangleIndex)
			{

				const TriangleCol * triangle = &triangleCols[triangleIndex];
				const float * P0 = (const float*)&triangle->v1;
				const float * v10 = (const float*)&triangle->e21;
				const float * v02 = (const float*)&triangle->e13;
				const float * v10x02 = (const float*)&triangle->e21x13;

				_mm_prefetch((char*)(triangle + 1), _MM_HINT_T0);

				__m256 vDx = D[0]._p;
				__m256 vDy = D[1]._p;
				__m256 vDz = D[2]._p;

#ifdef __AVX2__
				__m256 v0A[3] = {
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[0]), packet.Ox._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[1]), packet.Oy._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[2]), packet.Oz._p),
				};
				__m256 v0AxD[] = {
					_mm256_fmsub_ps(v0A[1], vDz, _mm256_mul_ps(v0A[2], vDy)),
					_mm256_fmsub_ps(v0A[2], vDx, _mm256_mul_ps(v0A[0], vDz)),
					_mm256_fmsub_ps(v0A[0], vDy, _mm256_mul_ps(v0A[1], vDx))
				};
#else
				__m256 v0A[3] = {
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[0]), packet.Ox._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[1]), packet.Oy._p),
					_mm256_sub_ps(_mm256_broadcast_ss(&P0[2]), packet.Oz._p),
				};
				__m256 v0AxD[] = {
					_mm256_sub_ps(_mm256_mul_ps(v0A[1], vDz), _mm256_mul_ps(v0A[2], vDy)),
					_mm256_sub_ps(_mm256_mul_ps(v0A[2], vDx), _mm256_mul_ps(v0A[0], vDz)),
					_mm256_sub_ps(_mm256_mul_ps(v0A[0], vDy), _mm256_mul_ps(v0A[1], vDx))
				};
#endif

				//V = ((p1 - p0)x(p0 -p2)).d
				//Va = ((p1 - p0)x(p0 -p2)).(p0 -a)
				//V1 = ((p0 - a)×d).(p0 -p2)
				//V2 = ((p0 - a)×d).(p1 -p0)

				__m256 v10x02_x = _mm256_broadcast_ss(&v10x02[0]);
				__m256 v10x02_y = _mm256_broadcast_ss(&v10x02[1]);
				__m256 v10x02_z = _mm256_broadcast_ss(&v10x02[2]);

#ifdef __AVX2__
				__m256 t = _mm256_fmadd_ps(v10x02_z, v0A[2], _mm256_fmadd_ps(v10x02_y, v0A[1], _mm256_mul_ps(v10x02_x, v0A[0])));
				__m256 v = _mm256_fmadd_ps(v10x02_z, vDz, _mm256_fmadd_ps(v10x02_y, vDy, _mm256_mul_ps(v10x02_x, vDx)));
#else
				__m256 t = _mm256_add_ps(_mm256_mul_ps(v10x02_z, v0A[2]), _mm256_add_ps(_mm256_mul_ps(v10x02_y, v0A[1]), _mm256_mul_ps(v10x02_x, v0A[0])));
				__m256 v = _mm256_add_ps(_mm256_mul_ps(v10x02_z, vDz), _mm256_add_ps(_mm256_mul_ps(v10x02_y, vDy), _mm256_mul_ps(v10x02_x, vDx)));
#endif

				v = RCPNR(v);

				__m256 v02_x = _mm256_broadcast_ss(&v02[0]);
				__m256 v02_y = _mm256_broadcast_ss(&v02[1]);
				__m256 v02_z = _mm256_broadcast_ss(&v02[2]);

#ifdef __AVX2__
				__m256 a = _mm256_fmadd_ps(v0AxD[2], v02_z, _mm256_fmadd_ps(v0AxD[1], v02_y, _mm256_mul_ps(v0AxD[0],  v02_x)));
#else
				__m256 a = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v02_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1], v02_y), _mm256_mul_ps(v0AxD[0],  v02_x)));
#endif

				__m256 v10_x = _mm256_broadcast_ss(&v10[0]);
				__m256 v10_y = _mm256_broadcast_ss(&v10[1]);
				__m256 v10_z = _mm256_broadcast_ss(&v10[2]);

#ifdef __AVX2__
				__m256 b = _mm256_fmadd_ps(v0AxD[2], v10_z, _mm256_fmadd_ps(v0AxD[1], v10_y, _mm256_mul_ps(v0AxD[0], v10_x)));
#else
				__m256 b = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v10_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1],  v10_y), _mm256_mul_ps(v0AxD[0],  v10_x)));
#endif

				a = _mm256_mul_ps(a, v);
				b = _mm256_mul_ps(b, v);
				t = _mm256_mul_ps(t, v);

				__m256 front = _mm256_and_ps(_mm256_cmp_ps(t, _mm256_setzero_ps(), _CMP_GT_OQ), _mm256_cmp_ps(tmax, t, _CMP_GT_OQ));
				__m256 in = _mm256_and_ps(_mm256_cmp_ps(a, _mm256_setzero_ps(), _CMP_GE_OQ ), _mm256_cmp_ps(b, _mm256_setzero_ps(), _CMP_GE_OQ));

				in = _mm256_and_ps(in, _mm256_cmp_ps(_mm256_add_ps(a, b), _mm256_set1_ps(1.0f), _CMP_LE_OQ));

				__m256 hit = _mm256_and_ps(in, front);

				h = _mm256_or_ps(h, hit);

				occluded |= _mm256_movemask_ps(h);
				if (occluded == FULL_HIT)
					break;

			}

			// force tmax=-1 for hit rays so that all subsequent intersection tests on them will miss
			__m256 MINUS_ONE = _mm256_set1_ps(-1.0f);
			tmax = _mm256_blendv_ps(tmax, MINUS_ONE, h);
			_mm256_store_ps((float*)&packet.t , tmax);

			// mark hits in output array and set t=0 in ray structs
			for (int i = 0; i < MAX_RAYS_IN_PACKET; ++i)
			{
				if ((1 << i) & occluded)
				{
					unsigned int id = packet.rayOffsets[i] / sizeof(Ray);
					rays[id].t = 0.0f;
					infos[id].dist = 0.0f;
				}
			}

			if (occluded < FULL_HIT)
				activePackets[unoccludedPackets++] = activePackets[packetIndex];

		}

		return unoccludedPackets;
	}

	inline void RayBundle::IntersectTrianglesShadow(const TriangleCol * triangleCols, int triangleBegin, int triangleEnd)
	{
		RayPacket & packet = *activePackets[0];
		__m256 tmax = packet.t._p;
		__m256 h = _mm256_setzero_ps();
		unsigned int occluded = intersectionCache.mask;

		for (int triangleIndex = triangleBegin; triangleIndex < triangleEnd; ++triangleIndex)
		{

			const TriangleCol * triangle = &triangleCols[triangleIndex];
			const float * P0 = (const float*)&triangle->v1;
			const float * v10 = (const float*)&triangle->e21;
			const float * v02 = (const float*)&triangle->e13;
			const float * v10x02 = (const float*)&triangle->e21x13;

			_mm_prefetch((char*)(triangle + 1), _MM_HINT_T0);

			__m256 vDx = _mm256_load_ps(intersectionCache.Dx.v);
			__m256 vDy = _mm256_load_ps(intersectionCache.Dy.v);
			__m256 vDz = _mm256_load_ps(intersectionCache.Dz.v);
			__m256 vOx = _mm256_load_ps(packet.Ox.v);
			__m256 vOy = _mm256_load_ps(packet.Oy.v);
			__m256 vOz = _mm256_load_ps(packet.Oz.v);

			__m256 v0A[3] = {
				_mm256_sub_ps(_mm256_broadcast_ss(&P0[0]), vOx),
				_mm256_sub_ps(_mm256_broadcast_ss(&P0[1]), vOy),
				_mm256_sub_ps(_mm256_broadcast_ss(&P0[2]), vOz),
			};

			//V = ((p1 - p0)x(p0 -p2)).d
			//Va = ((p1 - p0)x(p0 -p2)).(p0 -a)
			//V1 = ((p0 - a)×d).(p0 -p2)
			//V2 = ((p0 - a)×d).(p1 -p0)

			__m256 v10x02_x = _mm256_broadcast_ss(&v10x02[0]);
			__m256 v10x02_y = _mm256_broadcast_ss(&v10x02[1]);
			__m256 v10x02_z = _mm256_broadcast_ss(&v10x02[2]);

#ifdef __AVX2__
			__m256 t = _mm256_fmadd_ps(v10x02_z, v0A[2], _mm256_fmadd_ps(v10x02_y, v0A[1], _mm256_mul_ps(v10x02_x, v0A[0])));
			__m256 v = _mm256_fmadd_ps(v10x02_z, vDz, _mm256_fmadd_ps(v10x02_y, vDy, _mm256_mul_ps(v10x02_x, vDx)));
#else
			__m256 t = _mm256_add_ps(_mm256_mul_ps(v10x02_z, v0A[2]), _mm256_add_ps(_mm256_mul_ps(v10x02_y, v0A[1]), _mm256_mul_ps(v10x02_x, v0A[0])));
			__m256 v = _mm256_add_ps(_mm256_mul_ps(v10x02_z, vDz), _mm256_add_ps(_mm256_mul_ps(v10x02_y, vDy), _mm256_mul_ps(v10x02_x, vDx)));
#endif

			v = RCPNR(v);

			__m256 v02_x = _mm256_broadcast_ss(&v02[0]);
			__m256 v02_y = _mm256_broadcast_ss(&v02[1]);
			__m256 v02_z = _mm256_broadcast_ss(&v02[2]);

#ifdef __AVX2__
			__m256 v0AxD[] = {
				_mm256_fmsub_ps( v0A[1], vDz, _mm256_mul_ps( v0A[2], vDy ) ),
				_mm256_fmsub_ps( v0A[2], vDx, _mm256_mul_ps( v0A[0], vDz ) ),
				_mm256_fmsub_ps( v0A[0], vDy, _mm256_mul_ps( v0A[1], vDx ) )
			};
			__m256 a = _mm256_fmadd_ps(v0AxD[2], v02_z, _mm256_fmadd_ps(v0AxD[1],  v02_y, _mm256_mul_ps(v0AxD[0], v02_x)));
#else
			__m256 v0AxD[] = {
				_mm256_sub_ps(_mm256_mul_ps(v0A[1], vDz), _mm256_mul_ps(v0A[2], vDy)),
				_mm256_sub_ps(_mm256_mul_ps(v0A[2], vDx), _mm256_mul_ps(v0A[0], vDz)),
				_mm256_sub_ps(_mm256_mul_ps(v0A[0], vDy), _mm256_mul_ps(v0A[1], vDx))
			};
			__m256 a = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v02_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1], v02_y), _mm256_mul_ps(v0AxD[0], v02_x)));
#endif

			__m256 v10_x = _mm256_broadcast_ss(&v10[0]);
			__m256 v10_y = _mm256_broadcast_ss(&v10[1]);
			__m256 v10_z = _mm256_broadcast_ss(&v10[2]);

#ifdef __AVX2__
			__m256 b = _mm256_fmadd_ps(v0AxD[2], v10_z, _mm256_fmadd_ps( v0AxD[1],  v10_y, _mm256_mul_ps(v0AxD[0], v10_x)));
#else
			__m256 b = _mm256_add_ps(_mm256_mul_ps(v0AxD[2], v10_z), _mm256_add_ps(_mm256_mul_ps(v0AxD[1],  v10_y), _mm256_mul_ps(v0AxD[0], v10_x)));
#endif

			a = _mm256_mul_ps(a, v);
			b = _mm256_mul_ps(b, v);
			t = _mm256_mul_ps(t, v);

			__m256 front = _mm256_and_ps(_mm256_cmp_ps(t ,_mm256_setzero_ps(), _CMP_GT_OQ), _mm256_cmp_ps(tmax, t, _CMP_GT_OQ));
			__m256 in = _mm256_and_ps(_mm256_cmp_ps(a, _mm256_setzero_ps(),_CMP_GE_OQ), _mm256_cmp_ps(b, _mm256_setzero_ps(),_CMP_GE_OQ));

			in = _mm256_and_ps(in, _mm256_cmp_ps(_mm256_add_ps(a, b), _mm256_set1_ps(1.0f), _CMP_LE_OQ));

			__m256 hit = _mm256_and_ps(in, front);

			h = _mm256_or_ps(h,hit);

			occluded |= _mm256_movemask_ps(h);
			if (occluded == FULL_HIT)
				break;

		}

		// force tmax=-1 for hit rays so that all subsequent intersection tests on them will miss
		__m256 MINUS_ONE = _mm256_set1_ps(-1.0f);
		tmax = _mm256_blendv_ps(tmax, MINUS_ONE, h);
		_mm256_store_ps((float*)&packet.t , tmax);
		intersectionCache.mask = occluded;

	}

}	// namespace VRUT
