Comment puis-je multiplier et diviser avec précision les ints 64 bits?

J’ai une fonction C:

int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d) { /* should return (a * b * c)/d */ } 

Il est possible que a soit proche de INT64_MAX, mais que le résultat final ne déborde pas, par exemple si b = 1, c = d = 40. Cependant, je ne parviens pas à comprendre comment le calculer de manière à ne jamais perdre de données. arrondir (en faisant la division en premier) ou avoir un dépassement de résultat intermédiaire.

Si j’avais access à un type de données suffisamment grand pour contenir le produit entier de a, b et c, je ferais simplement les calculs dans ce type, puis les tronquerais, mais y a-t-il moyen de le faire sans grands entiers?

Ecrivez a = q*d + r avec |r| < |d| |r| < |d| (Je suppose que d != 0 , sinon le calcul n'a pas de sens quand même). Alors (a*b*c)/d = q*b*c + (r*b*c)/d . Si q*b*c déborde, le calcul entier débordera de toute façon, vous ne vous en souciez donc pas ou vous devez vérifier le débordement. r*b*c peut encore déborder, nous utilisons donc à nouveau la même méthode pour éviter les débordements,

 int64_t q = a/d, r = a%d; int64_t part1 = q*b*c; int64_t q1 = (r*b)/d, r1 = (r*b)%d; return part1 + q1*c + (r1*c)/d; 

Je suggérerais de trouver le plus grand commun diviseur de d avec chacun de a, b et c, en divisant les facteurs communs au fur et à mesure:

 common = gcd(a,d) // You can implement GCD using Euclid's algorithm a=a/common d=d/common common = gcd(b,d) b=b/common d=d/common common = gcd(c,d) c=c/common d=d/common 

Calculez ensuite a*b*c/d supprimant tous les facteurs communs. L’algorithme GCD d’Euclid s’exécute en temps logarithmique, donc cela devrait être assez efficace.

Il est facile de voir que certaines entrées produisent des sorties qui ne peuvent pas être représentées par le type de retour int64_t . Par exemple, fn(INT64_MAX, 2, 1, 1) . Cependant, l’approche suivante devrait vous permettre de renvoyer une réponse correcte pour toute combinaison d’entrées se int64_t dans la plage de int64_t .

 int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d) { /* find the integer and remainder portions of a/d */ int64_t leftI = a / d; int64_t leftR = a % d; /* multiply the integer portion of the result by b and c */ int64_t resultI = leftI * b * c; /* multiply the remainder portion by b */ int64_t resultR = leftR * b; resultI = resultI + (resultR / d) * c; /* multiply the remainder portion by c */ resultR = (resultR % d) * c; return resultI + (resultR / d); } 

Si vous travaillez sur x86_64, l’asm prend en charge les entiers 128 bits:

 int64_t fn(uint64_t a, uint64_t b, uint64_t c, uint64_t d) { asm ( "mulq %1\n" // a *= b "movq %%rbx, %%rdx\n"// rbx = upper 64 bit of the multiplication "mulq %2\n" // multiply the lower 64 bits by c "push %%rax\n" // temporarily save the lowest 64 bits on the stack "mov %%rcx, %%rdx\n" // rcx = upper 64 bits of the multiplication "movq %%rax, %%rbx\n"// "mulq %2\n" // multiply the upper 64 bits by c "addq %%rax, %%rcx\n"// combine the middle 64 bits "addcq %%rdx, $0\n" // transfer carry tp the higest 64 bits if present "divq %3\n" // divide the upper 128 (of 192) bits by d "mov %%rbx, %%rax\n" // rbx = result "pop %%rax\n" "divq %3\n" // divide remainder:lower 64 bits by d : "+a" (a) // assigns a to rax register as in/out , "+b" (b) // assigns b to rbx register : "g" (c) // assigns c to random register , "g" (d) // assigns d to random register : "edx", "rdx" // tells the comstackr that edx/rdx will be used internally, but does not need any input ); // b now holds the upper 64 bit if (a * b * c / d) > UINT64_MAX return a; } 

Veuillez noter que tous les entiers en entrée doivent avoir la même longueur. La longueur de travail sera le double de l’entrée. Fonctionne avec non signé seulement.

Les instructions natives div et mul sur x86 fonctionnent exactement en double longueur pour permettre les débordements. Malheureusement, je ne suis pas au courant d’un compilateur insortingnsèque pour les utiliser.