Geekly articles weekly

适用于Delphi的嵌入式ASM上的Fast Sin和Cos

大家好！

有必要对Sin和Cos进行快速计算。 计算的基础是泰勒级数的展开。 我在3D系统（OpenGL和我开发的图形库）中使用它。 不幸的是，不可能减少Double的“理想”数目，但这可以通过良好的加速度来抵消。 该代码在Delphi XE6中的汇编器中编写。 由SSE2使用。

不适合科学计算，但完全可以用于游戏中。
准确度足以覆盖使用数量Single的长度
与矩阵相乘。

总结：

  1. 达到的精度是：10.e-13
  2. 与CPU的最大差异为0.000000000000045。
  3. 与CPU相比，速度提高了4.75倍。
  4. 与Math.Sin和Math.Cos相比，速度提高了2.6倍。

为了进行测试，我使用了Intel Core-i7 6950X Extreme 3.0 GHz处理器。
Delphi源代码嵌入在汇编程序注释中。

源代码：

扰流板方向
var gSinTab: array of Double; gSinAddr: UInt64; const AbsMask64: UInt64 = ($7FFFFFFFFFFFFFFF); PI90: Double = (PI / 2.0); FSinTab: array[0..7] of Double = (1.0 / 6.0, //3! 1.0 / 120.0, //5! 1.0 / 5040.0, //7! 1.0 / 362880.0, //9! 1.0 / 39916800.0, //11! 1.0 / 6227020800.0, //13! 1.0 / 1307674368000.0, //15! 1.0 / 355687428096000.0); //17! //  Sin    QSinMaxP: array[0..3] of Double = (0.0, 1.0, 0.0, -1.0); //  Sin    QSinMaxM: array[0..3] of Double = (0.0, -1.0, 0.0, 1.0); //     QSinSignP: array[0..3] of Double = (1.0, 1.0, -1.0, -1.0); //     QSinSignM: array[0..3] of Double = (-1.0, -1.0, 1.0, 1.0); function Sinf(const A: Double): Double; asm .NOFRAME //   A  xmm0 //   // X := IntFMod(Abs(A), PI90, D); // Result := FNum - (Trunc(FNum / FDen) * FDen); pxor xmm3, xmm3 comisd A, xmm3 jnz @ANZ ret // Result := 0.0; @ANZ: //   //Minus := (A < 0.0); setc cl movsd xmm1, [AbsMask64] //Abs(A) pand A, xmm1 movsd xmm1, [PI90] //PI90 = FDen movsd xmm2, A // A = FNum divsd xmm2, xmm1 //(FNum / FDen) roundsd xmm2, xmm2, 11b //11b - Trunc cvtsd2si rax, xmm2 //  (X  EAX) mulsd xmm2, xmm1 subsd xmm0, xmm2 //----------------------------------- //  and rax, 3 // D := (D and 3); //----------------------------------- //     // if (X = 0.0) then // begin // if Minus then // Exit(QSinMaxM[D]) else // Exit(QSinMaxP[D]); // end; comisd xmm0, xmm3 jnz @XNZ shl rax, 3 //     Double (8 ) cmp cl, 1 //  ,   je @MaxMinus @MaxPlus: lea rdx, qword ptr [QSinMaxP] movsd xmm0, qword ptr [rdx + rax] ret @MaxMinus: lea rdx, qword ptr [QSinMaxM] movsd xmm0, qword ptr [rdx + rax] ret //----------------------------------- @XNZ: //        // if (D and 1) <> 0 then X := (PI90 - X); mov edx, eax and edx, 1 cmp edx, 0 je @DZ subsd xmm1, xmm0 movsd xmm0, xmm1 //----------------------------------- @DZ: // Result   xmm0 // X  xmm0 // N := (X * X)  xmm2 // F := (N * X)  xmm1 // N movsd xmm2, xmm0 //  X mulsd xmm2, xmm2 // N := (X * X) // F movsd xmm1, xmm2 //  N mulsd xmm1, xmm0 // F := (N * X) //     mov rdx, [gSinTab] movaps xmm4, dqword ptr [rdx] movaps xmm6, dqword ptr [rdx + 16] movaps xmm8, dqword ptr [rdx + 32] movaps xmm10, dqword ptr [rdx + 48] //    movhlps xmm5, xmm4 movhlps xmm7, xmm6 movhlps xmm9, xmm8 movhlps xmm11, xmm10 // Result := X - F * PDouble(gSinAddr)^; mulsd xmm4, xmm1 //FSinTab[0] subsd xmm0, xmm4 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 8)^; mulsd xmm1, xmm2 mulsd xmm5, xmm1 //FSinTab[1] addsd xmm0, xmm5 // F := (F * N); // Result := Result - F * PDouble(gSinAddr + 16)^; mulsd xmm1, xmm2 mulsd xmm6, xmm1 //FSinTab[2] subsd xmm0, xmm6 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 24)^; mulsd xmm1, xmm2 mulsd xmm7, xmm1 //FSinTab[3] addsd xmm0, xmm7 // F := (F * N); // Result := Result - F * PDouble(gSinAddr + 32)^; mulsd xmm1, xmm2 mulsd xmm8, xmm1 //FSinTab[4] subsd xmm0, xmm8 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 40)^; mulsd xmm1, xmm2 mulsd xmm9, xmm1 //FSinTab[5] addsd xmm0, xmm9 // F := (F * N); // Result := Result - F * PDouble(gSinAddr + 48)^; mulsd xmm1, xmm2 mulsd xmm10, xmm1 //FSinTab[6] subsd xmm0, xmm10 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 56)^; mulsd xmm1, xmm2 mulsd xmm11, xmm1 //FSinTab[7] addsd xmm0, xmm11 //----------------------------------- // if Minus then // Result := (Result * QSinSignM[D]) else // Result := (Result * QSinSignP[D]); shl rax, 3 //  8 cmp cl, 1 je @Minus lea rdx, qword ptr [QSinSignP] mulsd xmm0, qword ptr [rdx + rax] ret @Minus: lea rdx, qword ptr [QSinSignM] mulsd xmm0, qword ptr [rdx + rax] end; //------------------------------------ function Cosf(const A: Double): Double; asm .NOFRAME // A  xmm0 //   // X := IntFMod(AMod, PI90, D); // Result := FNum - (Trunc(FNum / FDen) * FDen); pxor xmm3, xmm3 comisd A, xmm3 jnz @ANZ movsd xmm0, [DoubleOne] ret // Result := 1.0; //----------------------------------- @ANZ: //   //Minus := (A < 0.0); setc cl movsd xmm1, [AbsMask64] //Abs(A) pand A, xmm1 //----------------------------------- movsd xmm1, [PI90] //PI90 = FDen //----------------------------------- // if Minus then // AMod := Abs(A) - PI90 else // AMod := Abs(A) + PI90; cmp cl, 1 je @SubPI90 addsd A, xmm1 jmp @IntFMod @SubPI90: subsd A, xmm1 //----------------------------------- @IntFMod: movsd xmm2, A // A = FNum divsd xmm2, xmm1 //(FNum / FDen) roundsd xmm2, xmm2, 11b //11b - Trunc cvtsd2si rax, xmm2 //  (X  EAX) mulsd xmm2, xmm1 subsd xmm0, xmm2 //----------------------------------- //  and rax, 3 // D := (D and 3); //----------------------------------- //     // if (X = 0.0) then // begin // if Minus then // Exit(QSinMaxM[D]) else // Exit(QSinMaxP[D]); // end; comisd xmm0, xmm3 jnz @XNZ shl rax, 3 //     Double (8 ) cmp cl, 1 //  ,   je @MaxMinus @MaxPlus: lea rdx, qword ptr [QSinMaxP] movsd xmm0, qword ptr [rdx + rax] ret @MaxMinus: lea rdx, qword ptr [QSinMaxM] movsd xmm0, qword ptr [rdx + rax] ret //----------------------------------- @XNZ: //        // if (D and 1) <> 0 then X := (PI90 - X); mov edx, eax and edx, 1 cmp edx, 0 je @DZ subsd xmm1, xmm0 movsd xmm0, xmm1 @DZ: // Result   xmm0 // X  xmm0 // N := (X * X)  xmm2 // F := (N * X)  xmm1 // N movsd xmm2, xmm0 //  X mulsd xmm2, xmm2 // N := (X * X) // F movsd xmm1, xmm2 //  N mulsd xmm1, xmm0 // F := (N * X) //     mov rdx, [gSinTab] movaps xmm4, dqword ptr [rdx] movaps xmm6, dqword ptr [rdx + 16] movaps xmm8, dqword ptr [rdx + 32] movaps xmm10, dqword ptr [rdx + 48] //    movhlps xmm5, xmm4 movhlps xmm7, xmm6 movhlps xmm9, xmm8 movhlps xmm11, xmm10 // Result := X - F * PDouble(gSinAddr)^; mulsd xmm4, xmm1 //FSinTab[0] subsd xmm0, xmm4 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 8)^; mulsd xmm1, xmm2 mulsd xmm5, xmm1 //FSinTab[1] addsd xmm0, xmm5 // F := (F * N); // Result := Result - F * PDouble(gSinAddr + 16)^; mulsd xmm1, xmm2 mulsd xmm6, xmm1 //FSinTab[2] subsd xmm0, xmm6 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 24)^; mulsd xmm1, xmm2 mulsd xmm7, xmm1 //FSinTab[3] addsd xmm0, xmm7 // F := (F * N); // Result := Result - F * PDouble(gSinAddr + 32)^; mulsd xmm1, xmm2 mulsd xmm8, xmm1 //FSinTab[4] subsd xmm0, xmm8 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 40)^; mulsd xmm1, xmm2 mulsd xmm9, xmm1 //FSinTab[5] addsd xmm0, xmm9 // F := (F * N); // Result := Result - F * PDouble(gSinAddr + 48)^; mulsd xmm1, xmm2 mulsd xmm10, xmm1 //FSinTab[6] subsd xmm0, xmm10 // F := (F * N); // Result := Result + F * PDouble(gSinAddr + 56)^; mulsd xmm1, xmm2 mulsd xmm11, xmm1 //FSinTab[7] addsd xmm0, xmm11 //----------------------------------- // if Minus then // Result := (Result * QSinSignM[D]) else // Result := (Result * QSinSignP[D]); shl rax, 3 //  8 cmp cl, 1 je @Minus lea rdx, qword ptr [QSinSignP] mulsd xmm0, qword ptr [rdx + rax] ret @Minus: lea rdx, qword ptr [QSinSignM] mulsd xmm0, qword ptr [rdx + rax] end; Initialization //  SetLength(gSinTab, 8); //  gSinAddr := UInt64(@FSinTab[0]); //     Move(FSinTab[0], gSinTab[0], SizeOf(Double) * 8); Finalization SetLength(gSinTab, 0);


这里的工作例子

欢迎提出建设性的建议和意见。

Source: https://habr.com/ru/post/zh-CN429786/

More articles:


All Articles