求三维向量组仿射变换的快速算法。

[pulley]

对一组三维向量进行仿射变换，也就是与仿射变换矩阵相乘(4*4)。我是每个向量分别与矩阵相乘的，已经考虑到了仿射变换的特殊性，减少了一些乘法运算，可是在应用中还是不够快。
有没有快速算法解决这个问题？由衷感谢。

[短歌如风]

快速算法我没见过，不过向量乘矩阵是一个高度可并行化的操作，可以通过用汇编代码的方式来充分利用FPU的流水线功能。通过指令顺序的排布，可以达到每条指令一个时钟周期。
我用的一段代码，在D6中实现，fastcall调用方式，第一个参数：EAX,第二个参数：EDX，第三个参数：ECX
矩阵使用80位临时实数格式，向量使用64位长实数格式，你按你的需要改一下吧。

type
TMatrix = array[0..3, 0..3] of Extended;
TVector3d = packed array [0..2] of Double;

procedure VectorsMulMatrix(var FirstVector: TVector3d; Count: Integer;
const AMatrix: TMatrix);
const
SizeOfMatrixElmt = 10;
SizeOfVectorElmt = 8;
asm
@LVECT_BEGIN:

FLD TBYTE PTR [ECX];
FMUL QWORD PTR [EAX];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt];
FMUL QWORD PTR [EAX];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * 2];
FMUL QWORD PTR [EAX];

FLD TBYTE PTR [ECX + SizeOfMatrixElmt * 4];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 + 1)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 + 2)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FXCH ST(2);
FADDP ST(5), ST(0);
FADDP ST(3), ST(0);
FADDP ST(1), ST(0);

FLD TBYTE PTR [ECX + SizeOfMatrixElmt * 4 * 2];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 2 + 1)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 2 + 2)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FXCH ST(2);
FADDP ST(5), ST(0);
FADDP ST(3), ST(0);
FADDP ST(1), ST(0);

FXCH ST(2);
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 3)];
FADDP ST(1), ST(0)
FXCH ST(1);
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 3 + 1)];
FADDP ST(1), ST(0);
FXCH ST(2);
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 3 + 2)];
FADDP ST(1), ST(0);
FXCH ST(1);

FSTP QWORD PTR [EAX];
FSTP QWORD PTR [EAX + SizeOfVectorElmt * 2];
FSTP QWORD PTR [EAX + SizeOfVectorElmt];
LEA EAX, [EAX + SizeOfVectorElmt * 3];

SUB EDX, 1;
JNZ @LVECT_BEGIN;
end;

[pulley]

速度确实很快。谢谢。

我改了一下，改用C，矩阵用double。我把我改的结果也贴上来：

typedef double TMatrix[4][4];
typedef double TVector3d[3];

enum {SizeOfMatrixElmt = 8, SizeOfVectorElmt = 8};

void VectorsMulMatrix(TVector3d* FirstVector, int Count, const TMatrix * AMatrix)
{
__asm
{
LVECT_BEGIN:
MOV EAX, FirstVector;
MOV EDX, Count;
MOV ECX, AMatrix;
FLD TBYTE PTR [ECX];
FMUL QWORD PTR [EAX];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt];
FMUL QWORD PTR [EAX];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * 2];
FMUL QWORD PTR [EAX];

FLD TBYTE PTR [ECX + SizeOfMatrixElmt * 4];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 + 1)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 + 2)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FXCH ST(2);
FADDP ST(5), ST(0);
FADDP ST(3), ST(0);
FADDP ST(1), ST(0);

FLD TBYTE PTR [ECX + SizeOfMatrixElmt * 4 * 2];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 2 + 1)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 2 + 2)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FXCH ST(2);
FADDP ST(5), ST(0);
FADDP ST(3), ST(0);
FADDP ST(1), ST(0);

FXCH ST(2);
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 3)];
FADDP ST(1), ST(0)
FXCH ST(1);
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 3 + 1)];
FADDP ST(1), ST(0);
FXCH ST(2);
FLD TBYTE PTR [ECX + SizeOfMatrixElmt * (4 * 3 + 2)];
FADDP ST(1), ST(0);
FXCH ST(1);

FSTP QWORD PTR [EAX];
FSTP QWORD PTR [EAX + SizeOfVectorElmt * 2];
FSTP QWORD PTR [EAX + SizeOfVectorElmt];
LEA EAX, [EAX + SizeOfVectorElmt * 3];

SUB EDX, 1;
JNZ LVECT_BEGIN;
}
}
在VC6中编译运行成功。

[pulley]

刚才贴的是错的，应该是这个：

typedef double TMatrix[4][4];
typedef double TVector3d[3];

enum {SizeOfMatrixElmt = 8, SizeOfVectorElmt = 8};

void VectorsMulMatrix(TVector3d* FirstVector, int Count, const TMatrix * AMatrix)
{
__asm
{
LVECT_BEGIN:
MOV EAX, FirstVector;
MOV EDX, Count;
MOV ECX, AMatrix;
FLD QWORD PTR [ECX];
FMUL QWORD PTR [EAX];
FLD QWORD PTR [ECX + SizeOfMatrixElmt];
FMUL QWORD PTR [EAX];
FLD QWORD PTR [ECX + SizeOfMatrixElmt * 2];
FMUL QWORD PTR [EAX];

FLD QWORD PTR [ECX + SizeOfMatrixElmt * 4];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FLD QWORD PTR [ECX + SizeOfMatrixElmt * (4 + 1)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FLD QWORD PTR [ECX + SizeOfMatrixElmt * (4 + 2)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt];
FXCH ST(2);
FADDP ST(5), ST(0);
FADDP ST(3), ST(0);
FADDP ST(1), ST(0);

FLD QWORD PTR [ECX + SizeOfMatrixElmt * 4 * 2];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FLD QWORD PTR [ECX + SizeOfMatrixElmt * (4 * 2 + 1)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FLD QWORD PTR [ECX + SizeOfMatrixElmt * (4 * 2 + 2)];
FMUL QWORD PTR [EAX + SizeOfVectorElmt * 2];
FXCH ST(2);
FADDP ST(5), ST(0);
FADDP ST(3), ST(0);
FADDP ST(1), ST(0);

FXCH ST(2);
FLD QWORD PTR [ECX + SizeOfMatrixElmt * (4 * 3)];
FADDP ST(1), ST(0)
FXCH ST(1);
FLD QWORD PTR [ECX + SizeOfMatrixElmt * (4 * 3 + 1)];
FADDP ST(1), ST(0);
FXCH ST(2);
FLD QWORD PTR [ECX + SizeOfMatrixElmt * (4 * 3 + 2)];
FADDP ST(1), ST(0);
FXCH ST(1);

FSTP QWORD PTR [EAX];
FSTP QWORD PTR [EAX + SizeOfVectorElmt * 2];
FSTP QWORD PTR [EAX + SizeOfVectorElmt];
LEA EAX, [EAX + SizeOfVectorElmt * 3];

SUB EDX, 1;
JNZ LVECT_BEGIN;
}
}