/********************************************************
Perlin Noise2 and Noise3 functions Optimised for 
Intel SIMD (PIII CeleronII) Instruction Set.

SIMD versions process 4 input vectors at once.

SIMD versions run more than twice as fast :)

This is my first attempt at using SIMD so constructive
comments welcome - any one from Intel reading ? 

Please excuse the gross use of prefetches! I just threw
them in there (they ARE needed)!

I'm SURE these are sub-optimal.

Compiled with Intel's excellent Performance Compiler V5.0
30 day trial version from Intel's Website
(Boy does this produce fast with vanilla C!!)

Rob James

email	robjames@europe.com
home	http://www.pocketmoon.com
April 2001

******************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <mathf.h>
#include <xmmintrin.h>  // intel simd header

//***  Note use of floats ***

//SIMD is basically 4 floats at once (or 2 doubles...)
typedef struct  {float x; float y; float z;} Vector;


// Standard Perlin 'optimised' noise functions
void	Noise2	(Vector  vec1,float *result);
void	Noise3	(Vector  vec1,float *result);


// Intel SIMD (PIII or CeleronII) optimised versions
void	Noise24	(Vector  vec1,Vector vec2,Vector vec3,Vector vec4,float *result);
void	Noise34	(Vector  vec1,Vector vec2,Vector vec3,Vector vec4,float *result);

//to initialise the 'noise' lookup arrays
void	PInitNoise();

//test functions for profiling
void	tnoise2();
void	tnoise3();


//from standard Noise functions
#define B	0x100
#define BM	0xff
//#define N	0x100000
#define N	 ((float)1048576.0)
#define NP	12   /* 2^N */
#define NM	0xfff
int		p[B + B + 2];

//alignment required for SIMD intrinsics to work
_declspec (align(32)) float g2[B + B + 2][2];//CHANGED FROM DOUBLE TO FLOAT and 3 to 2//
_declspec (align(32)) float g3[B + B + 2][3];//CHANGED FROM DOUBLE TO FLOAT and 3 to 2//


/* cubic spline interpolation */
#define s_curve(t)	( t * t * (3. - 2. * t) )

/* linear interpolation */
#define lerp(t, a, b)	( a + t * (b - a) )

// For noise3
#define at3(rx,ry,rz)	( rx * q[0] + ry * q[1] + rz * q[2] )

//for noise2
#define at2(rx,ry)	( rx * q[0] + ry * q[1] )

#define setup(u,b0,b1,r0,r1)\
	t = u + N;\
	b0 = ((int)t) & BM;\
	b1 = (b0+1) & BM;\
	r0 = t - (int)t;\
	r1 = r0 - 1.;



// anyone who's used Procedural noise should recognise these two functiomns ...
void Noise3(Vector vec, float*result)
{
	int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
	float rx0, rx1, ry0, ry1, rz0, rz1, *q, sx, sy, sz;
	float a, b, c, d, u, v,t;
	int i, j;

	setup(vec.x, bx0,bx1, rx0,rx1);
	setup(vec.y, by0,by1, ry0,ry1);
	setup(vec.z, bz0,bz1, rz0,rz1);

	i = p[ bx0 ];
	j = p[ bx1 ];

	b00 = p[ i + by0 ];
	b10 = p[ j + by0 ];
	b01 = p[ i + by1 ];
	b11 = p[ j + by1 ];

	sx = s_curve(rx0);
	sy = s_curve(ry0);
	sz = s_curve(rz0);

	q = g3[ b00 + bz0 ];	u = at3(rx0,ry0,rz0);
	q = g3[ b10 + bz0 ];	v = at3(rx1,ry0,rz0);
	a = lerp(sx, u, v);

	q = g3[ b01 + bz0 ];	u = at3(rx0,ry1,rz0);
	q = g3[ b11 + bz0 ];	v = at3(rx1,ry1,rz0);
	b = lerp(sx, u, v);

	c = lerp(sy, a, b);		/* interpolate in y at lo z */

	q = g3[ b00 + bz1 ];	u = at3(rx0,ry0,rz1);
	q = g3[ b10 + bz1 ];	v = at3(rx1,ry0,rz1);
	a = lerp(sx, u, v);

	q = g3[ b01 + bz1 ];	u = at3(rx0,ry1,rz1);
	q = g3[ b11 + bz1 ];	v = at3(rx1,ry1,rz1);
	b = lerp(sx, u, v);

	d = lerp(sy, a, b);		/* interpolate in y at hi z */
	*result =  1.5 * lerp(sz, c, d);	/* interpolate in z */

}


void Noise2(Vector vec, float *result)
{
	int	bx0, bx1, by0, by1, b00, b10, b01, b11;
	float	rx0, rx1, ry0, ry1, *q, sx, sy, a, b, u, v,t;
	int i, j;

	setup(vec.x, bx0,bx1, rx0,rx1);
	setup(vec.y, by0,by1, ry0,ry1);

	i = p[ bx0 ];
	j = p[ bx1 ];

	b00 = p[ i + by0 ];
	b10 = p[ j + by0 ];
	b01 = p[ i + by1 ];
	b11 = p[ j + by1 ];

	sx = s_curve(rx0);
	sy = s_curve(ry0);

	q = g2[ b00 ];		/* get random gradient */
	u = at2(rx0,ry0);		/* get weight on lo x side (lo y) */
	q = g2[ b10 ];
	v = at2(rx1,ry0);		/* get weight on hi x side (lo y) */
	a = lerp(sx, u, v);	/* get value at distance sx between u & v */

	/* similarly at hi y... */
	q = g2[ b01 ];	u = at2(rx0,ry1);
	q = g2[ b11 ];	v = at2(rx1,ry1);
	b = lerp(sx, u, v);

	*result = 1.5 * lerp(sy, a, b);	/* interpolate in y */	
}


// Now my 4-at-a-time functions using Intels SIMD intrinsics...
void Noise24(	Vector vec1,
				Vector vec2,
				Vector vec3,
				Vector vec4,
				float *result)/* MUST BE ALIGN(16)*/
{

	int	b000, b001, b002, b003,
		b100, b101, b102, b103,
		b010, b011, b012, b013,
		b110, b111, b112, b113;

	int	bx00, bx01, bx02, bx03,
		bx10, bx11, bx12, bx13,
		by00, by01, by02, by03,
		by10, by11, by12, by13;
	
	float	tx0,tx1,tx2,tx3,
			ty0,ty1,ty2,ty3;


	// first sight of an intrinsic!
	__m128	rx0, rx1, ry0, ry1,tmp,tmp2,  q0,q1;
	__m128  sx, sy, a, b, u, v , two, three;

	 int  i0, i1, i2, i3, j1,j2,j3,j0;

	// put (vec.x + 0x100000) into a temp var ??
	// ISSUE : bx0 etc are ints converted from floats
	// this is possible using __m64 _mm_cvtps_pi16( __m128 a) (composite)
    // but would entail using MMX (swich penalty) hmmm...

	//
	// When mixing MMX and SIMD we need to do a mm_empty - this could be slow
	// maybe not needed when using convert ???
	//

#define MM 1048576.0

	  int	itx1, itx2, itx3, itx0,
			 ity1, ity2, ity3, ity0;

/*
__m64 _mm_cvtps_pi16( __m128 a) (composite)
Convert the four single precision FP values in a to four signed 16-bit integer values.
r0 := (short)a0
r1 := (short)a1
r2 := (short)a2
r3 := (short)a3
*/

	tx0 = (vec1.x + MM);
	tx1 = (vec2.x + MM) ;
	tx2 = (vec3.x + MM) ;
	tx3 = (vec4.x + MM) ;

	itx0 = (int)(tx0); 
	itx1 = (int)(tx1);
	itx2 = (int)(tx2);
	itx3 = (int)(tx3);

	bx00 = itx0 & 0xff;
	bx10 = (bx00+1) & 0xff;

	bx01 = itx1 & 0xff;
	bx11 = (bx01+1) & 0xff;

	bx02 = itx2 & 0xff;
	bx12 = (bx02+1) & 0xff;
		
	bx03 = itx3 & 0xff;
	bx13 = (bx03+1) & 0xff;

	ty0 = (vec1.y + MM) ;
	ity0 = (int)(ty0);//-0.5
	by00 = ity0 & 0xff;
	by10 = (by00+1) & 0xff;
	
	ty1 = (vec2.y + MM) ;
	ity1 = (int)(ty1);//-.5
	by01 = ity1 & 0xff;
	by11 = (by01+1) & 0xff;

	ty2 = (vec3.y + MM) ;
	ity2 = (int)(ty2);//-0.5
	by02 = ity2 & 0xff;
	by12 = (by02+1) & 0xff;

	ty3 = (vec4.y + MM) ;	
	ity3 = (int)(ty3);//-0.5
	by03 = ity3 & 0xff;
	by13 = (by03+1) & 0xff;
	

	// could user  post incrementing pointers here...
	i0 = p[ bx00 ];
	j0 = p[ bx10 ];
	i1 = p[ bx01 ];
	j1 = p[ bx11 ];
	i2 = p[ bx02 ];
	j2 = p[ bx12 ];
	i3 = p[ bx03 ];
	j3 = p[ bx13 ];

	b000 = p[ i0 + by00 ];
	b010 = p[ i0 + by10 ];
	b100 = p[ j0 + by00 ];
	b110 = p[ j0 + by10 ];
	
	b001 = p[ i1 + by01 ];
	b011 = p[ i1 + by11 ];
	b101 = p[ j1 + by01 ];	
	b111 = p[ j1 + by11 ];
	
	b002 = p[ i2 + by02 ];
	b012 = p[ i2 + by12 ];
	b102 = p[ j2 + by02 ];	
	b112 = p[ j2 + by12 ];

	b003 = p[ i3 + by03 ];
	b013 = p[ i3 + by13 ];
	b103 = p[ j3 + by03 ];
	b113 = p[ j3 + by13 ];

	//!! try using FAST floattoint ... - tried it, slower :)
	rx0 =  _mm_set_ps(	tx0 - itx0,
						tx1 - itx1,
						tx2 - itx2,
						tx3 - itx3);


	// may be faster as another mm_set_ps with -1 at the end of each term
	tmp =  _mm_set_ps1(1.0 ) ;
	rx1 =  _mm_sub_ps(rx0, tmp);
	ry0 =  _mm_set_ps(	ty0 - ity0,
						ty1 - ity1,
						ty2 - ity2,
						ty3 - ity3);

	ry1 = _mm_sub_ps (ry0 ,tmp);


	//Now, these prefetches make a HUGE difference.
	//Without them the whole function is memory bound and no faster
	//than the original functions
	//
	//I have applied no real science in placement of these prefetches
	//I just chuck them in a few statements before I need them!
	_mm_prefetch((char*)&(g2[b000]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b001]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b002]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b003]),_MM_HINT_NTA);
	


	//you will now see original code as comments followed by simd equivalent...

	//	sx =  rx0 * rx0 * (3.0 - 2.0 * rx0);
	two  = _mm_set_ps1 (2.0)  ;	// Might be good to PRESTORE the constant arrays 1.0, 2.0 etc 
	sx   = _mm_mul_ps  (rx0,two);
	three  = _mm_set_ps1 (3.0);
	tmp2 = _mm_sub_ps  (three,sx);
	tmp  = _mm_mul_ps  (tmp2,rx0);
	sx   = _mm_mul_ps  (tmp, rx0);

	_mm_prefetch((char*)&(g2[b100]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b101]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b102]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b103]),_MM_HINT_NTA);

	
	//sy =  ry0 * ry0 * (3.0 - 2.0 * ry0);
	//tmp  = _mm_set_ps1 (2.0)  ;
	sy   = _mm_mul_ps  (ry0,two);

	//tmp  = _mm_set_ps1 (3.0);
	tmp2 = _mm_sub_ps  (three,sy);
	tmp  = _mm_mul_ps  (tmp2,ry0);
	sy   = _mm_mul_ps  (tmp, ry0);

	_mm_prefetch((char*)&(g2[b010]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b011]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b012]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b013]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b110]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b111]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b112]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g2[b113]),_MM_HINT_NTA);
	

//	q = g2[ b00 ];		/* get random gradient */      
//	u = rx0 * q[0] + ry0 * q[1];
	q0   = _mm_set_ps(g2[b000][0], g2[b001][0], g2[b002][0], g2[b003][0] );			
	q1   = _mm_set_ps(g2[b000][1], g2[b001][1], g2[b002][1], g2[b003][1] );			
	tmp  = _mm_mul_ps(q1, ry0);
	tmp2 = _mm_mul_ps(q0, rx0);
	u    = _mm_add_ps(tmp, tmp2);


//	q = g2[ b10 ];
//	v = rx1 * q[0] + ry0 * q[1] ;
	q0   = _mm_set_ps(g2[b100][0], g2[b101][0], g2[b102][0], g2[b103][0] );			
	q1   = _mm_set_ps(g2[b100][1], g2[b101][1], g2[b102][1], g2[b103][1] );			
	tmp  = _mm_mul_ps(q1, ry0);
	tmp2 = _mm_mul_ps(q0, rx1);
	v    = _mm_add_ps(tmp, tmp2);


//    a = u + sx * (v-u);  /* get value at distance sx between u & v */
	tmp  = _mm_sub_ps (v,u);
	tmp2 = _mm_mul_ps (tmp,sx);
	a    = _mm_add_ps (u,tmp2);

			
//	q = g2[ b01 ];
//	u = rx0*q[0] + ry1 * q[1];
	q0   = _mm_set_ps(g2[b010][0], g2[b011][0], g2[b012][0], g2[b013][0] );			
	q1   = _mm_set_ps(g2[b010][1], g2[b011][1], g2[b012][1], g2[b013][1] );			
	tmp  = _mm_mul_ps(q1, ry1);
	tmp2 = _mm_mul_ps(q0, rx0);
	u    = _mm_add_ps(tmp, tmp2);


//	q = g2[ b11 ];
//	v = rx1 * q[0] + ry1 * q[1];
	q0   = _mm_set_ps(g2[b110][0], g2[b111][0], g2[b112][0], g2[b113][0] );			
	q1   = _mm_set_ps(g2[b110][1], g2[b111][1], g2[b112][1], g2[b113][1] );			
	tmp  = _mm_mul_ps(q1, ry1);
	tmp2 = _mm_mul_ps(q0, rx1);
	v    = _mm_add_ps(tmp, tmp2);


//	b = u + sx * (v-u);
	tmp  = _mm_sub_ps (v,u);
	tmp2 = _mm_mul_ps (tmp,sx);
	b    = _mm_add_ps (u,tmp2);

//	result = 1.5 * (a + sy * (b-a))	/* interpolate in y */   why 1.5 ??
	tmp  = _mm_sub_ps(b,a);
	tmp2 = _mm_mul_ps(tmp,sy);
	tmp  = _mm_add_ps(tmp2,a);
	tmp2 = _mm_set_ps1(1.5);
	a    = _mm_mul_ps(tmp2,tmp);

	_mm_store_ps(result,  a );  // aligned


	//phew!!
	// As you see doing things 4-at-a-time does take some getting used to!
} 


void Noise34(Vector vec1,
		Vector vec2,
		Vector vec3,
		Vector vec4,
		float *result)/* MUST BE ALIGN(16)*/
{

	// same as Noise24 but with the extra Z component included...

	int	b000, b001, b002, b003,
		b100, b101, b102, b103,
		b010, b011, b012, b013,
		b110, b111, b112, b113;

	int	bx00, bx01, bx02, bx03,
		bx10, bx11, bx12, bx13,
		by00, by01, by02, by03,
		by10, by11, by12, by13,
		bz00, bz01, bz02, bz03,
		bz10, bz11, bz12, bz13;
	
	float	tx0,tx1,tx2,tx3,
		ty0,ty1,ty2,ty3,
		tz0,tz1,tz2,tz3;

	__m128	rx0, rx1, ry0, ry1, rz0, rz1, tmp,tmp2,tmp3,tmp4, q0,q1,q2;
	__m128	sx, sy, sz, a, b,c ,d,  u, v , two, three;

	 int  i0, i1, i2, i3, j1,j2,j3,j0;

	// put (vec.x + 0x100000) into a temp var ??
	// ISSUE : bx0 etc are ints converted from floats
	// this is possible using __m64 _mm_cvtps_pi16( __m128 a) (composite)
        // but would entail using MMX hmmm

	//
	// When mixing MMX and SIMD we need to do a mm_empty - this could be slow
	// maybe not needed when using convert ???
	//

#define MM 1048576.0

	int	itx1, itx2, itx3, itx0,
		ity1, ity2, ity3, ity0,
		itz1, itz2, itz3, itz0;

/*
__m64 _mm_cvtps_pi16( __m128 a) (composite)
Convert the four single precision FP values in a to four signed 16-bit integer values.
r0 := (short)a0
r1 := (short)a1
r2 := (short)a2
r3 := (short)a3
*/

	tx0 = (vec1.x + MM);
	tx1 = (vec2.x + MM) ;
	tx2 = (vec3.x + MM) ;
	tx3 = (vec4.x + MM) ;

	itx0 = (int)(tx0);    
	itx1 = (int)(tx1);
	itx2 = (int)(tx2);
	itx3 = (int)(tx3);

	bx00 = itx0 & 0xff;
	bx10 = (bx00+1) & 0xff;

	bx01 = itx1 & 0xff;
	bx11 = (bx01+1) & 0xff;

	bx02 = itx2 & 0xff;
	bx12 = (bx02+1) & 0xff;
		
	bx03 = itx3 & 0xff;
	bx13 = (bx03+1) & 0xff;



	ty0 = (vec1.y + MM) ;
	ity0 = (int)(ty0);
	by00 = ity0 & 0xff;
	by10 = (by00+1) & 0xff;
	
	ty1 = (vec2.y + MM) ;
	ity1 = (int)(ty1);
	by01 = ity1 & 0xff;
	by11 = (by01+1) & 0xff;

	ty2 = (vec3.y + MM) ;
	ity2 = (int)(ty2);
	by02 = ity2 & 0xff;
	by12 = (by02+1) & 0xff;

	ty3 = (vec4.y + MM) ;	
	ity3 = (int)(ty3);
	by03 = ity3 & 0xff;
	by13 = (by03+1) & 0xff;


	// ZED
	tz0 = (vec1.z + MM) ;
	itz0 = (int)(tz0);
	bz00 = itz0 & 0xff;
	bz10 = (bz00+1) & 0xff;
	
	tz1 = (vec2.z + MM) ;
	itz1 = (int)(tz1);
	bz01 = itz1 & 0xff;
	bz11 = (bz01+1) & 0xff;

	tz2 = (vec3.z + MM) ;
	itz2 = (int)(tz2);
	bz02 = itz2 & 0xff;
	bz12 = (bz02+1) & 0xff;

	tz3 = (vec4.z + MM) ;	
	itz3 = (int)(tz3);
	bz03 = itz3 & 0xff;
	bz13 = (bz03+1) & 0xff;

	

	// could user  post incrementing pointers here...
	i0 = p[ bx00 ];
	j0 = p[ bx10 ];
	i1 = p[ bx01 ];
	j1 = p[ bx11 ];
	i2 = p[ bx02 ];
	j2 = p[ bx12 ];
	i3 = p[ bx03 ];
	j3 = p[ bx13 ];

	b000 = p[ i0 + by00 ];
	b010 = p[ i0 + by10 ];
	b100 = p[ j0 + by00 ];
	b110 = p[ j0 + by10 ];
	
	b001 = p[ i1 + by01 ];
	b011 = p[ i1 + by11 ];
	b101 = p[ j1 + by01 ];	
	b111 = p[ j1 + by11 ];
	
	b002 = p[ i2 + by02 ];
	b012 = p[ i2 + by12 ];
	b102 = p[ j2 + by02 ];	
	b112 = p[ j2 + by12 ];

	b003 = p[ i3 + by03 ];
	b013 = p[ i3 + by13 ];
	b103 = p[ j3 + by03 ];
	b113 = p[ j3 + by13 ];


	// may be faster as another mm_set_ps with -1 at the end of each term
	tmp =  _mm_set_ps1(1.0) ;

	//!! try using FAST floattoint ...???
	rx0 =  _mm_set_ps(tx0 - itx0,
				tx1 - itx1,
				tx2 - itx2,
				tx3 - itx3);
	rx1 =  _mm_sub_ps(rx0, tmp);

	ry0 =  _mm_set_ps(ty0 - ity0,
				ty1 - ity1,
				ty2 - ity2,
				ty3 - ity3);
	ry1 = _mm_sub_ps (ry0 ,tmp);


	rz0 =  _mm_set_ps(	tz0 - itz0,
						tz1 - itz1,
						tz2 - itz2,
						tz3 - itz3);
	rz1 = _mm_sub_ps (rz0 ,tmp);


	_mm_prefetch((char*)&(g3[b000+bz00]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b001+bz01]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b002+bz02]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b003+bz03]),_MM_HINT_NTA);
	

//	sx =  rx0 * rx0 * (3.0 - 2.0 * rx0);
	two  = _mm_set_ps1 (2.0)  ;	// Might be good to PRESTORE the constant arrays 1.0, 2.0 etc 
	sx   = _mm_mul_ps  (rx0,two);
	three  = _mm_set_ps1 (3.0);
	tmp2 = _mm_sub_ps  (three,sx);
	tmp  = _mm_mul_ps  (tmp2,rx0);
	sx   = _mm_mul_ps  (tmp, rx0);

	_mm_prefetch((char*)&(g3[b100+bz00]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b101+bz01]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b102+bz02]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b103+bz03]),_MM_HINT_NTA);

	
	//sy =  ry0 * ry0 * (3.0 - 2.0 * ry0);
	//tmp  = _mm_set_ps1 (2.0)  ;
	sy   = _mm_mul_ps  (ry0,two);

	//tmp  = _mm_set_ps1 (3.0);
	tmp2 = _mm_sub_ps  (three,sy);
	tmp  = _mm_mul_ps  (tmp2,ry0);
	sy   = _mm_mul_ps  (tmp, ry0);

	_mm_prefetch((char*)&(g3[b010+bz00]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b011+bz01]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b012+bz02]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b013+bz03]),_MM_HINT_NTA);

	//sz =  rz0 * rz0 * (3.0 - 2.0 * rz0);
	//tmp  = _mm_set_ps1 (2.0)  ;
	sz   = _mm_mul_ps  (rz0,two);
	//tmp  = _mm_set_ps1 (3.0);
	tmp2 = _mm_sub_ps  (three,sz);
	tmp  = _mm_mul_ps  (tmp2,rz0);
	sz   = _mm_mul_ps  (tmp, rz0);

	_mm_prefetch((char*)&(g3[b110+bz00]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b111+bz01]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b112+bz02]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b113+bz03]),_MM_HINT_NTA);
	

//	q = g3[ b00 ];		/* get random gradient */      
//	u = rx0 * q[0] + ry0 * q[1] + rz0 * q[2];
	q0   = _mm_set_ps(g3[b000+bz00][0], g3[b001+bz01][0], g3[b002+bz02][0], g3[b003+bz03][0] );
	q1   = _mm_set_ps(g3[b000+bz00][1], g3[b001+bz01][1], g3[b002+bz02][1], g3[b003+bz03][1] );
	q2   = _mm_set_ps(g3[b000+bz00][2], g3[b001+bz01][2], g3[b002+bz02][2], g3[b003+bz03][2] );

	tmp3 = _mm_mul_ps(q2, rz0);
	tmp  = _mm_mul_ps(q1, ry0);
	tmp2 = _mm_mul_ps(q0, rx0);
	tmp4 = _mm_add_ps(tmp, tmp2);
	u    = _mm_add_ps(tmp4, tmp3);


	
	_mm_prefetch((char*)&(g3[b000+bz10]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b001+bz11]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b002+bz12]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b003+bz13]),_MM_HINT_NTA);
//	q = g3[ b10 + bz0 ];
//	v = rx1 * q[0] + ry0 * q[1] ;
	q0   = _mm_set_ps(g3[b100+bz00][0], g3[b101+bz01][0], g3[b102+bz02][0], g3[b103+bz03][0] );			
	q1   = _mm_set_ps(g3[b100+bz00][1], g3[b101+bz01][1], g3[b102+bz02][1], g3[b103+bz03][1] );			
	q2   = _mm_set_ps(g3[b100+bz00][2], g3[b101+bz01][2], g3[b102+bz02][2], g3[b103+bz03][2] );			

	tmp3 = _mm_mul_ps(q2, rz0);
	tmp  = _mm_mul_ps(q1, ry0);
	tmp2 = _mm_mul_ps(q0, rx1);
	tmp4 = _mm_add_ps(tmp, tmp2);
	v    = _mm_add_ps(tmp4, tmp3);



//    a = u + sx * (v-u);  /* get value at distance sx between u & v */
	tmp  = _mm_sub_ps (v,u);
	tmp2 = _mm_mul_ps (tmp,sx);
	a    = _mm_add_ps (u,tmp2);

			
	
	_mm_prefetch((char*)&(g3[b100+bz10]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b101+bz11]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b102+bz12]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b103+bz13]),_MM_HINT_NTA);
//	q = g3[ b01 ];
//	u = rx0*q[0] + ry1 * q[1];
	q0   = _mm_set_ps(g3[b010+bz00][0], g3[b011+bz01][0], g3[b012+bz02][0], g3[b013+bz03][0] );
	q1   = _mm_set_ps(g3[b010+bz00][1], g3[b011+bz01][1], g3[b012+bz02][1], g3[b013+bz03][1] );
	q2   = _mm_set_ps(g3[b010+bz00][2], g3[b011+bz01][2], g3[b012+bz02][2], g3[b013+bz03][2] );

	tmp3 = _mm_mul_ps(q2, rz0);
	tmp  = _mm_mul_ps(q1, ry1);
	tmp2 = _mm_mul_ps(q0, rx0);
	tmp4 = _mm_add_ps(tmp, tmp2);
	u    = _mm_add_ps(tmp4, tmp3);


	
	_mm_prefetch((char*)&(g3[b010+bz10]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b011+bz11]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b012+bz12]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b013+bz13]),_MM_HINT_NTA);
//	q = g3[ b11 ];
//	v = rx1 * q[0] + ry1 * q[1];
	q0   = _mm_set_ps(g3[b110+bz00][0], g3[b111+bz01][0], g3[b112+bz02][0], g3[b113+bz03][0] );
	q1   = _mm_set_ps(g3[b110+bz00][1], g3[b111+bz01][1], g3[b112+bz02][1], g3[b113+bz03][1] );
	q2   = _mm_set_ps(g3[b110+bz00][2], g3[b111+bz01][2], g3[b112+bz02][2], g3[b113+bz03][2] );

	tmp3 = _mm_mul_ps(q2, rz0);
	tmp  = _mm_mul_ps(q1, ry1);
	tmp2 = _mm_mul_ps(q0, rx1);
	tmp4 = _mm_add_ps(tmp, tmp2);
	v    = _mm_add_ps(tmp4, tmp3);



//	b = u + sx * (v-u);
	tmp  = _mm_sub_ps (v,u);
	tmp2 = _mm_mul_ps (tmp,sx);
	b    = _mm_add_ps (u,tmp2);


//	c = a + sy * (a-b);
	tmp  = _mm_sub_ps (b,a);
	tmp2 = _mm_mul_ps (tmp,sy);
	c    = _mm_add_ps (a,tmp2);

/*NEW */

	_mm_prefetch((char*)&(g3[b110+bz10]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b111+bz11]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b112+bz12]),_MM_HINT_NTA);
	_mm_prefetch((char*)&(g3[b113+bz13]),_MM_HINT_NTA);

//	q = g3[ b00 + bz1 ];
//	u = rx0*q[0] + ry1 * q[1];
	q0   = _mm_set_ps(g3[b000+bz10][0], g3[b001+bz11][0], g3[b002+bz12][0], g3[b003+bz13][0] );
	q1   = _mm_set_ps(g3[b000+bz10][1], g3[b001+bz11][1], g3[b002+bz12][1], g3[b003+bz13][1] );
	q2   = _mm_set_ps(g3[b000+bz10][2], g3[b001+bz11][2], g3[b002+bz12][2], g3[b003+bz13][2] );

	tmp3 = _mm_mul_ps(q2, rz1);
	tmp  = _mm_mul_ps(q1, ry0);
	tmp2 = _mm_mul_ps(q0, rx0);
	tmp4 = _mm_add_ps(tmp, tmp2);
	u    = _mm_add_ps(tmp4, tmp3);


//	q = g3[ b10 + bz1 ];
//	v = rx1 * q[0] + ry1 * q[1];
	q0   = _mm_set_ps(g3[b100+bz10][0], g3[b101+bz11][0], g3[b102+bz12][0], g3[b103+bz13][0] );
	q1   = _mm_set_ps(g3[b100+bz10][1], g3[b101+bz11][1], g3[b102+bz12][1], g3[b103+bz13][1] );
	q2   = _mm_set_ps(g3[b100+bz10][2], g3[b101+bz11][2], g3[b102+bz12][2], g3[b103+bz13][2] );

	tmp3 = _mm_mul_ps(q2, rz1);
	tmp  = _mm_mul_ps(q1, ry0);
	tmp2 = _mm_mul_ps(q0, rx1);
	tmp4 = _mm_add_ps(tmp, tmp2);
	v    = _mm_add_ps(tmp4, tmp3);

//	a = u + sx * (v-u);
	tmp  = _mm_sub_ps (v,u);
	tmp2 = _mm_mul_ps (tmp,sx);
	a    = _mm_add_ps (u,tmp2);




//	q = g3[ b01 + bz1 ];
//	v = rx1 * q[0] + ry1 * q[1];
	q0   = _mm_set_ps(g3[b010+bz10][0], g3[b011+bz11][0], g3[b012+bz12][0], g3[b013+bz13][0] );
	q1   = _mm_set_ps(g3[b010+bz10][1], g3[b011+bz11][1], g3[b012+bz12][1], g3[b013+bz13][1] );
	q2   = _mm_set_ps(g3[b010+bz10][2], g3[b011+bz11][2], g3[b012+bz12][2], g3[b013+bz13][2] );

	tmp3 = _mm_mul_ps(q2, rz1);
	tmp  = _mm_mul_ps(q1, ry1);
	tmp2 = _mm_mul_ps(q0, rx0);
	tmp4 = _mm_add_ps(tmp, tmp2);
	u    = _mm_add_ps(tmp4, tmp3);

	q0   = _mm_set_ps(g3[b110+bz10][0], g3[b111+bz11][0], g3[b112+bz12][0], g3[b113+bz13][0] );
	q1   = _mm_set_ps(g3[b110+bz10][1], g3[b111+bz11][1], g3[b112+bz12][1], g3[b113+bz13][1] );
	q2   = _mm_set_ps(g3[b110+bz10][2], g3[b111+bz11][2], g3[b112+bz12][2], g3[b113+bz13][2] );

	tmp3 = _mm_mul_ps(q2, rz1);
	tmp  = _mm_mul_ps(q1, ry1);
	tmp2 = _mm_mul_ps(q0, rx1);
	tmp4 = _mm_add_ps(tmp, tmp2);
	v    = _mm_add_ps(tmp4, tmp3);

	tmp  = _mm_sub_ps (v,u);
	tmp2 = _mm_mul_ps (tmp,sx);
	b    = _mm_add_ps (u,tmp2);

	tmp  = _mm_sub_ps (b,a);
	tmp2 = _mm_mul_ps (tmp,sy);
	d    = _mm_add_ps (a,tmp2);


//	result = 1.5 * (c + sz * (d-c))	/* interpolate in y */   why 1.5 ??
	tmp  = _mm_sub_ps(d,c);
	tmp2 = _mm_mul_ps(tmp,sz);
	tmp  = _mm_add_ps(tmp2,c);
	tmp2 = _mm_set_ps1(1.5);
	a    = _mm_mul_ps(tmp2,tmp);

	_mm_store_ps(result,  a );  // aligned
} 






#include "noisetab.c"
void PInitNoise()
{
int i;

//set up our two random number lookup tables
//one for 2d noise the other for 3d...

	for (i= 0; i < B+B+2; ++i)
	{
		p[i]= p_precomputed[i];
		g3[i][0]	= g2[i][0]	= g_precomputed[i][0];
		g3[i][1] 	= g2[i][1]	= g_precomputed[i][1];              
		g3[i][2]			= g_precomputed[i][2];              
	}
}


void tnoise2 ()
{

int x;
Vector v;

// some aligned storage for test reults
_declspec (align(16)) float a[1000] ;
_declspec (align(16)) float b[1000]  ;


	v.x = (float)rand()/10000.0;
	v.y = (float)rand()/10000.0;
	v.z = (float)rand()/10000.0;

	//printf("%f %f %f   :",v.x,v.y,v.z);

	for (x=0;x<1000;x++)
		Noise2(v, &(a[x]));

	for (x=0;x<1000;x+=4)
		Noise24(v,v,v,v, &(b[x]));

	//printf("%f   %f   %f \n",a[0]-b[0],a[456]-b[456],a[999]-b[999]);

}



void tnoise3 ()
{

_declspec (align(16)) float a[1000]  ;
_declspec (align(16)) float b[1000]  ;

int x;
Vector v;

	v.x = (float)rand()/10000.0;
	v.y = (float)rand()/10000.0;
	v.z = (float)rand()/10000.0;

//	printf("%f %f %f\t:",v.x,v.y,v.z);

	for (x=0;x<1000;x++)
		Noise3(v, &(a[x]));

	for (x=0;x<1000;x+=4)
		Noise34(v, v, v, v, &(b[x]));

//	printf("%f %f   %f %f\n",a[0],b[0],a[999],b[999]);
}


//Almost there...
int main(int argc, char** argv)
{
int x;
	

	PInitNoise();
	for (x=0;x < 1000;x++)
	{
		tnoise2();
		tnoise3();
	}
}

