Benchmarking interpolation functions for time domain signals.

interpolator_bench.c

/*
    James William Fletcher (github.com/mrbid)
        September 2021 - March 2024

    Benchmarking interpolation functions for time domain signals, such as audio.
    https://james-william-fletcher.medium.com/benchmarking-interpolation-functions-9f84de397109

    references of interest:
    http://www.ee.ic.ac.uk/pcheung/teaching/ee3_Study_Project/Sinewave%20Generation(708).pdf
    https://demonstrations.wolfram.com/SineWaveGenerationUsingAnUnstableIIRFilter/
    
    notes: All of the Olli Niemitalo interpolators will compile with double to float conversion ops in C.

    compile: gcc interpolator_bench.c -Ofast -lm -o interpolator_bench
*/
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <string.h>
#include <locale.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>
#include <errno.h>
#include <fcntl.h>
#include <x86intrin.h>

#define interval 32000000 // microseconds

// ----------------------------------------------------------
// https://www.musicdsp.org/en/latest/Other/93-hermite-interpollation.html
// 2002-01-17 03:14:54: Dave from Muon Software, originally from Josh Scholar

float order3Spline(const float x, const float L1, const float L0, const float H0, const float H1)
{
    return
    L0 +
    .5f*
    x*(H0-L1 +
       x*(H0 + L0*(-2.f) + L1 +
          x*( (H0 - L0)*9.f + (L1 - H1)*3.f +
             x*((L0 - H0)*15.f + (H1 -  L1)*5.f +
                x*((H0 - L0)*6.f + (L1 - H1)*2.f )))));
}

// ----------------------------------------------------------
// https://www.musicdsp.org/en/latest/Other/93-hermite-interpollation.html
// 2002-04-09 16:55:35

// original
float hermite1(float x, float y0, float y1, float y2, float y3)
{
	// 4-point, 3rd-order Hermite (x-form)
	float c0 = y1;
	float c1 = 0.5f * (y2 - y0);
	float c2 = y0 - 2.5f * y1 + 2.f * y2 - 0.5f * y3;
	float c3 = 1.5f * (y1 - y2) + 0.5f * (y3 - y0);

	return ((c3 * x + c2) * x + c1) * x + c0;
}

// james mccartney
float hermite2(float x, float y0, float y1, float y2, float y3)
{
	// 4-point, 3rd-order Hermite (x-form)
	float c0 = y1;
	float c1 = 0.5f * (y2 - y0);
	float c3 = 1.5f * (y1 - y2) + 0.5f * (y3 - y0);
	float c2 = y0 - y1 + c1 - c3;

	return ((c3 * x + c2) * x + c1) * x + c0;
}

// james mccartney
float hermite3(float x, float y0, float y1, float y2, float y3)
{
	// 4-point, 3rd-order Hermite (x-form)
	float c0 = y1;
	float c1 = 0.5f * (y2 - y0);
	float y0my1 = y0 - y1;
	float c3 = (y1 - y2) + 0.5f * (y3 - y0my1 - y2);
	float c2 = y0my1 + c1 - c3;

	return ((c3 * x + c2) * x + c1) * x + c0;
}

// laurent de soras
float hermite4(float frac_pos, float xm1, float x0, float x1, float x2)
{
	const float c = (x1 - xm1) * 0.5f;
	const float v = x0 - x1;
	const float w = c + v;
	const float a = w + v + (x2 - x0) * 0.5f;
	const float b_neg = w + a;

	return ((((a * frac_pos) - b_neg) * frac_pos + c) * frac_pos + x0);
}

// ----------------------------------------------------------
// Paul Bourke (December 1999)
// http://paulbourke.net/miscellaneous/interpolation/

// Each parameter begins with an 'm' or 'p' apart from the current sample which
// begins with 's', this means all parameters beginning with 'm' are samples into
// the past (-) before the current sample, and all samples beginning with 'p'
// it are samples into the future (+).

float LinearInterpolate(float a, float b, float i)
{
	return (b - a) * i + a;
}

float CosineInterpolate(float s1, float p2, float mu)
{
	const float mu2 = (1.f - cosf(mu * 3.141592741f)) * 0.5f;
	return s1 * (1.f - mu2) + p2 * mu2;
}

float CubicInterpolate(float m0, float s1, float p2, float p3, float mu)
{
	const float mu2 = mu * mu;
	const float a0 = (p3 - p2 - m0 + s1);
	return a0 * mu * mu2 + (m0 - s1 - a0) * mu2 + (p2 - m0) * mu + s1;
}

// coefficients by Paul Breeuwsma (paulinternet.nl)
float CatmullInterpolate(float m0, float s1, float p2, float p3, float mu)
{
	const float mu2 = mu * mu;
	const float a0 = -0.5f * m0 + 1.5f * s1 - 1.5f * p2 + 0.5f * p3;
	const float a1 = m0 - 2.5f * s1 + 2.f * p2 - 0.5f * p3;
	const float a2 = -0.5f * m0 + 0.5f * p2;
	const float a3 = s1;
	return a0 * mu * mu2 + a1 * mu2 + a2 * mu + a3;
}

float HermiteInterpolate(float m0, float s1, float p2, float p3, float mu, float tension, float bias)
{
	float mv0, mv1, mu2, mu3;
	float a0, a1, a2, a3;

	mu2 = mu * mu;
	mu3 = mu2 * mu;
	mv0 = (s1 - m0) * (1.f + bias) * (1.f - tension) / 2.f;
	mv0 += (p2 - s1) * (1.f - bias) * (1.f - tension) / 2.f;
	mv1 = (p2 - s1) * (1.f + bias) * (1.f - tension) / 2.f;
	mv1 += (p3 - p2) * (1.f - bias) * (1.f - tension) / 2.f;
	a0 = 2.f * mu3 - 3.f * mu2 + 1.f;
	a1 = mu3 - 2.f * mu2 + mu;
	a2 = mu3 - mu2;
	a3 = -2.f * mu3 + 3.f * mu2;

	return (a0 * s1 + a1 * mv0 + a2 * mv1 + a3 * p2);
}

// ----------------------------------------------------------
// Olli Niemitalo (October 2001)
// http://yehar.com/blog/wp-content/uploads/2009/08/deip.pdf

// I think the z-form uses Horner’s rule in an attempt to reduce multiplications.
// https://en.wikipedia.org/wiki/Horner%27s_method
/*
	6. Implementation 36

	To choose between the x- and the z-form, we find out the number of operations
	in both. First we find out which needs more multiplications. If multiplications are
	expensive on the target platform, then that is what matters. However, if multiplica-
	tions are as cheap as addition and subtraction operations, then the total number
	of operations counts. Optimized C language source code for the x- and z-form
	algorithms is given in the following, and the total number of operations is counted
	from each.
*/

// 4-point, 3rd-order B-spline (x-form)
float bspline_43x(const float *y, const float x)
{
	float ym1py1 = y[-1] + y[1];
	float c0 = 1 / 6.0 * ym1py1 + 2 / 3.0 * y[0];
	float c1 = 1 / 2.0 * (y[1] - y[-1]);
	float c2 = 1 / 2.0 * ym1py1 - y[0];
	float c3 = 1 / 2.0 * (y[0] - y[1]) + 1 / 6.0 * (y[2] - y[-1]);
	return ((c3 * x + c2) * x + c1) * x + c0;
}

// 4-point, 3rd-order B-spline (z-form)
float bspline_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-1] + y[2], modd1 = y[2] - y[-1];
	float even2 = y[0] + y[1], modd2 = y[1] - y[0];
	float c0 = 1 / 48.0 * even1 + 23 / 48.0 * even2;
	float c1 = 1 / 8.0 * modd1 + 5 / 8.0 * modd2;
	float c2 = 1 / 4.0 * (even1 - even2);
	float c3 = 1 / 6.0 * modd1 - 1 / 2.0 * modd2;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// 6-point, 5th-order B-spline (x-form)
float bspline_65x(const float *y, const float x)
{
	float ym2py2 = y[-2] + y[2], ym1py1 = y[-1] + y[1];
	float y2mym2 = y[2] - y[-2], y1mym1 = y[1] - y[-1];
	float sixthym1py1 = 1 / 6.0 * ym1py1;
	float c0 = 1 / 120.0 * ym2py2 + 13 / 60.0 * ym1py1 + 11 / 20.0 * y[0];
	float c1 = 1 / 24.0 * y2mym2 + 5 / 12.0 * y1mym1;
	float c2 = 1 / 12.0 * ym2py2 + sixthym1py1 - 1 / 2.0 * y[0];
	float c3 = 1 / 12.0 * y2mym2 - 1 / 6.0 * y1mym1;
	float c4 = 1 / 24.0 * ym2py2 - sixthym1py1 + 1 / 4.0 * y[0];
	float c5 = 1 / 120.0 * (y[3] - y[-2]) + 1 / 24.0 * (y[-1] - y[2]) +
			   1 / 12.0 * (y[1] - y[0]);
	return ((((c5 * x + c4) * x + c3) * x + c2) * x + c1) * x + c0;
}

// 4-point, 3rd-order Lagrange (x-form)
float lagrange_43x(const float *y, const float x)
{
	float c0 = y[0];
	float c1 = y[1] - 1 / 3.0 * y[-1] - 1 / 2.0 * y[0] - 1 / 6.0 * y[2];
	float c2 = 1 / 2.0 * (y[-1] + y[1]) - y[0];
	float c3 = 1 / 6.0 * (y[2] - y[-1]) + 1 / 2.0 * (y[0] - y[1]);
	return ((c3 * x + c2) * x + c1) * x + c0;
}

// 4-point, 3rd-order Lagrange (z-form)
float lagrange_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-1] + y[2], odd1 = y[-1] - y[2];
	float even2 = y[0] + y[1], odd2 = y[0] - y[1];
	float c0 = 9 / 16.0 * even2 - 1 / 16.0 * even1;
	float c1 = 1 / 24.0 * odd1 - 9 / 8.0 * odd2;
	float c2 = 1 / 4.0 * (even1 - even2);
	float c3 = 1 / 2.0 * odd2 - 1 / 6.0 * odd1;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// 6-point, 5th-order Lagrange (x-form)
float lagrange_65x(const float *y, const float x)
{
	float ym1py1 = y[-1] + y[1];
	float twentyfourthym2py2 = 1 / 24.0 * (y[-2] + y[2]);
	float c0 = y[0];
	float c1 = 1 / 20.0 * y[-2] - 1 / 2.0 * y[-1] - 1 / 3.0 * y[0] + y[1] - 1 / 4.0 * y[2] + 1 / 30.0 * y[3];
	float c2 = 2 / 3.0 * ym1py1 - 5 / 4.0 * y[0] - twentyfourthym2py2;
	float c3 = 5 / 12.0 * y[0] - 7 / 12.0 * y[1] + 7 / 24.0 * y[2] - 1 / 24.0 * (y[-2] + y[-1] + y[3]);
	float c4 = 1 / 4.0 * y[0] - 1 / 6.0 * ym1py1 + twentyfourthym2py2;
	float c5 = 1 / 120.0 * (y[3] - y[-2]) + 1 / 24.0 * (y[-1] - y[2]) + 1 / 12.0 * (y[1] - y[0]);
	return ((((c5 * x + c4) * x + c3) * x + c2) * x + c1) * x + c0;
}

// 6-point, 5th-order Lagrange (z-form)
float lagrange_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-2] + y[3], odd1 = y[-2] - y[3];
	float even2 = y[-1] + y[2], odd2 = y[-1] - y[2];
	float even3 = y[0] + y[1], odd3 = y[0] - y[1];
	float c0 = 3 / 256.0 * even1 - 25 / 256.0 * even2 + 75 / 128.0 * even3;
	float c1 = 25 / 384.0 * odd2 - 75 / 64.0 * odd3 - 3 / 640.0 * odd1;
	float c2 = 13 / 32.0 * even2 - 17 / 48.0 * even3 - 5 / 96.0 * even1;
	float c3 = 1 / 48.0 * odd1 - 13 / 48.0 * odd2 + 17 / 24.0 * odd3;
	float c4 = 1 / 48.0 * even1 - 1 / 16.0 * even2 + 1 / 24.0 * even3;
	float c5 = 1 / 24.0 * odd2 - 1 / 12.0 * odd3 - 1 / 120.0 * odd1;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}

// 4-point, 3rd-order Hermite (x-form)
float hermite_43x(const float *y, const float x)
{
	float c0 = y[0];
	float c1 = 1 / 2.0 * (y[1] - y[-1]);
	float c2 = y[-1] - 5 / 2.0 * y[0] + 2 * y[1] - 1 / 2.0 * y[2];
	float c3 = 1 / 2.0 * (y[2] - y[-1]) + 3 / 2.0 * (y[0] - y[1]);
	return ((c3 * x + c2) * x + c1) * x + c0;
}

// 4-point, 3rd-order Hermite (z-form)
float hermite_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-1] + y[2], odd1 = y[-1] - y[2];
	float even2 = y[0] + y[1], odd2 = y[0] - y[1];
	float c0 = 9 / 16.0 * even2 - 1 / 16.0 * even1;
	float c1 = 1 / 8.0 * odd1 - 11 / 8.0 * odd2;
	float c2 = 1 / 4.0 * (even1 - even2);
	float c3 = 3 / 2.0 * odd2 - 1 / 2.0 * odd1;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// 6-point, 3rd-order Hermite (x-form)
float hermite_63x(const float *y, const float x)
{
	float c0 = y[0];
	float c1 = 1 / 12.0 * (y[-2] - y[2]) + 2 / 3.0 * (y[1] - y[-1]);
	float c2 = 5 / 4.0 * y[-1] - 7 / 3.0 * y[0] + 5 / 3.0 * y[1] - 1 / 2.0 * y[2] + 1 / 12.0 * y[3] - 1 / 6.0 * y[-2];
	float c3 = 1 / 12.0 * (y[-2] - y[3]) + 7 / 12.0 * (y[2] - y[-1]) + 4 / 3.0 * (y[0] - y[1]);
	return ((c3 * x + c2) * x + c1) * x + c0;
}

// 6-point, 3rd-order Hermite (z-form)
float hermite_63z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-2] + y[3], odd1 = y[-2] - y[3];
	float even2 = y[-1] + y[2], odd2 = y[-1] - y[2];
	float even3 = y[0] + y[1], fourthirdthodd3 = 4 / 3.0 * (y[0] - y[1]);
	float c0 = 1 / 96.0 * even1 - 3 / 32.0 * even2 + 7 / 12.0 * even3;
	float c1 = 7 / 48.0 * odd2 - fourthirdthodd3 - 1 / 48.0 * odd1;
	float c2 = 3 / 8.0 * even2 - 1 / 3.0 * even3 - 1 / 24.0 * even1;
	float c3 = 1 / 12.0 * odd1 - 7 / 12.0 * odd2 + fourthirdthodd3;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// 6-point, 5th-order Hermite (x-form)
float hermite_65x(const float *y, const float x)
{
	float eighthym2 = 1 / 8.0 * y[-2];
	float eleventwentyfourthy2 = 11 / 24.0 * y[2];
	float twelfthy3 = 1 / 12.0 * y[3];
	float c0 = y[0];
	float c1 = 1 / 12.0 * (y[-2] - y[2]) + 2 / 3.0 * (y[1] - y[-1]);
	float c2 = 13 / 12.0 * y[-1] - 25 / 12.0 * y[0] + 3 / 2.0 * y[1] - eleventwentyfourthy2 + twelfthy3 - eighthym2;
	float c3 = 5 / 12.0 * y[0] - 7 / 12.0 * y[1] + 7 / 24.0 * y[2] - 1 / 24.0 * (y[-2] + y[-1] + y[3]);
	float c4 = eighthym2 - 7 / 12.0 * y[-1] + 13 / 12.0 * y[0] - y[1] + eleventwentyfourthy2 - twelfthy3;
	float c5 = 1 / 24.0 * (y[3] - y[-2]) + 5 / 24.0 * (y[-1] - y[2]) + 5 / 12.0 * (y[1] - y[0]);
	return ((((c5 * x + c4) * x + c3) * x + c2) * x + c1) * x + c0;
}

// 6-point, 5th-order Hermite (z-form)
float hermite_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-2] + y[3], odd1 = y[-2] - y[3];
	float even2 = y[-1] + y[2], odd2 = y[-1] - y[2];
	float even3 = y[0] + y[1], odd3 = y[0] - y[1];
	float c0 = 3 / 256.0 * even1 - 25 / 256.0 * even2 + 75 / 128.0 * even3;
	float c1 = -3 / 128.0 * odd1 + 61 / 384.0 * odd2 - 87 / 64.0 * odd3;
	float c2 = -5 / 96.0 * even1 + 13 / 32.0 * even2 - 17 / 48.0 * even3;
	float c3 = 5 / 48.0 * odd1 - 11 / 16.0 * odd2 + 37 / 24.0 * odd3;
	float c4 = 1 / 48.0 * even1 - 1 / 16.0 * even2 + 1 / 24.0 * even3;
	float c5 = -1 / 24.0 * odd1 + 5 / 24.0 * odd2 - 5 / 12.0 * odd3;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}

// 4-point, 5th-order 2nd-order-osculating (x-form)
float osculating_45x(const float *y, const float x)
{
	float y1my0 = y[1] - y[0];
	float y2mym1 = y[2] - y[-1];
	float c0 = y[0];
	float c1 = 1 / 2.0 * (y[1] - y[-1]);
	float c2 = 1 / 2.0 * (y[-1] + y[1]) - y[0];
	float c3 = 9 / 2.0 * y1my0 - 3 / 2.0 * y2mym1;
	float c4 = 5 / 2.0 * y2mym1 - 15 / 2.0 * y1my0;
	float c5 = y[-1] + 3 * y1my0 - y[2];
	return ((((c5 * x + c4) * x + c3) * x + c2) * x + c1) * x + c0;
}

// 4-point, 5th-order 2nd-order-osculating (z-form)
float osculating_45z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-1] + y[2], odd1 = y[-1] - y[2];
	float even2 = y[0] + y[1], odd2 = y[0] - y[1];
	float c0 = 9 / 16.0 * even2 - 1 / 16.0 * even1;
	float c1 = 3 / 16.0 * odd1 - 25 / 16.0 * odd2;
	float c2 = 1 / 4.0 * (even1 - even2);
	float c5 = odd1 - 3 * odd2;
	return (((c5 * z * z - c5) * z + c2) * z + c1) * z + c0;
}

// 6-point, 5th-order 2nd-order-osculating (x-form)
float osculating_65x(const float *y, const float x)
{
	float c0 = y[0];
	float c1 = 1 / 12.0 * (y[-2] - y[2]) + 2 / 3.0 * (y[1] - y[-1]);
	float c2 = 2 / 3.0 * (y[-1] + y[1]) - 1 / 24.0 * (y[-2] + y[2]) - 5 / 4.0 * y[0];
	float c3 = 13 / 8.0 * y[-1] - 35 / 12.0 * y[0] + 11 / 4.0 * y[1] -
			   11 / 8.0 * y[2] + 7 / 24.0 * y[3] - 3 / 8.0 * y[-2];
	float c4 = 13 / 24.0 * y[-2] - 8 / 3.0 * y[-1] + 21 / 4.0 * y[0] -
			   31 / 6.0 * y[1] + 61 / 24.0 * y[2] - 1 / 2.0 * y[3];
	float c5 = 5 / 24.0 * (y[3] - y[-2]) + 25 / 24.0 * (y[-1] - y[2]) +
			   25 / 12.0 * (y[1] - y[0]);
	return ((((c5 * x + c4) * x + c3) * x + c2) * x + c1) * x + c0;
}

// 6-point, 5th-order 2nd-order-osculating (z-form)
float osculating_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-2] + y[3], odd1 = y[-2] - y[3];
	float even2 = y[-1] + y[2], odd2 = y[-1] - y[2];
	float even3 = y[0] + y[1], odd3 = y[0] - y[1];
	float c0 = 3 / 256.0 * even1 - 25 / 256.0 * even2 + 75 / 128.0 * even3;
	float c1 = 27 / 128.0 * odd2 - 281 / 192.0 * odd3 - 13 / 384.0 * odd1;
	float c2 = 13 / 32.0 * even2 - 17 / 48.0 * even3 - 5 / 96.0 * even1;
	float c3 = 3 / 16.0 * odd1 - 53 / 48.0 * odd2 + 19 / 8.0 * odd3;
	float c4 = 1 / 48.0 * even1 - 1 / 16.0 * even2 + 1 / 24.0 * even3;
	float c5 = 25 / 24.0 * odd2 - 25 / 12.0 * odd3 - 5 / 24.0 * odd1;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}

// 4-point, 2nd-order Watte tri-linear (x-form)
float watte_42x(const float *y, const float x)
{
	float ym1py2 = y[-1] + y[2];
	float c0 = y[0];
	float c1 = 3 / 2.0 * y[1] - 1 / 2.0 * (y[0] + ym1py2);
	float c2 = 1 / 2.0 * (ym1py2 - y[0] - y[1]);
	return (c2 * x + c1) * x + c0;
}

// 4-point, 2nd-order Watte tri-linear (z-form)
float watte_42z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-1] + y[2], even2 = y[0] + y[1];
	float c0 = 5 / 8.0 * even2 - 1 / 8.0 * even1;
	float c1 = y[1] - y[0];
	float c2 = 1 / 2.0 * (even1 - even2);
	return (c2 * z + c1) * z + c0;
}

// 4-point, 2nd-order parabolic 2x (x-form)
float parabolic_42x(const float *y, const float x)
{
	float y1mym1 = y[1] - y[-1];
	float c0 = 1 / 2.0 * y[0] + 1 / 4.0 * (y[-1] + y[1]);
	float c1 = 1 / 2.0 * y1mym1;
	float c2 = 1 / 4.0 * (y[2] - y[0] - y1mym1);
	return (c2 * x + c1) * x + c0;
}

// 4-point, 2nd-order parabolic 2x (z-form)
float parabolic_42z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[-1] + y[2], even2 = y[0] + y[1];
	float c0 = 1 / 16.0 * even1 + 7 / 16.0 * even2;
	float c1 = 1 / 4.0 * (y[1] - y[0] + y[2] - y[-1]);
	float c2 = 1 / 4.0 * (even1 - even2);
	return (c2 * z + c1) * z + c0;
}

// Optimal 2x (2-point, 3rd-order) (z-form)
float optimal2_23z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float c0 = even1 * 0.50037842517188658;
	float c1 = odd1 * 1.00621089801788210;
	float c2 = even1 * -0.004541102062639801;
	float c3 = odd1 * -1.57015627178718420;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 4x (2-point, 3rd-order) (z-form)
float optimal4_23z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float c0 = even1 * 0.50013034073688023;
	float c1 = odd1 * 1.09617817497678520;
	float c2 = even1 * -0.001564088842561871;
	float c3 = odd1 * -1.32598918957298410;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 8x (2-point, 3rd-order) (z-form)
float optimal8_23z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float c0 = even1 * 0.50004007194083089;
	float c1 = odd1 * 1.06397659072500650;
	float c2 = even1 * -0.000480863289971321;
	float c3 = odd1 * -0.73514591836770027;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 16x (2-point, 3rd-order) (z-form)
float optimal16_23z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float c0 = even1 * 0.50001096675880796;
	float c1 = odd1 * 1.03585606328743830;
	float c2 = even1 * -0.000131601105693441;
	float c3 = odd1 * -0.38606621963374965;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 32x (2-point, 3rd-order) (z-form)
float optimal_32z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float c0 = even1 * 0.50000286037713559;
	float c1 = odd1 * 1.01889120864375270;
	float c2 = even1 * -0.000034324525627571;
	float c3 = odd1 * -0.19775766248673177;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 2x (4-point, 2nd-order) (z-form)
float optimal2_42z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.42334633257225274 + even2 * 0.07668732202139628;
	float c1 = odd1 * 0.26126047291143606 + odd2 * 0.24778879018226652;
	float c2 = even1 * -0.213439787561776841 + even2 * 0.21303593243799016;
	return (c2 * z + c1) * z + c0;
}

// Optimal 4x (4-point, 2nd-order) (z-form)
float optimal4_42z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.38676264891201206 + even2 * 0.11324319172521946;
	float c1 = odd1 * 0.01720901456660906 + odd2 * 0.32839294317251788;
	float c2 = even1 * -0.228653995318581881 + even2 * 0.22858390767180370;
	return (c2 * z + c1) * z + c0;
}

// Optimal 8x (4-point, 2nd-order) (z-form)
float optimal8_42z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.32852206663814043 + even2 * 0.17147870380790242;
	float c1 = odd1 * -0.35252373075274990 + odd2 * 0.45113687946292658;
	float c2 = even1 * -0.240052062078895181 + even2 * 0.24004281672637814;
	return (c2 * z + c1) * z + c0;
}

// Optimal 16x (4-point, 2nd-order) (z-form)
float optimal16_42z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.20204741371575463 + even2 * 0.29795268253813623;
	float c1 = odd1 * -1.11855475338366150 + odd2 * 0.70626377291054832;
	float c2 = even1 * -0.245061178654743641 + even2 * 0.24506002360805534;
	return (c2 * z + c1) * z + c0;
}

// Optimal 32x (4-point, 2nd-order) (z-form)
float optimal32_42z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * -0.04817865217726123 + even2 * 0.54817866412548932;
	float c1 = odd1 * -2.62328241292796620 + odd2 * 1.20778105913587620;
	float c2 = even1 * -0.247552438397138281 + even2 * 0.24755229501840223;
	return (c2 * z + c1) * z + c0;
}

// Optimal 2x (4-point, 3rd-order) (z-form)
float optimal2_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.45868970870461956 + even2 * 0.04131401926395584;
	float c1 = odd1 * 0.48068024766578432 + odd2 * 0.17577925564495955;
	float c2 = even1 * -0.246185007019907091 + even2 * 0.24614027139700284;
	float c3 = odd1 * -0.36030925263849456 + odd2 * 0.10174985775982505;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 4x (4-point, 3rd-order) (z-form)
float optimal4_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46209345013918979 + even2 * 0.03790693583186333;
	float c1 = odd1 * 0.51344507801315964 + odd2 * 0.16261507145522014;
	float c2 = even1 * -0.248540332990294211 + even2 * 0.24853570133765701;
	float c3 = odd1 * -0.42912649274763925 + odd2 * 0.13963062613760227;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 8x (4-point, 3rd-order) (z-form)
float optimal8_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46360002085841184 + even2 * 0.03640000638072349;
	float c1 = odd1 * 0.52776949859997280 + odd2 * 0.15746108253367153;
	float c2 = even1 * -0.249658121535793251 + even2 * 0.24965779466617388;
	float c3 = odd1 * -0.46789242171187317 + odd2 * 0.15551896027602030;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 16x (4-point, 3rd-order) (z-form)
float optimal16_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46436507349411416 + even2 * 0.03563492826010761;
	float c1 = odd1 * 0.53463126553787166 + odd2 * 0.15512856361039451;
	float c2 = even1 * -0.249923540967159741 + even2 * 0.24992351991649797;
	float c3 = odd1 * -0.48601256046234864 + odd2 * 0.16195131297091253;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 32x (4-point, 3rd-order) (z-form)
float optimal32_43z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46465589031535864 + even2 * 0.03534410979496938;
	float c1 = odd1 * 0.53726845877054186 + odd2 * 0.15424449410914165;
	float c2 = even1 * -0.249981930954029101 + even2 * 0.24998192963009191;
	float c3 = odd1 * -0.49369595780454456 + odd2 * 0.16455902278580614;
	return ((c3 * z + c2) * z + c1) * z + c0;
}

// Optimal 2x (4-point, 4th-order) (z-form)
float optimal2_44z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.45645918406487612 + even2 * 0.04354173901996461;
	float c1 = odd1 * 0.47236675362442071 + odd2 * 0.17686613581136501;
	float c2 = even1 * -0.253674794204558521 + even2 * 0.25371918651882464;
	float c3 = odd1 * -0.37917091811631082 + odd2 * 0.11952965967158000;
	float c4 = even1 * 0.04252164479749607 + even2 * -0.04289144034653719;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 4x (4-point, 4th-order) (z-form)
float optimal4_44z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46567255120778489 + even2 * 0.03432729708429672;
	float c1 = odd1 * 0.53743830753560162 + odd2 * 0.15429462557307461;
	float c2 = even1 * -0.251942101340217441 + even2 * 0.25194744935939062;
	float c3 = odd1 * -0.46896069955075126 + odd2 * 0.15578800670302476;
	float c4 = even1 * 0.00986988334359864 + even2 * -0.00989340017126506;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 8x (4-point, 4th-order) (z-form)
float optimal8_44z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46771532012068961 + even2 * 0.03228466824404497;
	float c1 = odd1 * 0.55448654344364423 + odd2 * 0.14851181120641987;
	float c2 = even1 * -0.250587283698110121 + even2 * 0.25058765188457821;
	float c3 = odd1 * -0.49209020939096676 + odd2 * 0.16399414834151946;
	float c4 = even1 * 0.00255074537015887 + even2 * -0.00255226912537286;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 16x (4-point, 4th-order) (z-form)
float optimal16_44z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46822774170144532 + even2 * 0.03177225758005808;
	float c1 = odd1 * 0.55890365706150436 + odd2 * 0.14703258836343669;
	float c2 = even1 * -0.250153411893796031 + even2 * 0.25015343462990891;
	float c3 = odd1 * -0.49800710906733769 + odd2 * 0.16600005174304033;
	float c4 = even1 * 0.00064264050033187 + even2 * -0.00064273459469381;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 32x (4-point, 4th-order) (z-form)
float optimal32_44z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float c0 = even1 * 0.46835497211269561 + even2 * 0.03164502784253309;
	float c1 = odd1 * 0.56001293337091440 + odd2 * 0.14666238593949288;
	float c2 = even1 * -0.250038759826233691 + even2 * 0.25003876124297131;
	float c3 = odd1 * -0.49949850957839148 + odd2 * 0.16649935475113800;
	float c4 = even1 * 0.00016095224137360 + even2 * -0.00016095810460478;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 2x (6-point, 4th-order) (z-form)
float optimal2_64z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.37484203669443822 + even2 * 0.11970939637439368 + even3 * 0.00544862268096358;
	float c1 = odd1 * 0.19253897284651597 + odd2 * 0.22555179040018719 + odd3 * 0.02621377625620669;
	float c2 = even1 * -0.154026006475653071 + even2 * 0.10546111301131367 + even3 * 0.04856757454258609;
	float c3 = odd1 * -0.06523685579716083 + odd2 * -0.04867197815057284 + odd3 * 0.04200764942718964;
	float c4 = even1 * 0.03134095684084392 + even2 * -0.04385804833432710 + even3 * 0.01249475765486819;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 4x (6-point, 4th-order) (z-form)
float optimal4_64z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.26148143200222657 + even2 * 0.22484494681472966 + even3 * 0.01367360612950508;
	float c1 = odd1 * -0.20245593827436142 + odd2 * 0.29354348112881601 + odd3 * 0.06436924057941607;
	float c2 = even1 * -0.022982104451679701 + even2 * -0.09068617668887535 + even3 * 0.11366875749521399;
	float c3 = odd1 * 0.36296419678970931 + odd2 * -0.26421064520663945 + odd3 * 0.08591542869416055;
	float c4 = even1 * 0.02881527997393852 + even2 * -0.04250898918476453 + even3 * 0.01369173779618459;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 8x (6-point, 4th-order) (z-form)
float optimal8_64z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.07571827673995030 + even2 * 0.39809419102537769 + even3 * 0.02618753167558019;
	float c1 = odd1 * -0.87079480370960549 + odd2 * 0.41706012247048818 + odd3 * 0.12392296259397995;
	float c2 = even1 * 0.186883718356452901 + even2 * -0.40535151498252686 + even3 * 0.21846781431808182;
	float c3 = odd1 * 1.09174419992174300 + odd2 * -0.62917625718809478 + odd3 * 0.15915674384870970;
	float c4 = even1 * 0.03401038103941584 + even2 * -0.05090907029392906 + even3 * 0.01689861603514873;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 16x (6-point, 4th-order) (z-form)
float optimal16_64z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * -0.30943127416213301 + even2 * 0.75611844407537543 + even3 * 0.05331283006820442;
	float c1 = odd1 * -2.23586327978235700 + odd2 * 0.66020840412562265 + odd3 * 0.25104761112921636;
	float c2 = even1 * 0.625420761014402691 + even2 * -1.06313460380183860 + even3 * 0.43771384337431529;
	float c3 = odd1 * 2.57088518304678090 + odd2 * -1.36878543609177150 + odd3 * 0.30709424868485174;
	float c4 = even1 * 0.03755086455339280 + even2 * -0.05631219122315393 + even3 * 0.01876132424143207;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 32x (6-point, 4th-order) (z-form)
float optimal32_64z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * -1.05730227922290790 + even2 * 1.45069541587021430 + even3 * 0.10660686335233649;
	float c1 = odd1 * -4.87455554035028720 + odd2 * 1.12509567592532630 + odd3 * 0.49985370215839708;
	float c2 = even1 * 1.479370435823112101 + even2 * -2.34405608915933780 + even3 * 0.86468565335070746;
	float c3 = odd1 * 5.42677291742286180 + odd2 * -2.79672428287565160 + odd3 * 0.59267998874843331;
	float c4 = even1 * 0.03957507923965987 + even2 * -0.05936083498715066 + even3 * 0.01978575568000696;
	return (((c4 * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 2x (6-point, 5th-order) (z-form)
float optimal2_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.40513396007145713 + even2 * 0.09251794438424393 + even3 * 0.00234806603570670;
	float c1 = odd1 * 0.28342806338906690 + odd2 * 0.21703277024054901 + odd3 * 0.01309294748731515;
	float c2 = even1 * -0.191337682540351941 + even2 * 0.16187844487943592 + even3 * 0.02946017143111912;
	float c3 = odd1 * -0.16471626190554542 + odd2 * -0.00154547203542499 + odd3 * 0.03399271444851909;
	float c4 = even1 * 0.03845798729588149 + even2 * -0.05712936104242644 + even3 * 0.01866750929921070;
	float c5 = odd1 * 0.04317950185225609 + odd2 * -0.01802814255926417 + odd3 * 0.00152170021558204;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 4x (6-point, 5th-order) (z-form)
float optimal4_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.41496902959240894 + even2 * 0.08343081932889224 + even3 * 0.00160015038681571;
	float c1 = odd1 * 0.31625515004859783 + odd2 * 0.21197848565176958 + odd3 * 0.00956166668408054;
	float c2 = even1 * -0.203271896548875371 + even2 * 0.17989908432249280 + even3 * 0.02337283412161328;
	float c3 = odd1 * -0.20209241069835732 + odd2 * 0.01760734419526000 + odd3 * 0.02985927012435252;
	float c4 = even1 * 0.04100948858761910 + even2 * -0.06147760875085254 + even3 * 0.02046802954581191;
	float c5 = odd1 * 0.06607747864416924 + odd2 * -0.03255079211953620 + odd3 * 0.00628989632244913;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 8x (6-point, 5th-order) (z-form)
float optimal8_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.41660797292569773 + even2 * 0.08188468587188069 + even3 * 0.00150734119050266;
	float c1 = odd1 * 0.32232780822726981 + odd2 * 0.21076321997422021 + odd3 * 0.00907649978070957;
	float c2 = even1 * -0.205219993961471501 + even2 * 0.18282942057327367 + even3 * 0.02239057377093268;
	float c3 = odd1 * -0.21022298520246224 + odd2 * 0.02176417471349534 + odd3 * 0.02898626924395209;
	float c4 = even1 * 0.04149963966704384 + even2 * -0.06224707096203808 + even3 * 0.02074742969707599;
	float c5 = odd1 * 0.07517133281176167 + odd2 * -0.03751837438141215 + odd3 * 0.00747588873055296;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 16x (6-point, 5th-order) (z-form)
float optimal16_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.41809989254549901 + even2 * 0.08049339946273310 + even3 * 0.00140670799165932;
	float c1 = odd1 * 0.32767596257424964 + odd2 * 0.20978189376640677 + odd3 * 0.00859567104974701;
	float c2 = even1 * -0.206944618112960001 + even2 * 0.18541689550861262 + even3 * 0.02152772260740132;
	float c3 = odd1 * -0.21686095413034051 + odd2 * 0.02509557922091643 + odd3 * 0.02831484751363800;
	float c4 = even1 * 0.04163046817137675 + even2 * -0.06244556931623735 + even3 * 0.02081510113314315;
	float c5 = odd1 * 0.07990500783668089 + odd2 * -0.03994519162531633 + odd3 * 0.00798609327859495;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}

// Optimal 32x (6-point, 5th-order) (z-form)
float optimal32_65z(const float *y, const float x)
{
	float z = x - 1 / 2.0;
	float even1 = y[1] + y[0], odd1 = y[1] - y[0];
	float even2 = y[2] + y[-1], odd2 = y[2] - y[-1];
	float even3 = y[3] + y[-2], odd3 = y[3] - y[-2];
	float c0 = even1 * 0.42685983409379380 + even2 * 0.07238123511170030 + even3 * 0.00075893079450573;
	float c1 = odd1 * 0.35831772348893259 + odd2 * 0.20451644554758297 + odd3 * 0.00562658797241955;
	float c2 = even1 * -0.217009177221292431 + even2 * 0.20051376594086157 + even3 * 0.01649541128040211;
	float c3 = odd1 * -0.25112715343740988 + odd2 * 0.04223025992200458 + odd3 * 0.02488727472995134;
	float c4 = even1 * 0.04166946673533273 + even2 * -0.06250420114356986 + even3 * 0.02083473440841799;
	float c5 = odd1 * 0.08349799235675044 + odd2 * -0.04174912841630993 + odd3 * 0.00834987866042734;
	return ((((c5 * z + c4) * z + c3) * z + c2) * z + c1) * z + c0;
}


//--------------------------------------------------------------------------------------------------------------
int srandfq = 74235;
void srandf(const int seed)
{
    srandfq = seed;
}
float randf()
{
    // https://www.musicdsp.org/en/latest/Other/273-fast-float-random-numbers.html
    // moc.liamg@seir.kinimod
    srandfq *= 16807;
    return (float)(srandfq & 0x7FFFFFFF) * 4.6566129e-010f;
}

// microtimer
uint64_t microtime()
{
	struct timeval tv;
	struct timezone tz;
	memset(&tz, 0, sizeof(struct timezone));
	gettimeofday(&tv, &tz);
	return 1000000 * tv.tv_sec + tv.tv_usec;
}

// sample buffer
float data[8];
float rndData()
{
	for(int i = 0; i < 8; i++)
		data[i] = randf();
	return data[0];
}

// main
int main()
{
	// set highest priority
	errno = 0;
	if(nice(-20) < 0)
	{
		while(errno != 0)
		{
			errno = 0;
			if(nice(-20) < 0)
				printf("Attempting to set process to nice of -20 (run with sudo)...\n");
			sleep(1);
		}
	}

	// init random seed
	srandf(time(0));

	// prep benchmarks
	setlocale(LC_NUMERIC, "");
	float ret = 0;
	unsigned long e = 0, ui = 0;
	uint64_t st = 0, et = 0, avg = 0;
	printf("\nExecutions in %'u seconds (higher is better) ...\n\n", interval / 1000000);

	///-------------------------------

	// rndData()
	e = 0; st = microtime();
	while(microtime() - st <= interval){ret += rndData(); e++;}
	printf("rndData(): %'lu\n", e);

	///-------------------------------

	// order3Spline()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += order3Spline(data[0], data[1], data[2], data[3], data[4]); e++;}
	printf("order3Spline(): %'lu\n", e);

	///-------------------------------

	// hermite1()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite1(data[0], data[1], data[2], data[3], data[4]); e++;}
	printf("hermite1(): %'lu\n", e);

	// hermite2()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite2(data[0], data[1], data[2], data[3], data[4]); e++;}
	printf("hermite2(): %'lu\n", e);

	// hermite3()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite3(data[0], data[1], data[2], data[3], data[4]); e++;}
	printf("hermite3(): %'lu\n", e);

	// hermite4()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite4(data[0], data[1], data[2], data[3], data[4]); e++;}
	printf("hermite4(): %'lu\n", e);

	///-------------------------------

	// LinearInterpolate()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += LinearInterpolate(data[0], data[1], data[2]); e++;}
	printf("LinearInterpolate(): %'lu\n", e);

	// CosineInterpolate()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += CosineInterpolate(data[0], data[1], data[2]); e++;}
	printf("CosineInterpolate(): %'lu\n", e);

	/// CubicInterpolate()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += CubicInterpolate(data[0], data[1], data[2], data[3], data[4]); e++;}
	printf("CubicInterpolate(): %'lu\n", e);

	// CatmullInterpolate()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += CatmullInterpolate(data[0], data[1], data[2], data[3], data[4]); e++;}
	printf("CatmullInterpolate(): %'lu\n", e);

	// HermiteInterpolate()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += HermiteInterpolate(data[0], data[1], data[2], data[3], data[4], data[5], data[6]); e++;}
	printf("HermiteInterpolate(): %'lu\n", e);

	///-------------------------------

	// bspline_43x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += bspline_43x(&data[3], data[0]); e++;}
	printf("bspline_43x(): %'lu\n", e);

	// bspline_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += bspline_43z(&data[3], data[0]); e++;}
	printf("bspline_43z(): %'lu\n", e);

	// bspline_65x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += bspline_65x(&data[3], data[0]); e++;}
	printf("bspline_65x(): %'lu\n", e);

	// lagrange_43x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += lagrange_43x(&data[3], data[0]); e++;}
	printf("lagrange_43x(): %'lu\n", e);

	// lagrange_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += lagrange_43z(&data[3], data[0]); e++;}
	printf("lagrange_43z(): %'lu\n", e);

	// lagrange_65x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += lagrange_65x(&data[3], data[0]); e++;}
	printf("lagrange_65x(): %'lu\n", e);

	// lagrange_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += lagrange_65z(&data[3], data[0]); e++;}
	printf("lagrange_65z(): %'lu\n", e);

	// hermite_43x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite_43x(&data[3], data[0]); e++;}
	printf("hermite_43x(): %'lu\n", e);

	// hermite_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite_43z(&data[3], data[0]); e++;}
	printf("hermite_43z(): %'lu\n", e);

	// hermite_63x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite_63x(&data[3], data[0]); e++;}
	printf("hermite_63x(): %'lu\n", e);

	// hermite_63z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite_63z(&data[3], data[0]); e++;}
	printf("hermite_63z(): %'lu\n", e);

	// hermite_65x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite_65x(&data[3], data[0]); e++;}
	printf("hermite_65x(): %'lu\n", e);

	// hermite_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += hermite_65z(&data[3], data[0]); e++;}
	printf("hermite_65z(): %'lu\n", e);

	// osculating_45x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += osculating_45x(&data[3], data[0]); e++;}
	printf("osculating_45x(): %'lu\n", e);

	// osculating_45z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += osculating_45z(&data[3], data[0]); e++;}
	printf("osculating_45z(): %'lu\n", e);

	// osculating_65x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += osculating_65x(&data[3], data[0]); e++;}
	printf("osculating_65x(): %'lu\n", e);

	// osculating_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += osculating_65z(&data[3], data[0]); e++;}
	printf("osculating_65z(): %'lu\n", e);

	// watte_42x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += watte_42x(&data[3], data[0]); e++;}
	printf("watte_42x(): %'lu\n", e);

	// watte_42z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += watte_42z(&data[3], data[0]); e++;}
	printf("watte_42z(): %'lu\n", e);

	// parabolic_42x()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += parabolic_42x(&data[3], data[0]); e++;}
	printf("parabolic_42x(): %'lu\n", e);

	// parabolic_42z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += parabolic_42z(&data[3], data[0]); e++;}
	printf("parabolic_42z(): %'lu\n", e);

	// optimal2_23z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal2_23z(&data[3], data[0]); e++;}
	printf("optimal2_23z(): %'lu\n", e);

	// optimal4_23z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal4_23z(&data[3], data[0]); e++;}
	printf("optimal4_23z(): %'lu\n", e);

	// optimal8_23z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal8_23z(&data[3], data[0]); e++;}
	printf("optimal8_23z(): %'lu\n", e);

	// optimal16_23z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal16_23z(&data[3], data[0]); e++;}
	printf("optimal16_23z(): %'lu\n", e);

	// optimal_32z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal_32z(&data[3], data[0]); e++;}
	printf("optimal_32z(): %'lu\n", e);

	// optimal2_42z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal2_42z(&data[3], data[0]); e++;}
	printf("optimal2_42z(): %'lu\n", e);

	// optimal4_42z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal4_42z(&data[3], data[0]); e++;}
	printf("optimal4_42z(): %'lu\n", e);

	// optimal8_42z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal8_42z(&data[3], data[0]); e++;}
	printf("optimal8_42z(): %'lu\n", e);

	// optimal16_42z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal16_42z(&data[3], data[0]); e++;}
	printf("optimal16_42z(): %'lu\n", e);

	// optimal32_42z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal32_42z(&data[3], data[0]); e++;}
	printf("optimal32_42z(): %'lu\n", e);

	// optimal2_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal2_43z(&data[3], data[0]); e++;}
	printf("optimal2_43z(): %'lu\n", e);

	// optimal4_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal4_43z(&data[3], data[0]); e++;}
	printf("optimal4_43z(): %'lu\n", e);

	// optimal8_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal8_43z(&data[3], data[0]); e++;}
	printf("optimal8_43z(): %'lu\n", e);

	// optimal16_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal16_43z(&data[3], data[0]); e++;}
	printf("optimal16_43z(): %'lu\n", e);

	// optimal32_43z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal32_43z(&data[3], data[0]); e++;}
	printf("optimal32_43z(): %'lu\n", e);

	// optimal2_44z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal2_44z(&data[3], data[0]); e++;}
	printf("optimal2_44z(): %'lu\n", e);

	// optimal4_44z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal4_44z(&data[3], data[0]); e++;}
	printf("optimal4_44z(): %'lu\n", e);

	// optimal8_44z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal8_44z(&data[3], data[0]); e++;}
	printf("optimal8_44z(): %'lu\n", e);

	// optimal16_44z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal16_44z(&data[3], data[0]); e++;}
	printf("optimal16_44z(): %'lu\n", e);

	// optimal32_44z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal32_44z(&data[3], data[0]); e++;}
	printf("optimal32_44z(): %'lu\n", e);

	// optimal2_64z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal2_64z(&data[3], data[0]); e++;}
	printf("optimal2_64z(): %'lu\n", e);

	// optimal4_64z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal4_64z(&data[3], data[0]); e++;}
	printf("optimal4_64z(): %'lu\n", e);

	// optimal8_64z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal8_64z(&data[3], data[0]); e++;}
	printf("optimal8_64z(): %'lu\n", e);

	// optimal16_64z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal16_64z(&data[3], data[0]); e++;}
	printf("optimal16_64z(): %'lu\n", e);

	// optimal32_64z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal32_64z(&data[3], data[0]); e++;}
	printf("optimal32_64z(): %'lu\n", e);

	// optimal2_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal2_65z(&data[3], data[0]); e++;}
	printf("optimal2_65z(): %'lu\n", e);

	// optimal4_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal4_65z(&data[3], data[0]); e++;}
	printf("optimal4_65z(): %'lu\n", e);

	// optimal8_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal8_65z(&data[3], data[0]); e++;}
	printf("optimal8_65z(): %'lu\n", e);

	// optimal16_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal16_65z(&data[3], data[0]); e++;}
	printf("optimal16_65z(): %'lu\n", e);

	// optimal32_65z()
	e = 0; st = microtime();
	while(microtime() - st <= interval){rndData(); ret += optimal32_65z(&data[3], data[0]); e++;}
	printf("optimal32_65z(): %'lu\n", e);

	// done
	printf("%c\n", (char)ret);  // forces the compiler to not disregard the functions we are testing
	return 0;
}