Sine performance test

sin_test.jai

/*
   Time LUT (512): 181ms error: 798.972624
   Time LUT (1024): 181ms error: 199.799809
   Time LUT (4096): 181ms error: 13.382117
   Time LUT (16MB): 1279ms error: 0.540578
   Time Root Factored 13: 141ms error: 33.625392
*/

#import "Basic";
#import "Math";
#import "Random";

main :: () {
    iterations :: 100_000_000;

    sine_lut_512  :: #bake_arguments sin_lookup_table(size = 512);
    sine_lut_1024 :: #bake_arguments sin_lookup_table(size = 1024);
    sine_lut_4096 :: #bake_arguments sin_lookup_table(size = 4096);
    sine_lut_huge :: #bake_arguments sin_lookup_table(size = 0x100_0000);

    result_lut_512  := test_performance(sine_lut_512,         iterations);
    result_lut_1024 := test_performance(sine_lut_1024,        iterations);
    result_lut_4096 := test_performance(sine_lut_4096,        iterations);
    result_lut_huge := test_performance(sine_lut_huge,        iterations);
    result_poly     := test_performance(sin_root_factored_13, iterations);

    report_test_result(result_lut_512,  "LUT (512)");
    report_test_result(result_lut_1024, "LUT (1024)");
    report_test_result(result_lut_4096, "LUT (4096)");
    report_test_result(result_lut_huge, "LUT (16MB)");
    report_test_result(result_poly,     "Root Factored 13");
}

test_performance :: ($proc: (phase: float) -> float, iterations: int) -> Test_Result {
    phase_values := NewArray(iterations, float);
    ground_truth := NewArray(iterations, float);
    results      := NewArray(iterations, float);

    for 0..iterations-1 {
        phase_values[it] = cast(float) random_get_zero_to_one_new() - 0.5;
        ground_truth[it] = sin(phase_values[it] * TAU);
    }

    start_time := current_time_monotonic();
    for 0..iterations-1 {
        results[it] = inline proc(phase_values[it]);
    }
    end_time   := current_time_monotonic();

    error: float64;
    for 0..iterations-1 {
        error += abs(cast(float64) results[it] - cast(float64) ground_truth[it]);
    }

    return {end_time - start_time, error};
}

sin_lookup_table :: (signed_phase: float, $size: int) -> float {
    #assert is_power_of_two(size);

    table :: #run -> [] float {
        result := NewArray(size, float);

        for *result {
            phase := cast(float) it_index / size;
            signed_phase := phase - 0.5;

            it.* = sin(TAU * signed_phase);
        }

        return result;
    }

    scaled_phase := (signed_phase + 0.5) * size;

    low  := cast(int)  scaled_phase;
    high := cast(int) (scaled_phase + 1);

    high &= size - 1;
    a    := scaled_phase - cast(float) low;

    return table[low] * (1.0 - a) + table[high] * a;
}

sin_root_factored_13 :: (signed_phase: float) -> float {
    phase2 := signed_phase * signed_phase;

    fma :: (a: float, b: float, c: float) -> float #expand {
        return a * b + c;
    }

    polynomial := fma(cast(float) 3.1616015434265137, phase2, cast(float) -14.09662590026855);
    polynomial  = fma(polynomial, phase2, cast(float)  38.495872497558594);
    polynomial  = fma(polynomial, phase2, cast(float) -67.07662200927734);
    polynomial  = fma(polynomial, phase2, cast(float)  64.83582305908203);
    polynomial  = fma(polynomial, phase2, cast(float) -25.132740020751953);

    return signed_phase * (phase2 - 0.25) * polynomial;
}

Test_Result :: struct {
    time_taken: Apollo_Time;
    error:      float64;
}

report_test_result :: (result: Test_Result, name: string) {
    log("Time %: %ms error: %", name, to_milliseconds(result.time_taken), result.error);
}

is_power_of_two :: (x: int) -> bool {
    return x & (x-1) == 0;
}



/*

Disassembly of sin_root_factored13 when inlined in the loop:

  222a20:	f3 41 0f 10 3e       	movss  (%r14),%xmm7

  222a25:	44 0f 28 c7          	movaps %xmm7,%xmm8
  222a29:	f3 44 0f 59 c7       	mulss  %xmm7,%xmm8
  222a2e:	45 0f 28 c8          	movaps %xmm8,%xmm9
  222a32:	f3 44 0f 59 c8       	mulss  %xmm0,%xmm9
  222a37:	f3 44 0f 58 c9       	addss  %xmm1,%xmm9
  222a3c:	f3 45 0f 59 c8       	mulss  %xmm8,%xmm9
  222a41:	f3 44 0f 58 ca       	addss  %xmm2,%xmm9
  222a46:	f3 45 0f 59 c8       	mulss  %xmm8,%xmm9
  222a4b:	f3 44 0f 58 cb       	addss  %xmm3,%xmm9
  222a50:	f3 45 0f 59 c8       	mulss  %xmm8,%xmm9
  222a55:	f3 44 0f 58 cc       	addss  %xmm4,%xmm9
  222a5a:	f3 45 0f 59 c8       	mulss  %xmm8,%xmm9
  222a5f:	f3 44 0f 58 cd       	addss  %xmm5,%xmm9
  222a64:	f3 44 0f 58 c6       	addss  %xmm6,%xmm8
  222a69:	f3 44 0f 59 c7       	mulss  %xmm7,%xmm8
  222a6e:	f3 45 0f 59 c1       	mulss  %xmm9,%xmm8
  222a73:	f3 44 0f 11 00       	movss  %xmm8,(%rax)
  222a78:	48 83 c0 04          	add    $0x4,%rax
  222a7c:	49 83 c6 04          	add    $0x4,%r14
  222a80:	49 ff cc             	dec    %r12
  222a83:	75 9b                	jne    222a20 <test_performance_1700000030+0x1b0>

Disassembly of sin_lookup_table when inlined in the loop:

  226f50:	f3 41 0f 10 1e       	movss  (%r14),%xmm3

  226f55:	f3 0f 58 d8          	addss  %xmm0,%xmm3
  226f59:	f3 0f 59 d9          	mulss  %xmm1,%xmm3
  226f5d:	f3 48 0f 2c cb       	cvttss2si %xmm3,%rcx

  226f62:	0f 28 e3             	movaps %xmm3,%xmm4
  226f65:	f3 0f 58 e2          	addss  %xmm2,%xmm4
  226f69:	f3 48 0f 2c d4       	cvttss2si %xmm4,%rdx

  226f6e:	0f 57 e4             	xorps  %xmm4,%xmm4
  226f71:	f3 48 0f 2a e1       	cvtsi2ss %rcx,%xmm4
  226f76:	f3 0f 5c dc          	subss  %xmm4,%xmm3
  226f7a:	0f 28 e2             	movaps %xmm2,%xmm4
  226f7d:	f3 0f 5c e3          	subss  %xmm3,%xmm4
  226f81:	f3 0f 59 24 8d d8 a3 	mulss  0x20a3d8(,%rcx,4),%xmm4
  226f88:	20 00
  226f8a:	81 e2 ff 01 00 00    	and    $0x1ff,%edx
  226f90:	f3 0f 59 1c 95 d8 a3 	mulss  0x20a3d8(,%rdx,4),%xmm3
  226f97:	20 00
  226f99:	f3 0f 58 dc          	addss  %xmm4,%xmm3
  226f9d:	f3 0f 11 18          	movss  %xmm3,(%rax)
  226fa1:	48 83 c0 04          	add    $0x4,%rax
  226fa5:	49 83 c6 04          	add    $0x4,%r14
  226fa9:	49 ff cc             	dec    %r12
  226fac:	75 a2                	jne    226f50 <test_performance_1700000023+0x190>

Each loop has been 4-wide autovectorised by LLVM, but without utilizing specialised
instruction sets.

Both of these could be hand-vectorised to be most-likely faster - the FMA one using
Intel AVX FMA instructions, the lookup table using AVX2 Gather instructions.

*/