Geekly articles weekly

Fast Sin y Cos en ASM integrado para Delphi

Hola a todos!

Era necesario escribir un cálculo rápido de Sin y Cos. La base para los cálculos fue la expansión de la serie Taylor. Lo uso en sistemas 3D (OpenGL y la biblioteca gráfica de mi desarrollo). Desafortunadamente, es imposible reducir el número de "ideal" para Double, pero esto se compensa con una buena aceleración. El código está escrito en el ensamblador en Delphi XE6. Usado por SSE2.

No apto para la informática científica, pero completamente para usar en juegos.
La precisión es suficiente para cubrir la longitud del número de Single, que se utiliza
para multiplicar por una matriz.

En resumen:

  1. La precisión alcanzada es: 10.e-13
  2. La discrepancia máxima con la CPU es 0.000000000000045.
  3. La velocidad aumenta en comparación con la CPU en 4.75 veces.
  4. La velocidad aumenta en comparación con Math.Sin y Math.Cos en 2.6 veces.

Para la prueba utilicé un procesador Intel Core-i7 6950X Extreme 3.0 GHz.
El código fuente de Delphi está incrustado en los comentarios del ensamblador.

Código fuente:

Encabezado de spoiler
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);


Ejemplo de trabajo aquí

Sugerencias constructivas y comentarios son bienvenidos.

Source: https://habr.com/ru/post/es429786/

More articles:


All Articles