The published documents on how to generate sse instructions via gcc's vector
arithmetic support are generally poor and contain conflicting information. I've
been experimenting with different ways to generate sse, sse2, and mmx code
on the x86_64. The x86_64 is notable in that it has 16 sse registers (rather
than 8) and that by default it does all math via those registers by default. It
seems to do a halfway decent job of scheduling non-vectorized code but I wanted
to find a way to exercise more of the instruction set.

The case I decided to play with first was cubic interpolation in Linuxsampler. I
am NOT trying to optimize for a specific case for any particular reason, just
trying to get to where I understand SSE well enough to apply it when there *is*
a reason.  Rather than let the knowledge rot, I thought I'd stick it in here in
the hope it might be useful to others, especially those that are hand optimizing
something where getting the compiler to do some of the heavy lifting is helpful. 


Here's some C code that at least compiles - and may not otherwise be correct -
and only generates 37 sse instructions (vs the original stereo code which
generated over 60). My note on this is that a stereo version would not be much
more complex (basically just ignore the upper half of the sse registers) and
still result in less instructions than the default stereo cubic interpolation
code. It would probably still be outperformed by the original stereo code,
however, and perhaps the 3dnow code would be the only thing that could
outperform that.

Note the v4sf declaration in the union is second. Declared that way, you can
still do structure assignment of const variables, and just specify the vector
part of the union during the math.

Also, constants have to be declared as per below in order to retain the
expressiveness of vector = vector + vector.

Telling the compiler it's a v2sf (for stereo) is hell on sse code generation,
but does nice stuff when 3dnow is specified....

Compilation line (on gcc4): g++ -O3 -msse -msse2  -fverbose-asm -S t.cpp

#include <iostream>
#include <stdio.h>

typedef struct {
        float left;
        float right;
        float front;
        float rear;
} quad_samples;

typedef float v4sf __attribute__ ((vector_size (16)));

typedef union {
        quad_samples s;
        v4sf v;
} quad_sample_t;

const quad_sample_t c3 = {3.0f,3.0f,3.0f,3.0f};
const quad_sample_t cpoint5 = {.5f,.5f,.5f,.5f};
const quad_sample_t c2 = { 2.0f,2.0f,2.0f,2.0f};
const quad_sample_t c5 = { 5.0f,5.0f,5.0f,5.0f};
static quad_sample_t pos_fract;

typedef float sample_t;

static quad_sample_t Interpolate1StepQuadCPP(quad_sample_t* pSrc, double* Pos,
float& Pitch) {
     int   pos_int   = (long) *Pos;  // integer position
//     float pos_fract = *Pos - pos_int;     // fractional part of position
     pos_int <<= 1;
     quad_sample_t samplePoint;
                    quad_sample_t xm1 = pSrc[pos_int];
                    quad_sample_t x0  = pSrc[pos_int+1];
                    quad_sample_t x1  = pSrc[pos_int+2];
                    quad_sample_t x2  = pSrc[pos_int+3];
                    quad_sample_t a;
                    a.v  = (c3.v * (x0.v - x1.v) - xm1.v + x2.v) * cpoint5.v;
                    quad_sample_t b;
                    b.v  = c2.v * x1.v + xm1.v- (c5.v* x0.v + x2.v) * cpoint5.v;
                    quad_sample_t c;
                    c.v  = (x1.v - xm1.v) * cpoint5.v;
                    samplePoint.v = (((a.v * pos_fract.v) + b.v) * pos_fract.v +
c.v) * pos_fract.v + x0.v;

                *Pos += Pitch;
                return samplePoint;
            }

static double t = 1.0;
static float t1 = .5;
extern void b (quad_sample_t a) ;

main(char *argv[], int argc) {
        quad_sample_t result;
        quad_sample_t samples[1024];
        result = Interpolate1StepQuadCPP(samples,&t,t1);
        b (result);
//quad_sample_t results = (quad_sample_t) result;
//      printf("%f:,%f:",results.t.left,results.t.right);
}