/* NOTE simon (28/04/25 17:37:11): cl cubic.c -Fecubic.exe -nologo -GR- -EHa- -Oi -W4 -wd4100 -wd4189 -wd4201 -FC -Zi -Zo -diagnostics:caret -diagnostics:column -O2 -MT -GS- -WX -std:c11 -link -INCREMENTAL:NO -opt:ref,icf */ #include #include #include #include #define _USE_MATH_DEFINES #include typedef float f32; typedef double f64; typedef uintptr_t umm; typedef uint8_t u8; typedef uint32_t u32; typedef uint64_t u64; typedef int32_t s32; typedef uint32_t b32; #define cast( type, expr ) ( type ) ( expr ) #define array_count( array ) ( ( umm ) ( sizeof( array ) / sizeof( ( array )[ 0 ] ) ) ) #define _assert( expr ) if ( !( expr ) ) { __debugbreak( ); } #define random_32_multiplier 6364136223846793005ull typedef struct random_32_t { u64 state; u64 increment; } random_32_t; random_32_t random_32_make( u64 seed, u64 sequence ) { random_32_t result = { 0 }; result.increment = ( sequence << 1 ) | 1; result.state = result.increment + seed; result.state = ( result.state * random_32_multiplier ) + result.increment; return result; } u32 random_32_get( random_32_t* random_32 ) { /* NOTE simon: 64 xsh rr 32 */ u64 temp = random_32->state; random_32->state = ( random_32->state * random_32_multiplier ) + random_32->increment; /* NOTE simon: 2^5 => 32 bit rotation. */ u8 rotation = ( u8 ) ( temp >> ( 64 - 5 ) ); temp = ( ( temp >> 18 ) ^ temp ) >> 27; u32 result = ( u32 ) ( temp & 0xffffffff ); result = ( result >> rotation ) | ( result << ( 32 - rotation ) ); return result; } f32 random_32_get_f32( random_32_t* random_32 ) { u32 n = random_32_get( random_32 ); f32 f = cast( f32, n >> 8 ) * 0x1.0p-24f; return f; } static f32 interpolation_linear( f32 start, f32 end, f32 t ) { f32 result = start * ( 1 - t ) + end * t; return result; } static f32 interpolation_cubic( f32 start, f32 cp_1, f32 cp_2, f32 end, f32 t ) { f32 u = 1 - t; f32 result = ( u * u * u * start ) + ( 3 * u * u * t * cp_1 ) + ( 3 * u * t * t * cp_2 ) + ( t * t * t * end ); return result; } #include #include #include _Alignas( 64 ) u8 shuffle_lut[ 16 * 16 ] = { 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b0000 0, 1, 2, 3, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b0001 4, 5, 6, 7, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b0010 0, 1, 2, 3, 4, 5, 6, 7, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b0011 8, 9, 10, 11, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b0100 0, 1, 2, 3, 8, 9, 10, 11, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b0101 4, 5, 6, 7, 8, 9, 10, 11, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b0110 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0x80, 0x80, 0x80, 0x80, // 0b0111 12, 13, 14, 15, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b1000 0, 1, 2, 3, 12, 13, 14, 15, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b1001 4, 5, 6, 7, 12, 13, 14, 15, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b1010 0, 1, 2, 3, 4, 5, 6, 7, 12, 13, 14, 15, 0x80, 0x80, 0x80, 0x80, // 0b1011 8, 9, 10, 11, 12, 13, 14, 15, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, 0x80, // 0b1100 0, 1, 2, 3, 8, 9, 10, 11, 12, 13, 14, 15, 0x80, 0x80, 0x80, 0x80, // 0b1101 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 0x80, 0x80, 0x80, 0x80, // 0b1110 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 // 0b1111 }; _Alignas( 64 ) u8 shuffle_lut_2[ 16 * 16 ] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // 0b0000 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b0001 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, // 0b0010 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, // 0b0011 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, // 0b0100 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, // 0b0101 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, // 0b0110 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, // 0b0111 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1000 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1001 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1010 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1011 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1100 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1101 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1110 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1111 }; _Alignas( 64 ) u8 shuffle_lut_3[ 16 * 16 ] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // 0b0000 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b0001 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, // 0b0010 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b0011 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, // 0b0100 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b0101 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, // 0b0110 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b0111 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, 12, 13, 14, 15, // 0b1000 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b1001 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, // 0b1010 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b1011 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, 8, 9, 10, 11, // 0b1100 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b1101 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7, // 0b1110 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, // 0b1111 }; typedef struct cubic_state_t { f32 a, b, c; f32 start; #if 1 f32 cp_1, cp_2, end; #endif } cubic_state_t; static void cubic_state_initialize( cubic_state_t* state, f32 start, f32 cp_1, f32 cp_2, f32 end ) { *state = ( cubic_state_t ) { 0 }; state->a = ( -start + 3 * cp_1 - 3 * cp_2 + end ); state->b = ( 3 * start - 6 * cp_1 + 3 * cp_2 ); state->c = ( -3 * start + 3 * cp_1 ); state->start = start; #if 1 state->cp_1 = cp_1; state->cp_2 = cp_2; state->end = end; #endif } static f32 interpolation_cubic_t_from_x_bisection_simd_1( cubic_state_t* state, f32 x, f32 epsilon ) { f32 d = state->start; f32 result = 0; /* NOTE simon (24/04/25 18:41:29): The order is inverted to my mental model. */ __m128 factors_1 = _mm_set_ps( 0.2f, 0.4f, 0.6f, 0.8f ); __m128 factors_2 = _mm_set_ps( 0.8f, 0.6f, 0.4f, 0.2f ); __m128 ts = factors_2; __m128 as = _mm_set1_ps( state->a ); __m128 bs = _mm_set1_ps( state->b ); __m128 cs = _mm_set1_ps( state->c ); __m128 ds = _mm_set1_ps( d ); __m128 epsilons = _mm_set_ps1( epsilon ); __m128 xs = _mm_set1_ps( x ); __m128 ls = _mm_set1_ps( 0 ); __m128 rs = _mm_set1_ps( 1 ); while ( 1 ) { /* NOTE simon (26/04/25 13:57:04): Compute cubic. */ __m128 tt = _mm_mul_ps( ts, ts ); __m128 ttt = _mm_mul_ps( tt, ts ); __m128 attt = _mm_mul_ps( as, ttt ); __m128 btt = _mm_mul_ps( bs, tt ); __m128 ct = _mm_mul_ps( cs, ts ); __m128 f1 = _mm_add_ps( attt, btt ); __m128 f2 = _mm_add_ps( ct, ds ); __m128 f = _mm_add_ps( f1, f2 ); __m128 f_minus_x = _mm_sub_ps( f, xs ); /* NOTE simon (26/04/25 13:57:12): Absolute value */ __m128 n_zero = _mm_set1_ps( -0.0f ); __m128 abs = _mm_andnot_ps( n_zero, f_minus_x ); __m128 cmp = _mm_cmple_ps( abs, epsilons ); /* NOTE simon (26/04/25 13:57:37): Have we found the result ? */ u32 mask = _mm_movemask_ps( cmp ); if ( mask ) { /* NOTE simon (26/04/25 13:58:08): Extract the result. */ /* NOTE simon (06/05/25 16:33:16): We can use the same look up table for extraction and packing below, so we do to save 256 bytes in the cache. */ __m128i shuffle_mask = *( __m128i* ) ( shuffle_lut_2 + ( mask * 16 ) ); __m128i shuffled = _mm_shuffle_epi8( _mm_castps_si128( ts ), shuffle_mask ); result = _mm_cvtss_f32( _mm_castsi128_ps( shuffled ) ); break; } /* NOTE simon (26/04/25 13:58:20): Update left and right limits. */ __m128 setl = _mm_cmplt_ps( f, xs ); u32 lmask = _mm_movemask_ps( setl ); u32 rmask = 15 ^ lmask; { __m128i shuffle_mask = *( __m128i* ) ( shuffle_lut_2 + ( lmask * 16 ) ); /* NOTE simon (29/04/25 16:20:07): Storing in a variable isn't great. The compiler can still do the right thing and use the shuffle with the memory operand, but it's not guaranteed. */ __m128i shuffled = _mm_shuffle_epi8( _mm_castps_si128( ts ), shuffle_mask ); __m128i l_blend_mask = _mm_set1_epi8( lmask ? 0xff : 0 ); ls = _mm_blendv_ps( ls, _mm_castsi128_ps( shuffled ), _mm_castsi128_ps( l_blend_mask ) ); } { __m128i shuffle_mask = *( __m128i* ) ( shuffle_lut_3 + ( rmask * 16 ) ); __m128i shuffled = _mm_shuffle_epi8( _mm_castps_si128( ts ), shuffle_mask ); __m128i r_blend_mask = _mm_set1_epi8( rmask ? 0xff : 0 ); rs = _mm_blendv_ps( rs, _mm_castsi128_ps( shuffled ), _mm_castsi128_ps( r_blend_mask ) ); } /* NOTE simon (24/04/25 18:12:28): Lerp */ __m128 left = _mm_mul_ps( ls, factors_1 ); __m128 right = _mm_mul_ps( rs, factors_2 ); ts = _mm_add_ps( left, right ); } return result; } static f32 interpolation_cubic_t_from_x_bisection_simd_2( cubic_state_t* state, f32 x, f32 epsilon ) { f32 d = state->start; f32 result = 0; /* NOTE simon (24/04/25 18:41:29): The order is inverted to my mental model. */ __m128 factors_1 = _mm_set_ps( 0.2f, 0.4f, 0.6f, 0.8f ); __m128 factors_2 = _mm_set_ps( 0.8f, 0.6f, 0.4f, 0.2f ); __m128 ts = factors_2; __m128 as = _mm_set1_ps( state->a ); __m128 bs = _mm_set1_ps( state->b ); __m128 cs = _mm_set1_ps( state->c ); __m128 ds = _mm_set1_ps( d ); __m128 epsilons = _mm_set_ps1( epsilon ); __m128 xs = _mm_set1_ps( x ); __m128 ls = _mm_set1_ps( 0 ); __m128 rs = _mm_set1_ps( 1 ); while ( 1 ) { /* NOTE simon (26/04/25 13:57:04): Compute cubic. */ __m128 tt = _mm_mul_ps( ts, ts ); __m128 ttt = _mm_mul_ps( tt, ts ); __m128 attt = _mm_mul_ps( as, ttt ); __m128 btt = _mm_mul_ps( bs, tt ); __m128 ct = _mm_mul_ps( cs, ts ); __m128 f1 = _mm_add_ps( attt, btt ); __m128 f2 = _mm_add_ps( ct, ds ); __m128 f = _mm_add_ps( f1, f2 ); __m128 f_minus_x = _mm_sub_ps( f, xs ); /* NOTE simon (26/04/25 13:57:12): Absolute value */ __m128 n_zero = _mm_set1_ps( -0.0f ); __m128 abs = _mm_andnot_ps( n_zero, f_minus_x ); /* NOTE simon (26/04/25 13:57:37): Have we found the result ? */ __m128 cmp = _mm_cmple_ps( abs, epsilons ); u32 mask = _mm_movemask_ps( cmp ); if ( mask ) { /* NOTE simon (26/04/25 13:58:08): Extract the result. */ __m128 shuffled; if ( mask & 0x1 ) { shuffled = _mm_shuffle_ps( ts, ts, 0 ); } else if ( mask & 0x2 ) { shuffled = _mm_shuffle_ps( ts, ts, 1 ); } else if ( mask & 0x4 ) { shuffled = _mm_shuffle_ps( ts, ts, 2 ); } else { // _assert( mask & 0x8 ); shuffled = _mm_shuffle_ps( ts, ts, 3 ); } result = _mm_cvtss_f32( shuffled ); break; } /* NOTE simon (26/04/25 13:58:20): Update left and right limits. */ __m128 setl = _mm_cmplt_ps( f, xs ); u32 lmask = _mm_movemask_ps( setl ); u32 rmask = ~lmask; if ( lmask & 0x8 ) { // 0b11111111 ls = _mm_shuffle_ps( ts, ts, 0xff ); } else if ( lmask & 0x4 ) { // 0b10101010 ls = _mm_shuffle_ps( ts, ts, 0xaa ); } else if ( lmask & 0x2 ) { // 0b01010101 ls = _mm_shuffle_ps( ts, ts, 0x55 ); } else if ( lmask & 0x1 ) { // 0b00000000 ls = _mm_shuffle_ps( ts, ts, 0 ); } if ( rmask & 0x1 ) { rs = _mm_shuffle_ps( ts, ts, 0 ); } else if ( rmask & 0x2 ) { rs = _mm_shuffle_ps( ts, ts, 0x55 ); } else if ( rmask & 0x4 ) { rs = _mm_shuffle_ps( ts, ts, 0xaa ); } else if ( rmask & 0x8 ) { rs = _mm_shuffle_ps( ts, ts, 0xff ); } /* NOTE simon (24/04/25 18:12:28): Lerp */ __m128 left = _mm_mul_ps( ls, factors_1 ); __m128 right = _mm_mul_ps( rs, factors_2 ); ts = _mm_add_ps( left, right ); } return result; } static f32 interpolation_cubic_t_from_x_bisection_simd_3( cubic_state_t* state, f32 x, f32 epsilon ) { f32 d = state->start; f32 result = 0; /* NOTE simon (24/04/25 18:41:29): The order is inverted to my mental model. */ __m128 factors_1 = _mm_set_ps( 0.2f, 0.4f, 0.6f, 0.8f ); __m128 factors_2 = _mm_set_ps( 0.8f, 0.6f, 0.4f, 0.2f ); __m128 ts = factors_2; __m128 as = _mm_set1_ps( state->a ); __m128 bs = _mm_set1_ps( state->b ); __m128 cs = _mm_set1_ps( state->c ); __m128 ds = _mm_set1_ps( d ); __m128 epsilons = _mm_set_ps1( epsilon ); __m128 xs = _mm_set1_ps( x ); __m128 ls = _mm_set1_ps( 0 ); __m128 rs = _mm_set1_ps( 1 ); while ( 1 ) { /* NOTE simon (26/04/25 13:57:04): Compute cubic. */ __m128 tt = _mm_mul_ps( ts, ts ); __m128 ttt = _mm_mul_ps( tt, ts ); __m128 attt = _mm_mul_ps( as, ttt ); __m128 btt = _mm_mul_ps( bs, tt ); __m128 ct = _mm_mul_ps( cs, ts ); __m128 f1 = _mm_add_ps( attt, btt ); __m128 f2 = _mm_add_ps( ct, ds ); __m128 f = _mm_add_ps( f1, f2 ); __m128 f_minus_x = _mm_sub_ps( f, xs ); /* NOTE simon (26/04/25 13:57:12): Absolute value */ __m128 n_zero = _mm_set1_ps( -0.0f ); __m128 abs = _mm_andnot_ps( n_zero, f_minus_x ); __m128 cmp = _mm_cmple_ps( abs, epsilons ); /* NOTE simon (26/04/25 13:57:37): Have we found the result ? */ u32 mask = _mm_movemask_ps( cmp ); if ( mask ) { /* NOTE simon (26/04/25 13:58:08): Extract the result. */ __m128 shuffled; if ( mask & 0x1 ) { shuffled = _mm_shuffle_ps( ts, ts, 0 ); } else if ( mask & 0x2 ) { shuffled = _mm_shuffle_ps( ts, ts, 1 ); } else if ( mask & 0x4 ) { shuffled = _mm_shuffle_ps( ts, ts, 2 ); } else { // _assert( mask & 0x8 ); shuffled = _mm_shuffle_ps( ts, ts, 3 ); } result = _mm_cvtss_f32( shuffled ); break; } /* NOTE simon (26/04/25 13:58:20): Update left and right limits. */ __m128 setl = _mm_cmplt_ps( f, xs ); u32 lmask = _mm_movemask_ps( setl ); u32 rmask = ~lmask; if ( lmask & 0x8 ) { // 0b11111111 ls = _mm_shuffle_ps( ts, ts, 0xff ); } else if ( lmask & 0x4 ) { // 0b10101010 ls = _mm_shuffle_ps( ts, ts, 0xaa ); } else if ( lmask & 0x2 ) { // 0b01010101 ls = _mm_shuffle_ps( ts, ts, 0x55 ); } else if ( lmask & 0x1 ) { // 0b00000000 ls = _mm_shuffle_ps( ts, ts, 0 ); } if ( rmask & 0x1 ) { rs = _mm_shuffle_ps( ts, ts, 0 ); } else if ( rmask & 0x2 ) { rs = _mm_shuffle_ps( ts, ts, 0x55 ); } else if ( rmask & 0x4 ) { rs = _mm_shuffle_ps( ts, ts, 0xaa ); } else if ( rmask & 0x8 ) { rs = _mm_shuffle_ps( ts, ts, 0xff ); } /* NOTE simon (24/04/25 18:12:28): Lerp */ __m128 left = _mm_mul_ps( ls, factors_1 ); __m128 right = _mm_mul_ps( rs, factors_2 ); ts = _mm_add_ps( left, right ); } return result; } static f32 interpolation_cubic_t_from_x_bisection( cubic_state_t* state, f32 x, f32 epsilon ) { f32 d = state->start; f32 l = 0; f32 w = 1; f32 t = 0; while ( 1 ) { w *= 0.5f; t = l + w; f32 f = state->a * t * t * t + state->b * t * t + state->c * t + d; f32 f_minus_x = f - x; if ( fabsf( f_minus_x ) < epsilon ) { break; } l = ( f > x ) ? l : t; } return t; } static f32 interpolation_cubic_t_from_x_bisection_simd_4( cubic_state_t* state, f32 x, f32 epsilon ) { f32 t = 0; f32 d = state->start; __m128 as = _mm_set1_ps( state->a ); __m128 bs = _mm_set1_ps( state->b ); __m128 cs = _mm_set1_ps( state->c ); __m128 ds = _mm_set1_ps( d ); __m128 epsilons = _mm_set_ps1( epsilon ); __m128 xs = _mm_set_ps1( x ); __m128 ls = _mm_set_ps1( 0 ); __m128 factors = _mm_set_ps( 0.8f, 0.6f, 0.4f, 0.2f ); __m128 distance = _mm_set_ps1( 1.0f ); __m128 scale = _mm_set_ps1( 0.2f ); while ( 1 ) { __m128 ws = _mm_mul_ps( distance, factors ); __m128 ts = _mm_add_ps( ls, ws ); /* NOTE simon (26/04/25 13:57:04): Compute cubic. */ __m128 tt = _mm_mul_ps( ts, ts ); __m128 ttt = _mm_mul_ps( tt, ts ); __m128 attt = _mm_mul_ps( as, ttt ); __m128 btt = _mm_mul_ps( bs, tt ); __m128 ct = _mm_mul_ps( cs, ts ); __m128 f1 = _mm_add_ps( attt, btt ); __m128 f2 = _mm_add_ps( ct, ds ); __m128 f = _mm_add_ps( f1, f2 ); __m128 f_minus_x = _mm_sub_ps( f, xs ); /* NOTE simon (26/04/25 13:57:12): Absolute value */ __m128 n_zero = _mm_set1_ps( -0.0f ); __m128 abs = _mm_andnot_ps( n_zero, f_minus_x ); /* NOTE simon (26/04/25 13:57:37): Have we found the result ? */ __m128 cmp = _mm_cmple_ps( abs, epsilons ); u32 mask = _mm_movemask_ps( cmp ); if ( mask ) { /* NOTE simon (26/04/25 13:58:08): Extract the result. */ __m128 shuffled; if ( mask & 0x1 ) { shuffled = _mm_shuffle_ps( ts, ts, 0 ); } else if ( mask & 0x2 ) { shuffled = _mm_shuffle_ps( ts, ts, 1 ); } else if ( mask & 0x4 ) { shuffled = _mm_shuffle_ps( ts, ts, 2 ); } else { // _assert( mask & 0x8 ); shuffled = _mm_shuffle_ps( ts, ts, 3 ); } t = _mm_cvtss_f32( shuffled ); break; } /* NOTE simon (26/04/25 13:58:20): Update left and distance. */ __m128 setl = _mm_cmplt_ps( f, xs ); u32 lmask = _mm_movemask_ps( setl ); if ( lmask & 0x8 ) { // 0b11111111 ls = _mm_shuffle_ps( ts, ts, 0xff ); } else if ( lmask & 0x4 ) { // 0b10101010 ls = _mm_shuffle_ps( ts, ts, 0xaa ); } else if ( lmask & 0x2 ) { // 0b01010101 ls = _mm_shuffle_ps( ts, ts, 0x55 ); } else if ( lmask & 0x1 ) { // 0b00000000 ls = _mm_shuffle_ps( ts, ts, 0 ); } distance = _mm_mul_ps( distance, scale ); } return t; } static f32 interpolation_cubic_t_from_x_bisection_simd_5( cubic_state_t* state, f32 x, f32 epsilon ) { f32 t = 0; f32 d = state->start; __m128 as = _mm_set1_ps( state->a ); __m128 bs = _mm_set1_ps( state->b ); __m128 cs = _mm_set1_ps( state->c ); __m128 ds = _mm_set1_ps( d ); __m128 epsilons = _mm_set_ps1( epsilon ); __m128 xs = _mm_set_ps1( x ); __m128 ls = _mm_set_ps1( 0 ); __m128 factors = _mm_set_ps( 0.8f, 0.6f, 0.4f, 0.2f ); __m128 distance = _mm_set_ps1( 1.0f ); __m128 scale = _mm_set_ps1( 0.2f ); while ( 1 ) { __m128 ws = _mm_mul_ps( distance, factors ); __m128 ts = _mm_add_ps( ls, ws ); /* NOTE simon (26/04/25 13:57:04): Compute cubic. */ __m128 tt = _mm_mul_ps( ts, ts ); __m128 ttt = _mm_mul_ps( tt, ts ); __m128 attt = _mm_mul_ps( as, ttt ); __m128 btt = _mm_mul_ps( bs, tt ); __m128 ct = _mm_mul_ps( cs, ts ); __m128 f1 = _mm_add_ps( attt, btt ); __m128 f2 = _mm_add_ps( ct, ds ); __m128 f = _mm_add_ps( f1, f2 ); __m128 f_minus_x = _mm_sub_ps( f, xs ); /* NOTE simon (26/04/25 13:57:12): Absolute value */ __m128 n_zero = _mm_set1_ps( -0.0f ); __m128 abs = _mm_andnot_ps( n_zero, f_minus_x ); /* NOTE simon (26/04/25 13:57:37): Have we found the result ? */ __m128 cmp = _mm_cmple_ps( abs, epsilons ); u32 mask = _mm_movemask_ps( cmp ); if ( mask ) { /* NOTE simon (26/04/25 13:58:08): Extract the result. */ __m128i shuffle_mask = *( __m128i* ) ( shuffle_lut_2 + ( mask * 16 ) ); __m128i shuffled = _mm_shuffle_epi8( _mm_castps_si128( ts ), shuffle_mask ); t = _mm_cvtss_f32( _mm_castsi128_ps( shuffled ) ); break; } /* NOTE simon (26/04/25 13:58:20): Update left and distance. */ __m128 setl = _mm_cmplt_ps( f, xs ); u32 lmask = _mm_movemask_ps( setl ); { __m128i shuffle_mask = *( __m128i* ) ( shuffle_lut_2 + ( lmask * 16 ) ); /* NOTE simon (29/04/25 16:20:07): Storing in a variable isn't great. The compiler can still do the right thing and use the suffle with the memory operand, but it's not guaranteed. */ __m128i shuffled = _mm_shuffle_epi8( _mm_castps_si128( ts ), shuffle_mask ); __m128i l_blend_mask = _mm_set1_epi8( lmask ? 0xff : 0 ); ls = _mm_blendv_ps( ls, _mm_castsi128_ps( shuffled ), _mm_castsi128_ps( l_blend_mask ) ); } distance = _mm_mul_ps( distance, scale ); } return t; } static f32 interpolation_cubic_t_from_x_bisection_simd_6( cubic_state_t* state, f32 x, f32 epsilon ) { f32 t = 0; f32 d = state->start; __m128 as = _mm_set1_ps( state->a ); __m128 bs = _mm_set1_ps( state->b ); __m128 cs = _mm_set1_ps( state->c ); __m128 ds = _mm_set1_ps( d ); // __m128 epsilons = _mm_set_ps1( epsilon ); __m128 xs = _mm_set_ps1( x ); __m128 ls = _mm_set_ps1( 0 ); __m128 factors = _mm_set_ps( 0.8f, 0.6f, 0.4f, 0.2f ); __m128 distance = _mm_set_ps1( 1.0f ); __m128 scale = _mm_set_ps1( 0.2f ); __m128 ts; __m128 f; for ( u32 i = 0; i < 9; i++ ) { __m128 ws = _mm_mul_ps( distance, factors ); ts = _mm_add_ps( ls, ws ); /* NOTE simon (26/04/25 13:57:04): Compute cubic. */ __m128 tt = _mm_mul_ps( ts, ts ); __m128 ttt = _mm_mul_ps( tt, ts ); __m128 attt = _mm_mul_ps( as, ttt ); __m128 btt = _mm_mul_ps( bs, tt ); __m128 ct = _mm_mul_ps( cs, ts ); __m128 f1 = _mm_add_ps( attt, btt ); __m128 f2 = _mm_add_ps( ct, ds ); f = _mm_add_ps( f1, f2 ); /* NOTE simon (26/04/25 13:58:20): Update left and distance. */ __m128 setl = _mm_cmplt_ps( f, xs ); u32 lmask = _mm_movemask_ps( setl ); { __m128i shuffle_mask = *( __m128i* ) ( shuffle_lut_2 + ( lmask * 16 ) ); /* NOTE simon (29/04/25 16:20:07): Storing in a variable isn't great. The compiler can still do the right thing and use the suffle with the memory operand, but it's not guaranteed. */ __m128i shuffled = _mm_shuffle_epi8( _mm_castps_si128( ts ), shuffle_mask ); __m128i l_blend_mask = _mm_set1_epi8( lmask ? 0xff : 0 ); ls = _mm_blendv_ps( ls, _mm_castsi128_ps( shuffled ), _mm_castsi128_ps( l_blend_mask ) ); } distance = _mm_mul_ps( distance, scale ); } #if 0 __m128 f_minus_x = _mm_sub_ps( f, xs ); /* NOTE simon (26/04/25 13:57:12): Absolute value */ __m128 n_zero = _mm_set1_ps( -0.0f ); __m128 abs = _mm_andnot_ps( n_zero, f_minus_x ); /* NOTE simon (26/04/25 13:57:37): Have we found the result ? */ __m128 cmp = _mm_cmple_ps( abs, epsilons ); u32 mask = _mm_movemask_ps( cmp ); /* NOTE simon (26/04/25 13:58:08): Extract the result. */ __m128 shuffled; if ( mask & 0x1 ) { shuffled = _mm_shuffle_ps( ts, ts, 0 ); } else if ( mask & 0x2 ) { shuffled = _mm_shuffle_ps( ts, ts, 1 ); } else if ( mask & 0x4 ) { shuffled = _mm_shuffle_ps( ts, ts, 2 ); } else { // _assert( mask & 0x8 ); shuffled = _mm_shuffle_ps( ts, ts, 3 ); } t = _mm_cvtss_f32( shuffled ); #else t = _mm_cvtss_f32( ts ); #endif return t; } #if 0 static f32 interpolation_cubic_t_from_x_NR_BI( cubic_state_t* state, f32 x, f32 guess, f32 epsilon ) { f32 a = ( -start + 3 * cp_1 - 3 * cp_2 + end ); f32 b = ( 3 * start - 6 * cp_1 + 3 * cp_2 ); f32 c = ( -3 * start + 3 * cp_1 ); f32 d = ( start - x ); f32 previous[ 2 ] = { FLT_MAX, FLT_MAX }; b32 do_bisection = 1; f32 bisection_l = 0; f32 bisection_r = 1; f32 t = guess; for ( umm i = 0; i < 10; i++ ) { f32 f = a * t * t * t + b * t * t + c * t + d; if ( fabsf( f ) < epsilon ) { do_bisection = 0; break; } f32 derivative = 3 * a * t * t + 2 * b * t + c; if ( derivative == 0 ) { break; } t = t - ( f / derivative ); if ( t == previous[ 0 ] || t == previous[ 1 ] ) { #if 0 if ( previous[ 0 ] >= 0 && previous[ 1 ] <= 1 ) { if ( previous[ 0 ] < previous[ 1 ] ) { bisection_l = previous[ 0 ]; bisection_r = previous[ 1 ]; } else { bisection_l = previous[ 1 ]; bisection_r = previous[ 0 ]; } guess = interpolation_linear( bisection_l, bisection_r, 0.5f ); } #endif break; } previous[ 0 ] = previous[ 1 ]; previous[ 1 ] = t; } if ( t < 0 || t > 1 || do_bisection ) { t = guess; while ( bisection_l < bisection_r ) { f32 f = a * t * t * t + b * t * t + c * t + start; f32 f_minus_x = f - x; if ( fabsf( f_minus_x ) < epsilon ) { break; } if ( f - d > x ) { bisection_r = t; } else { bisection_l = t; } t = interpolation_linear( bisection_l, bisection_r, 0.5f ); } } return t; } #elif 1 static f32 interpolation_cubic_t_from_x_NR_BI( cubic_state_t* state, f32 x, f32 guess, f32 epsilon ) { f32 d = ( state->start - x ); b32 do_bisection = 1; f32 bisection_l = 0; f32 bisection_r = 1; f32 t = guess; for ( umm i = 0; i < 10; i++ ) { f32 f = state->a * t * t * t + state->b * t * t + state->c * t + d; if ( fabsf( f ) < epsilon ) { do_bisection = 0; break; } f32 derivative = 3 * state->a * t * t + 2 * state->b * t + state->c; if ( derivative == 0 ) { break; } t = t - ( f / derivative ); } if ( t < 0 || t > 1 || do_bisection ) { #if 0 t = guess; while ( bisection_l < bisection_r ) { f32 f = state->a * t * t * t + state->b * t * t + state->c * t + state->start; f32 f_minus_x = f - x; if ( fabsf( f_minus_x ) < epsilon ) { break; } if ( f > x ) { bisection_r = t; } else { bisection_l = t; } t = bisection_l * 0.5f + bisection_r * 0.5f; } #elif 1 t = 0; f32 l = 0; f32 w = 1; while ( 1 ) { w *= 0.5f; t = l + w; f32 f = state->a * t * t * t + state->b * t * t + state->c * t + state->start; f32 f_minus_x = f - x; if ( fabsf( f_minus_x ) < epsilon ) { break; } l = ( f > x ) ? l : t; } #else t = interpolation_cubic_t_from_x_bisection_simd_5( state, x, epsilon ); #endif } return t; } #else /* NOTE simon (07/05/25 15:03:52): The generated code by MSVC doesn't look goods, lots of unpck ? */ static f32 interpolation_cubic_t_from_x_NR_BI( cubic_state_t* state, f32 x, f32 guess, f32 epsilon ) { f32 d = ( state->start - x ); b32 do_bisection = 1; f32 bisection_l = 0; f32 bisection_r = 1; f32 t = guess; __m128 ts = _mm_set_ps1( t ); __m128 as = _mm_set_ps1( state->a ); __m128 bs = _mm_set_ps1( state->b ); __m128 cs = _mm_set_ps1( state->c ); __m128 abs_mask = _mm_castsi128_ps( _mm_set_epi32( 0, 0, 0x80000000, 0x80000000 ) ); __m128 epsilons = _mm_set_ps1( epsilon ); for ( umm i = 0; i < 10; i++ ) { #if 0 f32 f = state->a * t * t * t + state->b * t * t + state->c * t + d; if ( fabsf( f ) < epsilon ) { do_bisection = 0; break; } f32 derivative = 3 * state->a * t * t + 2 * state->b * t + state->c; if ( derivative == 0 ) { break; } t = t - ( f / derivative ); #else /* NOTE simon (06/05/25 18:29:22): Trying to compute function and derivative at the same time. */ __m128 y = _mm_mul_ps( as, ts ); // at | at y = _mm_mul_ps( y, ts ); // att | att __m128 t_3 = _mm_blend_ps( ts, _mm_set_ps1( 3 ), 0x2 ); // t | 3 y = _mm_mul_ps( y, t_3 ); // attt | 3att __m128 z = _mm_mul_ps( bs, ts ); // bt | bt __m128 t_2 = _mm_blend_ps( ts, _mm_set_ps1( 2 ), 0x2 ); // t | 2 z = _mm_mul_ps( z, t_2 ); // btt | 2bt __m128 t_1 = _mm_blend_ps( ts, _mm_set_ps1( 1 ), 0x2 ); // t | 1 __m128 w = _mm_mul_ps( cs, t_1 ); // ct | c __m128 f1 = _mm_add_ps( y, z ); __m128 f2 = _mm_add_ps( w, _mm_set_ps( 0, 0, 0, d ) ); __m128 f = _mm_add_ps( f1, f2 ); __m128 abs = _mm_andnot_ps( abs_mask, f ); __m128 cmp = _mm_cmple_ps( abs, epsilons ); u32 mask = _mm_movemask_ps( cmp ); if ( mask & 0x1 ) { t = _mm_cvtss_f32( ts ); do_bisection = 0; break; } #if 0 /* NOTE simon (07/05/25 16:01:11): Derivative is close to zero. Not doing this test is faster on average. */ else if ( mask & 0x2 ) { break; } #endif __m128 a = _mm_shuffle_ps( f, f, 0x0 ); __m128 b = _mm_shuffle_ps( f, f, 0x5 ); __m128 div = _mm_div_ps( a, b ); ts = _mm_sub_ps( ts, div ); #endif } if ( t < 0 || t > 1 || do_bisection ) { #if 0 t = guess; while ( bisection_l < bisection_r ) { f32 f = state->a * t * t * t + state->b * t * t + state->c * t + state->start; f32 f_minus_x = f - x; if ( fabsf( f_minus_x ) < epsilon ) { break; } if ( f > x ) { bisection_r = t; } else { bisection_l = t; } t = bisection_l * 0.5f + bisection_r * 0.5f; } #elif 1 t = 0; f32 l = 0; f32 w = 1; while ( 1 ) { w *= 0.5f; t = l + w; f32 f = state->a * t * t * t + state->b * t * t + state->c * t + state->start; f32 f_minus_x = f - x; if ( fabsf( f_minus_x ) < epsilon ) { break; } l = ( f > x ) ? l : t; } #else t = interpolation_cubic_t_from_x_bisection_simd_5( state, x, epsilon ); #endif } return t; } #endif static f32 interpolation_cubic_t_from_x_NR_BI_SIMD( cubic_state_t* state, f32 x, f32 epsilon ) { f32 d = ( state->start - x ); __m128 as = _mm_set1_ps( state->a ); __m128 bs = _mm_set1_ps( state->b ); __m128 cs = _mm_set1_ps( state->c ); __m128 ds = _mm_set1_ps( d ); __m128 epsilons = _mm_set_ps1( epsilon ); // __m128 zeros = _mm_set_ps1( 0 ); // __m128 ones = _mm_set_ps1( 1 ); __m128 ts = _mm_set_ps( 0.8f, 0.6f, 0.4f, 0.2f ); b32 do_bisection = 1; f32 t = 0; u32 active_lanes = 0xf; for ( umm i = 0; i < 10; i++ ) { // f32 f = a * t * t * t + b * t * t + c * t + d; __m128 tt = _mm_mul_ps( ts, ts ); __m128 at = _mm_mul_ps( as, ts ); __m128 attt = _mm_mul_ps( at, tt ); __m128 bt = _mm_mul_ps( bs, ts ); __m128 btt = _mm_mul_ps( bt, ts ); __m128 ct = _mm_mul_ps( cs, ts ); __m128 f1 = _mm_add_ps( attt, btt ); __m128 f2 = _mm_add_ps( ct, ds ); __m128 fr = _mm_add_ps( f1, f2 ); /* NOTE simon (26/04/25 13:57:12): Absolute value */ __m128 n_zero = _mm_set1_ps( -0.0f ); __m128 abs = _mm_andnot_ps( n_zero, fr ); /* NOTE simon (26/04/25 13:57:37): Have we found the result ? */ __m128 cmp = _mm_cmple_ps( abs, epsilons ); u32 mask = _mm_movemask_ps( cmp ); mask &= active_lanes; if ( mask ) { /* NOTE simon (05/05/25 14:42:12): We can converge outside the 0 to 1 range, so verify that the result is between 0 and 1 */ __m128 gez = _mm_cmpge_ps( ts, _mm_set_ps1( 0 ) ); __m128 leo = _mm_cmple_ps( ts, _mm_set_ps1( 1 ) ); __m128 valid = _mm_and_ps( gez, leo ); u32 valid_mask = _mm_movemask_ps( valid ); u32 extract = mask & valid_mask; if ( extract ) { /* NOTE simon (26/04/25 13:58:08): Extract the result. */ __m128i shuffle_mask = *( __m128i* ) ( shuffle_lut + ( extract * 16 ) ); __m128i shuffled = _mm_shuffle_epi8( _mm_castps_si128( ts ), shuffle_mask ); t = _mm_cvtss_f32( _mm_castsi128_ps( shuffled ) ); do_bisection = 0; break; } active_lanes ^= ( mask & ~valid_mask ); } /* NOTE simon (06/05/25 15:16:01): f32 derivative = 3 * a * t * t + 2 * b * t + c; */ __m128 t3 = _mm_mul_ps( ts, _mm_set_ps1( 3 ) ); __m128 att3 = _mm_mul_ps( at, t3 ); __m128 bt2 = _mm_mul_ps( bt, _mm_set_ps1( 2 ) ); __m128 derivative = _mm_add_ps( att3, bt2 ); derivative = _mm_add_ps( derivative, cs ); /* NOTE simon (02/05/25 16:22:59): Can produce inf as we might divide by zero. */ __m128 divs = _mm_div_ps( fr, derivative ); ts = _mm_sub_ps( ts, divs ); } if ( do_bisection ) { t = interpolation_cubic_t_from_x_bisection_simd_5( state, x, epsilon ); } return t; } int main( int argc, char** argv ) { f32 precision = 0.00001f; u32 iteration_count = 100; f32 increment = 1.0f / cast( f32, iteration_count ); u32 branch_prediction_test = 1; #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < branch_prediction_test; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_bisection_simd_1( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "left pack total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < branch_prediction_test; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_bisection_simd_2( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "ifs total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < branch_prediction_test; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_bisection_simd_3( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "ifs 2 total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < branch_prediction_test; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_bisection( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "left dist total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < branch_prediction_test; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_bisection_simd_4( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "simd4 total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < branch_prediction_test; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_bisection_simd_5( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "simd5 total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < branch_prediction_test; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_bisection_simd_6( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "simd6 total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); f32 last = 0; for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < 1; l++ ) { #if 1 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); last = 0.5f; // f32 x = k * increment; // f32 guess = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_NR_BI( &state, x, last, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); last = t; // printf( "%lld\n", end - start ); total += ( end - start ); hit++; } } } } printf( "NR BI total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif #if 1 { random_32_t rng = random_32_make( 1234, 5678 ); u64 total = 0; u64 hit = 0; for ( u32 i = 0; i < iteration_count + 1; i++ ) { for ( u32 j = 0; j < iteration_count + 1; j++ ) { #if 0 f32 cp_1 = i * increment; f32 cp_2 = j * increment; #else f32 cp_1 = random_32_get_f32( &rng ); f32 cp_2 = random_32_get_f32( &rng ); #endif cubic_state_t state; cubic_state_initialize( &state, 0, cp_1, cp_2, 1 ); for ( u32 k = 0; k < iteration_count + 1; k++ ) { for ( u32 l = 0; l < 1; l++ ) { #if 0 f32 x = k * increment; #else f32 x = random_32_get_f32( &rng ); #endif u64 start = __rdtsc( ); f32 t = interpolation_cubic_t_from_x_NR_BI_SIMD( &state, x, precision * 0.75f ); u64 end = __rdtsc( ); f32 check = interpolation_cubic( 0, cp_1, cp_2, 1, t ); _assert( fabsf( check - x ) < precision ); total += ( end - start ); hit++; } } } } printf( "NR BI SIMD total: %lld | hit: %lld | avg: %lld\n", total, hit, total / hit ); } #endif return 0; }