ejemplos - pow c++
¿La forma más precisa de hacer una operación combinada de multiplicar y dividir en 64 bits? (9)
¿Cuál es la forma más precisa que puedo hacer una operación de multiplicar y dividir para enteros de 64 bits que funciona tanto en programas de 32 bits como de 64 bits (en Visual C ++)? (En caso de desbordamiento, necesito el resultado mod 2 64 ).
(Estoy buscando algo como MulDiv64 , excepto que este usa ensamblado en línea, que solo funciona en programas de 32 bits).
Obviamente, lanzar al double y hacia atrás es posible, pero me pregunto si hay una manera más precisa que no sea demasiado complicada. (es decir, no estoy buscando bibliotecas aritméticas de precisión arbitraria aquí).
¿Tiene el tipo COMP (tipo entero de 64 bits basado en x87) a su disposición en VC ++? Lo he usado ocasionalmente en Delphi en el pasado cuando necesitaba matemáticas enteras de 64 bits. Durante años, fue mucho más rápido que la matemática entera de 64 bits basada en la biblioteca, sin duda cuando se trataba de división.
En Delphi 2007 (el último que he instalado - 32 bits), implementaría MulDiv64 de esta manera:
function MulDiv64(const a1, a2, a3: int64): int64;
var
c1: comp absolute a1;
c2: comp absolute a2;
c3: comp absolute a3;
r: comp absolute result;
begin
r := c1*c2/c3;
end;
(Esas extrañas declaraciones absolutas superponen las variables de comp sobre sus contrapartes de enteros de 64 bits. Hubiera utilizado moldes de tipos simples, excepto que el compilador Delphi se confunde con eso, probablemente porque el lenguaje Delphi (o como lo llamen ahora) tiene no hay una clara distinción sintáctica entre la conversión de tipo (reinterpretación) y la conversión de tipo de valor).
De todos modos, Delphi 2007 presenta lo anterior de la siguiente manera:
0046129C 55 push ebp
0046129D 8BEC mov ebp,esp
0046129F 83C4F8 add esp,-$08
004612A2 DF6D18 fild qword ptr [ebp+$18]
004612A5 DF6D10 fild qword ptr [ebp+$10]
004612A8 DEC9 fmulp st(1)
004612AA DF6D08 fild qword ptr [ebp+$08]
004612AD DEF9 fdivp st(1)
004612AF DF7DF8 fistp qword ptr [ebp-$08]
004612B2 9B wait
004612B3 8B45F8 mov eax,[ebp-$08]
004612B6 8B55FC mov edx,[ebp-$04]
004612B9 59 pop ecx
004612BA 59 pop ecx
004612BB 5D pop ebp
004612BC C21800 ret $0018
La siguiente instrucción arroja 256204778801521550, que parece ser correcta.
writeln(MulDiv64($aaaaaaaaaaaaaaa, $555555555555555, $1000000000000000));
Si desea implementar esto como un ensamblado en línea de VC ++, es posible que necesite hacer algunos ajustes en los indicadores de redondeo predeterminados para lograr lo mismo, no lo sé, no he tenido la necesidad de averiguarlo, todavía :)
Aquí hay un método de aproximación que puede usar: (precisión completa a menos que a> 0x7FFFFFFF o b> 0x7FFFFFFF yc sea mayor que aob)
constexpr int64_t muldiv(int64_t a, int64_t b, int64_t c, unsigned n = 0) {
return (a < 0x7FFFFFFF && b < 0x7FFFFFFF) ? (a * b) / c : (n != 2) ? (c <= a) ? ((a / c) * b + muldiv(b, a % c, c, n + 1)) : muldiv(a, b, c / 2) / 2 : 0;
}
El módulo se usa para encontrar la pérdida de precisión, que luego se vuelve a conectar a la función. Esto es similar al algoritmo de división clásico.
2 fue elegido porque (x % x) % x = x % x .
Bueno, puedes cortar tus operandos de 64 bits en trozos de 32 bits (parte baja y alta). Que hacer la operación en ellos que quieras. Todos los resultados intermedios serán de menos de 64 bits y, por lo tanto, pueden almacenarse en los tipos de datos que tenga.
Considere que quiere multiplicar a by b y luego dividir por d :
uint64_t LossyMulDiv64(uint64_t a, uint64_t b, uint64_t d)
{
long double f = long double(b)/d;
uint64_t highPart = uint64_t((a & ~0xffffffff) * f + 0.5);
uint64_t lowPart = uint64_t((a & 0xffffffff) * f + 0.5);
return highPart + lowPart;
}
Este código divide el valor de a en partes más altas y más bajas de 32 bits, luego multiplica las partes de 32 bits por separado por una relación precisa de 52 bits de b a d , redondea las multiplicaciones parciales y las resume de nuevo en un entero. Todavía se pierde cierta precisión, pero el resultado es más preciso que el simple return a * double(b) / d; .
Esta es una respuesta wiki de la comunidad, ya que en realidad se trata de un puñado de indicadores de otros documentos / referencias (no puedo publicar el código relevante).
La multiplicación de dos entradas de 64 bits a un resultado de 128 bits es bastante fácil mediante una aplicación directa de la técnica de lápiz y papel que todo el mundo aprende en la escuela primaria.
El comentario de GregS es correcto: Knuth cubre la división en "El Arte de la Programación de Computadora, Segunda Edición, Volumen 2 / Algoritmos Seminumerical" al final de la Sección 4.3.1 Aritmética de Precisión Múltiple / Los Algoritmos Clásicos (páginas 255 - 265 en mi copia). No es una lectura fácil, al menos no para alguien como yo que se olvidó de la mayoría de las matemáticas más allá del Álgebra del 7º grado. Justo antes, Knuth también cubre el lado de la multiplicación de las cosas.
Algunas otras opciones para ideas (estas notas son para algoritmos de división, pero la mayoría también habla de multiplicación):
- Jack Crenshaw cubre los algoritmos de la división Knuth de una manera más legible en una serie de artículos de la revista Embedded System Programming 1997 (desafortunadamente, mis notas no tienen los problemas exactos). Lamentablemente, los artículos de viejos problemas ESP no son fáciles de encontrar en línea. Si tiene acceso a una biblioteca de la Universidad, tal vez haya algunos números atrasados o una copia de la Biblioteca de CD-ROM de ESP disponible para usted.
- Thomas Rodeheffer de investigación de Microsoft tiene un documento sobre la división de software entero: http://research.microsoft.com/pubs/70645/tr-2008-141.pdf
- El artículo de Karl Hasselström sobre "División rápida de enteros grandes": http://www.treskal.com/kalle/exjobb/original-report.pdf
- El "Arte del lenguaje ensamblador" de Randall Hyde (http://webster.cs.ucr.edu/AoA/Windows/HTML/AoATOC.html), específicamente el Volumen cuatro Sección 4.2.5 (División de precisión extendida): http://webster.cs.ucr.edu/AoA/Windows/HTML/AdvancedArithmetica2.html#998729 esto está en la variante de Hyde del lenguaje ensamblador x86, pero también hay algún pseudocódigo y suficiente explicación para transferir el algoritmo a C. También es lento: realizar el división bit a bit ...
No necesita aritmética de precisión arbitraria para eso. Solo necesitas aritmética de 128 bits. Es decir, necesita una multiplicación de 64 * 64 = 128 y una división de 128/64 = 64 (con un comportamiento de desbordamiento adecuado). Esto no es tan difícil de implementar manualmente.
Para el código de modo de 64 bits puede implementar la multiplicación 64 * 64 = 128 de forma similar a la implementación de la división 128/64 = 64: 64 aquí .
Para el código de 32 bits será más complejo porque no hay instrucciones de la CPU que harían la multiplicación o división de dichos operandos largos en el modo de 32 bits y tendrás que combinar varias multiplicaciones más pequeñas en una más grande y reimplementar la división larga.
Puede usar el código de esta respuesta como marco para construir divisiones largas.
Por supuesto, si tus divisores son siempre inferiores a 2 32 (o mejor aún 2 16 ), puedes dividir mucho más rápido en modo de 32 bits encadenando varias divisiones (divide los 64 (o 32) bits más significativos del dividendo por el divisor de 32 bits (o 16 bits), luego combina el resto de esta división con los siguientes 64 (32) bits del dividendo y divídelo por el divisor y sigue haciéndolo hasta que consumas todo el dividendo. Además, si el divisor es grande pero puede tenerse en cuenta en números suficientemente pequeños, dividir por sus factores usando esta división encadenada será mejor que la solución de bucle clásica.
Solo necesitas enteros de 64 bits. Hay algunas operaciones redundantes, pero eso permite usar 10 como base y paso en el depurador.
uint64_t const base = 1ULL<<32;
uint64_t const maxdiv = (base-1)*base + (base-1);
uint64_t multdiv(uint64_t a, uint64_t b, uint64_t c)
{
// First get the easy thing
uint64_t res = (a/c) * b + (a%c) * (b/c);
a %= c;
b %= c;
// Are we done?
if (a == 0 || b == 0)
return res;
// Is it easy to compute what remain to be added?
if (c < base)
return res + (a*b/c);
// Now 0 < a < c, 0 < b < c, c >= 1ULL
// Normalize
uint64_t norm = maxdiv/c;
c *= norm;
a *= norm;
// split into 2 digits
uint64_t ah = a / base, al = a % base;
uint64_t bh = b / base, bl = b % base;
uint64_t ch = c / base, cl = c % base;
// compute the product
uint64_t p0 = al*bl;
uint64_t p1 = p0 / base + al*bh;
p0 %= base;
uint64_t p2 = p1 / base + ah*bh;
p1 = (p1 % base) + ah * bl;
p2 += p1 / base;
p1 %= base;
// p2 holds 2 digits, p1 and p0 one
// first digit is easy, not null only in case of overflow
uint64_t q2 = p2 / c;
p2 = p2 % c;
// second digit, estimate
uint64_t q1 = p2 / ch;
// and now adjust
uint64_t rhat = p2 % ch;
// the loop can be unrolled, it will be executed at most twice for
// even bases -- three times for odd one -- due to the normalisation above
while (q1 >= base || (rhat < base && q1*cl > rhat*base+p1)) {
q1--;
rhat += ch;
}
// subtract
p1 = ((p2 % base) * base + p1) - q1 * cl;
p2 = (p2 / base * base + p1 / base) - q1 * ch;
p1 = p1 % base + (p2 % base) * base;
// now p1 hold 2 digits, p0 one and p2 is to be ignored
uint64_t q0 = p1 / ch;
rhat = p1 % ch;
while (q0 >= base || (rhat < base && q0*cl > rhat*base+p0)) {
q0--;
rhat += ch;
}
// we don''t need to do the subtraction (needed only to get the remainder,
// in which case we have to divide it by norm)
return res + q0 + q1 * base; // + q2 *base*base
}
Como esto está etiquetado con Visual C ++, daré una solución que abuse de los intrínsecos específicos de MSVC.
Este ejemplo es bastante complicado. Es una versión altamente simplificada del mismo algoritmo que utilizan GMP y java.math.BigInteger para grandes divisiones.
Aunque tengo un algoritmo más simple en mente, probablemente sea 30 veces más lento.
Esta solución tiene las siguientes limitaciones / comportamiento:
- Requiere x64. No se compilará en x86.
- El cociente no es cero.
- El cociente se satura si desborda 64 bits.
Tenga en cuenta que esto es para el caso de entero sin signo. Es trivial construir una envoltura alrededor de esto para que funcione también para casos firmados. Este ejemplo también debería producir resultados truncados correctamente.
Este código no está completamente probado. Sin embargo, ha superado todos los casos de prueba que he lanzado.
(Incluso casos que intencionalmente he construido para intentar romper el algoritmo).
#include <intrin.h>
uint64_t muldiv2(uint64_t a, uint64_t b, uint64_t c){
// Normalize divisor
unsigned long shift;
_BitScanReverse64(&shift,c);
shift = 63 - shift;
c <<= shift;
// Multiply
a = _umul128(a,b,&b);
if (((b << shift) >> shift) != b){
cout << "Overflow" << endl;
return 0xffffffffffffffff;
}
b = __shiftleft128(a,b,shift);
a <<= shift;
uint32_t div;
uint32_t q0,q1;
uint64_t t0,t1;
// 1st Reduction
div = (uint32_t)(c >> 32);
t0 = b / div;
if (t0 > 0xffffffff)
t0 = 0xffffffff;
q1 = (uint32_t)t0;
while (1){
t0 = _umul128(c,(uint64_t)q1 << 32,&t1);
if (t1 < b || (t1 == b && t0 <= a))
break;
q1--;
// cout << "correction 0" << endl;
}
b -= t1;
if (t0 > a) b--;
a -= t0;
if (b > 0xffffffff){
cout << "Overflow" << endl;
return 0xffffffffffffffff;
}
// 2nd reduction
t0 = ((b << 32) | (a >> 32)) / div;
if (t0 > 0xffffffff)
t0 = 0xffffffff;
q0 = (uint32_t)t0;
while (1){
t0 = _umul128(c,q0,&t1);
if (t1 < b || (t1 == b && t0 <= a))
break;
q0--;
// cout << "correction 1" << endl;
}
// // (a - t0) gives the modulus.
// a -= t0;
return ((uint64_t)q1 << 32) | q0;
}
Tenga en cuenta que si no necesita un resultado perfectamente truncado, puede eliminar completamente el último ciclo. Si hace esto, la respuesta será no más de 2 más grande que el cociente correcto.
Casos de prueba:
cout << muldiv2(4984198405165151231,6132198419878046132,9156498145135109843) << endl;
cout << muldiv2(11540173641653250113, 10150593219136339683, 13592284235543989460) << endl;
cout << muldiv2(449033535071450778, 3155170653582908051, 4945421831474875872) << endl;
cout << muldiv2(303601908757, 829267376026, 659820219978) << endl;
cout << muldiv2(449033535071450778, 829267376026, 659820219978) << endl;
cout << muldiv2(1234568, 829267376026, 1) << endl;
cout << muldiv2(6991754535226557229, 7798003721120799096, 4923601287520449332) << endl;
cout << muldiv2(9223372036854775808, 2147483648, 18446744073709551615) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 9223372036854775807) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 12) << endl;
cout << muldiv2(18446744073709551615, 18446744073709551615, 9223372036854775808) << endl;
Salida:
3337967539561099935
8618095846487663363
286482625873293138
381569328444
564348969767547451
1023786965885666768
11073546515850664288
1073741824
9223372032559808512
Overflow
18446744073709551615
Overflow
18446744073709551615