What follows are a pair of code example for finding an series of objects' speeds along the vector between them and the origin given the objects' 3D velocities and positions. It is written to use the AltiVec unit. This example does not handle arrays that are not 16 byte aligned and object counts that are not an even multiple of eight.

void CalcSpeed( float *velocities, float *distances, float *results, int count )
{

register float vX1, vY1, vZ1; //velocity 1
register float vX2, vY2, vZ2; //velocity 2
register float vX3, vY3, vZ3; //velocity 3
register float dX1, dY1, dZ1; //distance 1
register float dX2, dY2, dZ2; //distance 2
register float dX3, dY3, dZ3; //distance 3
register float DD1, DD2, DD3; //D dot D, 1 through 3
register float DV1, DV2, DV3; //D dot V, 1 through 3
register float rsqrt1, rsqrt2, rsqrt3;
register float one = 1.0;
register float oneHalf = 0.5;
register float *vel = (float*) velocities;
register float *pos = (float*) distances;

//Process three sounds per loop
for( ; count >= 3; count -= 3 )
{

//Load X, Y and Z velocity components for three
//input distances from the listener
dX1 = pos[0];
dY1 = pos[1];
dZ1 = pos[2];
dX2 = pos[3];
dY2 = pos[4];
dZ2 = pos[5];
dX3 = pos[6];
dY3 = pos[7];
dZ3 = pos[8];
pos += 9;

//Load X, Y and Z velocity components for
//three input velocities
vX1 = vel[0];
vY1 = vel[1];
vZ1 = vel[2];
vX2 = vel[3];
vY2 = vel[4];
vZ2 = vel[5];
vX3 = vel[6];
vY3 = vel[7];
vZ3 = vel[8];
vel += 9;

//Calculate D dot V, and store in DV#
DV1 = vX1 * dX1 + vY1 * dY1 + vZ1 * dZ1;
DV2 = vX2 * dX2 + vY2 * dY2 + vZ2 * dZ2;
DV3 = vX3 * dX3 + vY3 * dY3 + vZ3 * dZ3;

//Calculate D dot D for each of the three vectors
DD1 = dX1 * dX1 + dY1 * dY1 + dZ1 * dZ1;
DD2 = dX2 * dX2 + dY2 * dY2 + dZ2 * dZ2;
DD3 = dX3 * dX3 + dY3 * dY3 + dZ3 * dZ3;

//Calculate (D dot D)-0.5
rsqrt1 = __frsqrte( DD1 );
rsqrt2 = __frsqrte( DD2 );
rsqrt3 = __frsqrte( DD3 );

//Do one round of Newton-Raphson refinement
//because the results from frsqrte are
//only accurate to one part in 32. This should
//give us about one part in 1024 accuracy, which
//is enough for audio volume, for which we only need one
//part in 256 for the Sound Manager
rsqrt1 += oneHalf * rsqrt1 * ( one - DD1 * rsqrt1 * rsqrt1 );
rsqrt2 += oneHalf * rsqrt2 * ( one - DD2 * rsqrt2 * rsqrt2 );
rsqrt3 += oneHalf * rsqrt3 * ( one - DD3 * rsqrt3 * rsqrt3 );

//store out the result
//result = (D dot V) * 1.0 / sqrt( D dot D)
(results++)[0] = DV1 * rsqrt1;
(results++)[0] = DV2 * rsqrt2;
(results++)[0] = DV3 * rsqrt3;

}

//Process the remaining sounds
//This is the same algorithm as is used above
for( ; count > 0; count-- )
{

//Load X, Y and Z velocity components for three
//input distances from the listener
dX1 = pos[0];
dY1 = pos[1];
dZ1 = pos[2];
pos += 3;

//Load X, Y and Z velocity components for
//three input velocities
vX1 = vel[0];
vY1 = vel[1];
vZ1 = vel[2];
vel += 3;

//Calculate D dot V, and store in DV#
DV1 = vX1 * dX1 + vY1 * dY1 + vZ1 * dZ1;

//Calculate D dot D for each of the three vectors
DD1 = dX1 * dX1 + dY1 * dY1 + dZ1 * dZ1;

//Calculate (D dot D)-0.5
rsqrt1 = __frsqrte( DD1 );

//Do one round of Newton-Raphson refinement
//because the results from frsqrte are
//only accurate to one part in 32. This should
//give us about one part in 1024 accuracy, which
//is enough for audio volume, for which we only need one
//part in 256 for the Sound Manager
rsqrt1 += oneHalf * rsqrt1 * ( one - DD1 * rsqrt1 * rsqrt1 );

//store out the result
//result = (D dot V) * 1.0 / sqrt( D dot D)
(results++)[0] = DV1 * rsqrt1;

}

}

//velocities and distances are input as {x, y, z} interleaved 3-tuples
//This simplified function requires that velocities, distances
//and results are all 16 byte aligned. In addition, count must
//be a even multiple of 8.
void CalcSpeedVec( vector float *velocities, vector float *distances, vector float *results, int count )
{

int loopCount = count / ( 2 * vec_step( vector float ));
vector float zero = (vector float) vec_splat_u32(0);
vector unsigned char mergeLowHigh = (vector unsigned char)
( 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15, 20, 21, 22, 23 );
vector unsigned char mergeHighLow = vec_xor( mergeLowHigh, vec_splat_u8(8));

while( loopCount-- )
{

vector float vLoad1, vLoad2, vLoad3, vLoad4, vLoad5, vLoad6;
vector float dLoad1, dLoad2, dLoad3, dLoad4, dLoad5, dLoad6;
vector float dX1, dY1, dZ1, dX2, dY2, dZ2; //x, y, and z distances
vector float vX1, vY1, vZ1, vX2, vY2, vZ2; //x, y, and z velocities
vector float D_dot_D1, D_dot_V1, D_dot_D2, D_dot_V2;

//Load eight vector[3]'s from the distance and velocity arrays
dLoad1 = vec_ld( 0 * sizeof( vector float ), distances );
dLoad2 = vec_ld( 1 * sizeof( vector float ), distances );
dLoad3 = vec_ld( 2 * sizeof( vector float ), distances );
dLoad4 = vec_ld( 3 * sizeof( vector float ), distances );
dLoad5 = vec_ld( 4 * sizeof( vector float ), distances );
dLoad6 = vec_ld( 5 * sizeof( vector float ), distances );
distances += 6;

vLoad1 = vec_ld( 0 * sizeof( vector float ), velocities );
vLoad2 = vec_ld( 1 * sizeof( vector float ), velocities );
vLoad3 = vec_ld( 2 * sizeof( vector float ), velocities );
vLoad4 = vec_ld( 3 * sizeof( vector float ), velocities );
vLoad5 = vec_ld( 4 * sizeof( vector float ), velocities );
vLoad6 = vec_ld( 5 * sizeof( vector float ), velocities );
velocities += 6;

//Now, we transpose the format using a 4x3 matrix transpose
//to arrive at uniform vectors containing the four x's,
//four y's or four z's
//
//If you store your data in a planar format, rather
//than interleaved like this then you dont have to
//do this part.
//
// { X1, Y1, Z1, X2 } { X1, X2, X3, X4 }
// { Y2, Z2, X3, Y3 } -----> { Y1, Y2, Y3, Y4 }
// { Z3, X4, Y4, Z4 } { Z1, Z2, Z3, Z4 }
//
//It is possible to do each transpose in just five permutes
//rather than six.
//That would save a few cycles. Can you figure out how?
{

vector float temp1, temp2, temp3;

temp1 = vec_perm( dLoad1, dLoad2, mergeHighLow ); // X1, X3, Y1, Y3
temp2 = vec_perm( dLoad1, dLoad3, mergeLowHigh ); // Z1, Z3, X2, X4
temp3 = vec_perm( dLoad2, dLoad3, mergeHighLow ); // Y2, Y4, Z2, Z4
dX1 = vec_perm( temp1, temp2, mergeHighLow ); // X1, X2, X3, X4
dY1 = vec_perm( temp1, temp3, mergeLowHigh ); // Y1, Y2, Y3, Y4
dZ1 = vec_perm( temp2, temp3, mergeHighLow ); // Z1, Z2, Z3, Z4

temp1 = vec_perm( dLoad4, dLoad5, mergeHighLow ); // X1, X3, Y1, Y3
temp2 = vec_perm( dLoad4, dLoad6, mergeLowHigh ); // Z1, Z3, X2, X4
temp3 = vec_perm( dLoad5, dLoad6, mergeHighLow ); // Y2, Y4, Z2, Z4
dX2 = vec_perm( temp1, temp2, mergeHighLow ); // X1, X2, X3, X4
dY2 = vec_perm( temp1, temp3, mergeLowHigh ); // Y1, Y2, Y3, Y4
dZ2 = vec_perm( temp2, temp3, mergeHighLow ); // Z1, Z2, Z3, Z4

temp1 = vec_perm( vLoad1, vLoad2, mergeHighLow ); // X1, X3, Y1, Y3
temp2 = vec_perm( vLoad1, vLoad3, mergeLowHigh ); // Z1, Z3, X2, X4
temp3 = vec_perm( vLoad2, vLoad3, mergeHighLow ); // Y2, Y4, Z2, Z4
vX1 = vec_perm( temp1, temp2, mergeHighLow ); // X1, X2, X3, X4
vY1 = vec_perm( temp1, temp3, mergeLowHigh ); // Y1, Y2, Y3, Y4
vZ1 = vec_perm( temp2, temp3, mergeHighLow ); // Z1, Z2, Z3, Z4

temp1 = vec_perm( vLoad4, vLoad5, mergeHighLow ); // X1, X3, Y1, Y3
temp2 = vec_perm( vLoad4, vLoad6, mergeLowHigh ); // Z1, Z3, X2, X4
temp3 = vec_perm( vLoad5, vLoad6, mergeHighLow ); // Y2, Y4, Z2, Z4
vX2 = vec_perm( temp1, temp2, mergeHighLow ); // X1, X2, X3, X4
vY2 = vec_perm( temp1, temp3, mergeLowHigh ); // Y1, Y2, Y3, Y4
vZ2 = vec_perm( temp2, temp3, mergeHighLow ); // Z1, Z2, Z3, Z4

}

//Do the dot products: D dot D and D dot V
//Some of these will happen concurrently with the transpose operations above
//X part of the dot product
D_dot_D1 = vec_madd( dX1, dX1, zero );
D_dot_D2 = vec_madd( dX2, dX2, zero );
D_dot_V1 = vec_madd( dX1, vX1, zero );
D_dot_V2 = vec_madd( dX2, vX2, zero );

//Y part of the dot product
D_dot_D1 = vec_madd( dY1, dY1, D_dot_D1 );
D_dot_D2 = vec_madd( dY2, dY2, D_dot_D2 );
D_dot_V1 = vec_madd( dY1, vY1, D_dot_V1 );
D_dot_V2 = vec_madd( dY2, vY2, D_dot_V2 );

//Z part of the dot product
D_dot_D1 = vec_madd( dZ1, dZ1, D_dot_D1 );
D_dot_D2 = vec_madd( dZ2, dZ2, D_dot_D2 );
D_dot_V1 = vec_madd( dZ1, vZ1, D_dot_V1 );
D_dot_V2 = vec_madd( dZ2, vZ2, D_dot_V2 );

// Find the reciprocal square root estimate for D dot D
// We lose a four cycles here due to data dependency stalls.
// Unrolling this loop to do 16 at a time rather than 8 would
// allow you to avoid these.
//
// result = (D dot V) * (D dot D)-0.5
results[0] = vec_madd( D_dot_V1, vec_rsqrte( D_dot_D1), zero );
results[1] = vec_madd( D_dot_V2, vec_rsqrte( D_dot_D2), zero );
results += 2;

}

}

Here is the SimG4 trace for the function: