Principe¶
Dans le papyrus Rhind, le scribe Ahmès liste 87 problèmes résolus d'arithmétique, certains depuis 2000 ans avant notre ère.
Parmi eux, une technique de multiplication qui n'a besoin que de
- multiplication par deux
- addition
Elle consiste à
- décomposer un nombre en une somme de puissances de 2
- générer une table de puissances de 2 de l'autre nombre
- additionner les nombres de cette table qui font partie de la décomposition du premier nombre
La méthode du paysan russe améliore cette approche en calculant la décomposition en ligne plutôt qu'en utilisant une table de puissance
On veut calculer le produit $a \cdot b$.
Soit la décomposition binaire de $b = \sum d_k 2^k $ avec $d_k=0$ ou $1$.
Par exemple, $42 = 2 + 8 + 32 = 0 \cdot 2^0 + 1 \cdot 2^1 + 0 \cdot 2^2 + 1 \cdot 2^3 + 0 \cdot 2^4 + 1 \cdot 2^5$
Le produit devient
$\begin{aligned} a \cdot b & = \sum a \cdot d_k \cdot 2^k \\ & = \sum d_k \cdot a_k \text{ , avec } a_k = a \cdot 2^k \end{aligned}$
Les suites $a_k$ et $d_k$ se calculent aisément récursivement:
- $\left\{ \begin{aligned} a_0 & = a \\ a_k & = 2 \cdot a_{k-1} \end{aligned} \right. $
- $\left\{ \begin{aligned} b_0 & = b \\ b_k & = b_{k-1}\ \mathbf{div}\ 2 \end{aligned} \right.$
- $d_k = b_k\ \mathbf{mod}\ 2$
Mise en oeuvre¶
Les relations récursives nous donnent le cas général. Le cas trivial est simplement $a \cdot b = 0 \ \mathbf{si}\ b = 0$. On en déduit le code
def multiplication_russe(a,b):
if b == 0: return 0
else:
if b % 2 :
return a + multiplication_russe(a+a,b//2)
else:
return multiplication_russe(a+a,b//2)
a = 743
b = 42
print("{} * {} = {} \n".format(a,b,a*b))
multiplication_russe(a,b)
Sur un ordinateur utilisant une représentation binaire des nombres, les shifts et opérations logiques bit à bit permettent de calculer plus efficacement la multiplication par 2 (a << 1), la division entière par 2 (b >> 1) et le modulo 2 (b & 1)
def multiplication_russe(a,b):
if b == 1: return a
else:
if b & 1 :
return a + multiplication_russe(a << 1,b >> 1)
else:
return multiplication_russe(a << 1,b >> 1)
a = 743
b = 42
print("{} * {} = {} \n".format(a,b,a*b))
multiplication_russe(a,b)
Cette procédure peut aussi être mise en oeuvre itérativement. Pour cela, on ajoute une troisième suite de nombres qui accumulent les éléments à sommer plutôt que d'effectuer la somme en retour de récursion.
$\left\{ \begin{aligned} r_0 & = 0 \\ r_k & = r_{k-1} + d_{k-1} \cdot a_{k-1} \\ & = \begin{cases} r_{k-1} & \mathbf{si}\ d_{k-1} = 0, \\ r_{k-1} + a_{k-1} & \mathbf{sinon} \end{cases} \end{aligned} \right.$
def multiplication_russe(a,b):
if b == 0: return 0
r = 0
while b > 0:
if b&1: r += a
a <<= 1
b >>= 1
return r
a = 743
b = 42
print("{} * {} = {} \n".format(a,b,a*b))
multiplication_russe(a,b)
Reste de la division entière¶
Cette procédure peut aussi être utilisée pour calculer les modulos en un temps logarithmique en n'utilisant que la soustraction et la multiplication par 2.
Le cas général de la récursion calcule $u = a\ \mathbf{mod}\ b$ à partir de $v = a\ \mathbf{mod}\ 2 \cdot b$. On montre aisément que
$u = \begin{cases} v & \text{ si } 0 \leq v < b \\ v - b & \text{ si } b \leq v < 2 \cdot b \end{cases}$
ce qui couvre toutes les valeurs possibles de $v$.
Le cas trivial utilise quant à lui
$u = a - b \ \ \text{ si } b \leq a < 2 \cdot b$
def reste_division_entiere(a,b):
assert a >= b and b > 0
if a - b >= b:
a = reste_division_entiere(a,b << 1) # b*2
if a < b: return a
return a-b
reste_division_entiere(743,25)
Puissance de n¶
La même technique s'utilise aussi pour calculer $b^n$ pour $n$ entier quelconque.
def puissance(b,n):
r = 1
a = b
while n > 0:
if n&1: r = r * a
a = a * a
n >>= 1
return r
for i in range(10):
print("2^{} = {}".format(i,puissance(2,i)))
Fibonacci en temps logarithmique¶
Notons enfin que cette même procédure peut être utilisée pour calculer un nombre de Fibonacci en un temps logarithmique. Pour cela, on utilise l'algèbre matriciel pour calculer
$\left( \begin{aligned} & F_{n-1} \\ & F_n \end{aligned} \right)$ à partir de $\left( \begin{aligned} & F_{n-2} \\ & F_{n-1} \end{aligned} \right)$
On obtient la relation
$\left( \begin{aligned} & F_{n-1} \\ & F_n \end{aligned} \right) = \left( \begin{aligned} 0 \ 1 \\ 1 \ 1 \end{aligned} \right) \cdot \left( \begin{aligned} & F_{n-2} \\ & F_{n-1} \end{aligned} \right)$
On peut la réécrire
$\begin{aligned} \left( \begin{aligned} & F_{n-1} \\ & F_n \end{aligned} \right) & = \left( \begin{aligned} 0 \ 1 \\ 1 \ 1 \end{aligned} \right)^{2} \cdot \left( \begin{aligned} & F_{n-3} \\ & F_{n-2} \end{aligned} \right) = \left( \begin{aligned} 0 \ 1 \\ 1 \ 1 \end{aligned} \right)^{3} \cdot \left( \begin{aligned} & F_{n-4} \\ & F_{n-3} \end{aligned} \right) = ... \\ & = \left( \begin{aligned} 0 \ 1 \\ 1 \ 1 \end{aligned} \right)^{n-1} \cdot \left( \begin{aligned} & F_{0} \\ & F_{1} \end{aligned} \right) = \left( \begin{aligned} 0 \ 1 \\ 1 \ 1 \end{aligned} \right)^{n-1} \cdot \left( \begin{aligned} & 1 \\ & 1 \end{aligned} \right) \end{aligned} $
Calculer le $n^{ieme}$ nombre de Fibonacci consiste donc simplement à calculer la $(n-1)^{ieme}$ puissance d'une matrice 2x2, ce que nous venons de voir à une complexité logarithmique.
import numpy as np
def fibonacci(n):
a = np.array([[0, 1], [1, 1]])
r = np.array([[1, 0], [0, 1]])
while n > 0:
if n&1: r = np.matmul(r,a)
a = np.matmul(a,a)
n >>= 1
return r[1,1]
for i in range(10):
print(fibonacci(i))
|
© Olivier Cuisenaire, 2018 |