求助！D3DXMatrixRotationAxis对应的数学公式是什么？

Question

由于某些原因，目前不能使用D3D中现成的API，但是现在需要用到D3DXMatrixRotationAxis这个函数的功能，不过几代只记住了些基本的东西，这个带上旋转完全把自己搞晕了。
D3DXMATRIX *D3DXMatrixRotationAxis(          D3DXMATRIX *pOut,
    CONST D3DXVECTOR3 *pV,
    FLOAT Angle
);第一个参数是4x4的矩阵，第二个参数是三维向量，第三个参数是弧度。功能是根据pV轴旋转Angle角度，得到pOut这个旋转矩阵。
如果有哪位兄台数学比较强，希望能解答一下。另外我把它反汇编了，已经可以用汇编代码搞定，但是还是想知道公式。由于对浮点计算部分汇编可以说是完全不懂，无法从代码看出公式，下面是我反汇编后的代码，如果哪位大侠的汇编比较强悍，可以看出公式，也请不吝赐教。
float* __stdcall Vec3Normalize(float *pOut, float *pV)  //这个是三维向量单位化，这个就不汇编了，它会被下面的函数调用到
{
float rr = 1.0f/sqrt(pV[0] * pV[0] + pV[1] * pV[1] + pV[2] * pV[2]);
pOut[0] = pV[0] * rr, pOut[1] = pV[1] * rr, pOut[2] = pV[2] * rr;
return pOut;
}__declspec(naked) Matrix16* __stdcall MatrixRotationAxis(Matrix16 *pOut, float *pV, float Angle)  //这个是D3DXMatrixRotationAxis反汇编的代码
{
_asm
{
push        ebp  
mov         ebp,esp 
sub         esp,1Ch 
lea         eax,[ebp-4] 
mov         dword ptr [ebp-0Ch],eax 
push        esi  
lea         eax,[ebp-8] 
push        edi  
mov         dword ptr [ebp-10h],eax 
mov         eax,dword ptr [ebp-10h] 
mov         edx,dword ptr [ebp-0Ch] 
fld         dword ptr [ebp+10h] 
fsincos          
fstp        dword ptr [edx] 
fstp        dword ptr [eax] 
fld1             
mov         esi,dword ptr [ebp+0Ch] 
fsub        dword ptr [ebp-4] 
lea         edi,[ebp-1Ch] 
movs        dword ptr [edi],dword ptr [esi] 
lea         eax,[ebp-1Ch] 
fstp        dword ptr [ebp+10h] 
movs        dword ptr [edi],dword ptr [esi] 
push        eax  
lea         eax,[ebp-1Ch] 
push        eax  
movs        dword ptr [edi],dword ptr [esi] 
call        Vec3Normalize         //本来是call D3DXVec3Normalize的
mov         eax,dword ptr [ebp+8] 
fld         dword ptr [ebp+10h] 
fld         dword ptr [ebp-18h] 
fmul        dword ptr [ebp-1Ch] 
fmul        st,st(1) 
fstp        dword ptr [ebp+10h] 
fld         dword ptr [ebp-14h] 
fld         dword ptr [ebp-18h] 
fmul        st,st(1) 
fmul        st,st(2) 
fld         st(1) 
fmul        dword ptr [ebp-1Ch] 
fmul        st,st(3) 
fld         dword ptr [ebp-1Ch] 
fmul        dword ptr [ebp-1Ch] 
fmul        st,st(4) 
fadd        dword ptr [ebp-4] 
fstp        dword ptr [eax] 
fld         st(2) 
fmul        dword ptr [ebp-8] 
fld         dword ptr [ebp+10h] 
fadd        st,st(1) 
fstp        dword ptr [eax+4] 
fld         dword ptr [ebp-18h] 
fmul        dword ptr [ebp-8] 
fst         dword ptr [ebp+0Ch] 
fld         st(2) 
fsub        st,st(1) 
fstp        dword ptr [eax+8] 
fstp        st(0) 
fldz             
fstp        dword ptr [eax+0Ch] 
fld         dword ptr [ebp+10h] 
fsub        st,st(1) 
fstp        dword ptr [eax+10h] 
fstp        st(0) 
fld         dword ptr [ebp-18h] 
fmul        dword ptr [ebp-18h] 
fmul        st,st(4) 
fadd        dword ptr [ebp-4] 
fstp        dword ptr [eax+14h] 
fld         dword ptr [ebp-1Ch] 
fmul        dword ptr [ebp-8] 
fld         st(0) 
fadd        st,st(3) 
fstp        dword ptr [eax+18h] 
fldz             
fstp        dword ptr [eax+1Ch] 
fld         dword ptr [ebp+0Ch] 
fadd        st,st(2) 
fstp        dword ptr [eax+20h] 
fld         st(2) 
fsub        st,st(1) 
fstp        dword ptr [eax+24h] 
fstp        st(0) 
fstp        st(0) 
fstp        st(0) 
fld         st(0) 
fmul        st,st(1) 
fmul        st,st(2) 
fadd        dword ptr [ebp-4] 
fstp        dword ptr [eax+28h] 
pop         edi  
pop         esi  
fstp        st(0) 
fstp        st(0) 
fldz             
fstp        dword ptr [eax+2Ch] 
fldz             
fstp        dword ptr [eax+30h] 
fldz             
fstp        dword ptr [eax+34h] 
fldz             
fstp        dword ptr [eax+38h] 
fld1             
fstp        dword ptr [eax+3Ch] 
leave            
ret         0Ch
}
}

Answer 1

算了，我已经不对这个算法报希望了。如果有人能提供gcc和vc内嵌汇编代码的转换工具，或者更简单一点，告诉我__declspec(naked)在gcc里面对应的写法的话，一样给分。谢谢了

Answer 2

在gameres的whb9633的帮助下，已经知道了计算公式：D3DXMATRIXA16* MyMatrixRotationAxis(D3DXMATRIXA16* mat, D3DXVECTOR3* vec, float t)
{
float cost = cos(t), sint = sin(t), one_sub_cost = 1 - cost;
D3DXVECTOR3 v(*vec);
D3DXVec3Normalize(&v, &v);
*mat = D3DXMATRIXA16(v.x * v.x * one_sub_cost + cost, v.x * v.y * one_sub_cost + v.z * sint, v.x * v.z * one_sub_cost - v.y * sint, 0, \
v.x * v.y * one_sub_cost - v.z * sint, v.y * v.y * one_sub_cost + cost, v.y * v.z * one_sub_cost + v.x * sint, 0, \
v.x * v.z * one_sub_cost + v.y * sint, v.y * v.z * one_sub_cost - v.x * sint, v.z * v.z * one_sub_cost + cost, 0, \
0, 0, 0, 1);
return mat;
}