1 / 35

Multiprécision

Multiprécision. Intérêt. Pour les algorithmes de cryptographie à clés publiques (RSA, ECC, LUC, El-Gamal …), on a besoin de nombres premiers de plusieurs centaines de chiffres Une vérification même approchée de primalité nécessite des additions, multiplications et divisions en multiprécision.

clare-gilliam
Download Presentation

Multiprécision

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Multiprécision

  2. Intérêt • Pour les algorithmes de cryptographie à clés publiques (RSA, ECC, LUC, El-Gamal …), on a besoin de nombres premiers de plusieurs centaines de chiffres • Une vérification même approchée de primalité nécessite des additions, multiplications et divisions en multiprécision. • Les processeurs possède des unités arithmétiques entières ou flottantes de taille fixe

  3. Comment éviter ce cours ? • La plupart des algorithmes « scolaires » sont suffisant pour des nombres de quelques dizaines de chiffres et deviennent vraiment obsolètes au delà de 100 chiffres • Il existe des librairies de multiprécision toutes faites dans la plupart des languages de programmation

  4. Notion de complexité • On suppose que l’on dispose d’une machine capable de faire des opérations dites élémentaires. • Il s’agit d’opérations en temps fixe (Exemple: addition, soustraction en 32 ou 64 bits). • On suppose aussi que cette machine dispose d’une mémoire où est stockée les données d’un problème à résoudre. • La taille du problème, notée T, est le nombre de bits nécessaires au stockage. • La complexité du problème, notée C, est une fonction telle qu’un problème de taille T, nécessite au plus C(T) opérations élémentaires pour être résolu.

  5. n 1 5 10 50 n2/3.log(n) 0 4,66 10,53 51,72 n.log(n) 0 8,05 23,03 195,60 n1,58 1 12,72 38,02 483,48 2n 2 32 1024 1,12e15 Notion de complexité • Plutôt que de calculer exactement la fonction C, on se contente le plus souvent d’en donner un équivalent. • Exemple: l’addition scolaire a une complexité linéaire, le tri de n entiers a une complexité en n.log(n), l’algorithme de Meissel-Lehmer a une complexité en n2/3.log(n), l’algorithme de Karatsuba a une complexité en n1,58, l’algorithme de Chernikova a une complexité exponentielle. Comparaison de complexité

  6. Addition scolaire • A=(an-1an-2 … a1a0) représente A= ai Dioù D=232 ou 264 (suivant la taille des opérateurs) • B=(bn-1bn-2 … b1b0). Soit S=(snsn-1 … s1s0) avec S=A+B • On calcule par récurrence c0=0, sk=ak+bk+ck-1(mod. D) et ck=1 si bk+ck-1> D-1-ak et 0 sinon. • Le modulo se calcule automatiquement et le test se calcule en nombres inférieurs à D • La complexité est linéaire (le nombre d’opérations est inférieur à 4n)

  7. Addition: complexité • L’addition scolaire est linéaire • Quand on ajoute A=(an-1an-2 … a1a0) et B=(bn-1bn-2 … b1b0), la moindre modification d’un ak ou un bk change le résultat final • Un algorithme d’addition doit « connaître » tous les ak et bk. • L’algorithme ne peut pas être moins que linéaire

  8. Multiplication scolaire • On veut multiplier A=(an-1an-2 … a1a0) et B=(bn-1bn-2 … b1b0) • On forme les produits partiels A.bk et on les ajoute en les décalant (les produits albk sont à calculer en précision double  4 multiplications et 3 additions simples) • Le coût d’un calcul A.bk est linéaire. • On doit en calculer n et on doit les additionner • Le coût global est en n2.

  9. Méthode de Karatsuba (1962) • On veut multiplier A=(a2n-1a2n-2 … a1a0) et B=(b2n-1b2n-2 … b1b0) • On a A=A1.Dn+A0 et B=B1Dn+B0 avec A1=(a2n-1a2n-2 … an+1an), A0=(an-1an-2 … a1a0), B1=(b2n-1b2n-2 … bn+1bn), B0=(bn-1bn-2 … b1b0). • M=A.B=M2D2n+M1Dn+M0 avec M0=A0.B0, M1=A0.B1+A1.B0, M2=A1.B1. • On écrit M1=(A0+A1).(B0+B1)-M0-M2 • On a besoin de 3 multiplications au lieu de 4 dans le cas de la multiplication scolaire. • On continue récursivement (si le nombre de mots est impair, on coupe au mieux). • La complexité de la multiplication de deux chiffres de n mots est proportionnel à nlog2(3)n1.58

  10. Karatsuba: Exemple • A=171 659,B=891 512 • On prend D=1000A1=171, A0=659, B1=891,B0=512 • M2= A1.B1= 171891= 152 361, M0= A0.B0= 659512= 337 408 • M1= (A0+A1).(B0+B1)-M0-M1= (8301 403)-M0-M1= 1 243 602-M0-M1= 674 721 • R=D2.M2+D.M1+M0=153 306 058 408

  11. Méthode de Toom-Cook • On coupe en trois parties • A=A2D2n+A1Dn+A0, B=B2D2n+B1Dn+B0 • M=A.B= M4D4n+M3D3n+M2D2n+M1Dn+M0 • On pose P1=(A2+A1+A0).(B2+B1+B0), P-1=(A2+(-1)nA1+A0).(B2+(-1)nB1+B0), P2=(22nA2+2nA1+A0).(22nB2+2nB1+B0) • On a M0=A0B0, M1=2A2B2-(1/2)A0B0-(1/3)P-1+P1-(1/6)P2, M2=-A2B2-A0B0+(1/2)P-1+(1/2)P1et M3= - 2A2B2+(1/2)A0B0 -(1/6)P1-(1/2)P1+(1/6)P2, M4=A2B2 • On a 5 multiplications seulement. • La complexité est proportionnel à nlog3(5)n1.46

  12. Toom-Cook: Exemple • A=171 659, B=891 512 • D=100, A2=17, A1=16, A0=59, B2=89, B1=15, B0=12 • P1=(A2+A1+A0).(B2+B1+B0)=10 672 • P-1=(A2-A1+A0).(B2-B1+B0)=5 160 • P2=(4A2+2A1+A0).(4B2+2B1+B0)=63 282 • M0=A0B0=708, M4=A2B2=1 513, M1=2M4-(1/2)M0-(1/3)P-1+P1-(1/6)P2=1 077, M2=-M4-M0+(1/2)P-1+(1/2)P1=5 695, M3=-2M4+(1/2)M0-(1/6)P-1-(1/2)P1+(1/6)P2=1 679 • R=D4.M4+D3.M3+D2M2+D.M1+M0=153 036 058 408

  13. Généralisation • On peut continuer et couper en r=4,5,6… parties • La complexité est proportionnelle à nlog(2r-1)/log(r) • Théoriquement on peut donc obtenir des algorithmes de multiplications en n(1+) avec >0 aussi petit que possible • Mais les calculs intermédiaires sont de plus en plus lourds • Pratiquement on n’implémente même pas l’algorithme de Toom-Cook: si besoin, on lui préfère les algorithmes à transformées de Fourier

  14. Transformée de Fourier • Multiplier des nombres ou des polynômes, c’est pareil. • A=(an-1an-2 … a1a0) et B=(bn-1bn-2 … b1b0) • A(X)= an-1Xn-1+an-2Xn-2+ … +a1X1+a0,, B(X)= bn-1Xn-1+bn-2Xn-2+ … +b1X1+b0 • SiP(X)=A(X).B(X), alors A.B=P(D) • P est de degré 2n-2, si on connaît 2n valeurs distinctes de P, on peut déterminer P (par exemple par interpolation de Lagrange) • La transformée de Fourier se sert des points wk=exp(2ik/2n), 0k<2n • On définit Aj=A(wj), Bj=B(wj) et Pj=P(wj) • On a Pj=Aj.Bj

  15. Transformée de Fourier • On calcule les Aj et Bj par transformée de Fourier • On calcule les Pj = Aj . Bj • On fait une transformée de Fourier inverse pour déterminer les Pi • Le calcul des transformées de Fourier et de Fourier inverse se fait en une complexité n.log(n) • Le calcul des Pj est linéaire • Le coût global est donc en n.log(n)

  16. Multiplication: Complexité • L’algorithme scolaire a une complexité en n2 • L’algorithme de Karatsuba est en n1.58 • L’algorithme de Toom-Cook est en n1.46 • L’algorithme par transformée de Fourier de Schonage a une complexité en n.log(n) • Mais en fait, Winograd a démontré que la complexité de la multiplication est juste linéaire

  17. Equation f(x)=0 • On suppose f fonction réelle, de classe C1 sur un intervalle ]a,b[ qui s’annule dans l’intervalle • On définit par récurrence x0 ]a,b[, et xn+1=xn-f(xn)/f’(xn) pour n1 • Suivant f et x0, la suite xn converge vers une limite l telle que f(l)=0 • Méthode de Newton • Convergence généralement quadratique (le nombre de décimales exactes double à chaque itération)

  18. Calcul de l’inverse • Soit a un nombre réel positif, on cherche à calculer 1/a • On pose f(x)=1/x-a • L’itération devient xn+1=xn(2-a.xn) • Si xn-1/a= alors xn+1-1/a=a.2 • Si <1/a, alors a.2 <  • Donc si x0]0,2/a[, la suite converge vers 1/a • On a donc besoin d’une estimation grossière de 1/a pour initialiser l’algorithme • La convergence est quadratique

  19. Inverse: Exemple • a=0.52369817 • En regardant les premiers bits de a, on se rend compte que 0.5 < a < 1 donc 1 < 1/a < 2 on peut initialiser avec x0=1.5 • x1=x0(2-a.x0)=1.82167911 • x2=x1(2-a.x1)=1.90545810 • x3=x2(2-a.x2)=1.90948830 • x4=x3(2-a.x3)=1.90949683

  20. Division entière • Soit A et B deux entiers, on veut calculer A/B • On calcule q=1/B avec une précision inférieure à 1/A par excès (resp. par défaut) • On calcule A.q et on l’arrondit à l’entier par défaut (resp. par excès)

  21. Division entière (Exemple) • A=309 641 782 117, B=143 557 • En regardant les premiers bits de A, on sait que A<51011, on a donc à calculer q avec une précision de 210-12 • On a q=0.000 006 965 874 à 210-12 près • La valeur de q est par excès (l’itération suivante donne une valeur légèrement inférieure) • A.q=2 156 925.639 à moins de 1 près • On arrondit à l’entier inférieur et on trouve A/B=2 156 925

  22. Calcul d’un modulo • Soient A et B deux entiers, on veut calculer r=A mod B • On calcule le quotient entier q=A/B par l’algorithme précédent • On calcule alors r=A-q.B • Si B est petit par rapport à A, il existe d’autres méthodes

  23. Racine carrée • Soit a un nombre réel positif et soit à calculer a • On pose f(x)=x2-a et on utilise la méthode de Newton • xn+1=(xn+a/xn)/2 • Connue depuis l’antiquité sous le nom de méthode de Hénon • Problème: nécessite une division

  24. Racine carrée: 2ème méthode • f(x)=1/x2-a • L’itération devient xn+1=(3xn-a.xn3)/2 • La suite converge vers 1/a, on a plus qu’à multiplier par a pour trouver a. • Plus de divisions • Si x0]0, 2/a[, la suite converge quadratiquement vers 1/a

  25. Racine carrée: Exemple • a=5, 4<a<16, 2<a<4, 1/4<1/a <1/2 x0=0.375 • x1=(3x0-a.x03)/2=0.43066406 • x2=(3x1-a.x13)/2=0.44630628 • x3=(3x2-a.x23)/2=0.44721083 • x4=(3x3-a.x33)/2=0.44721359 • x5=(3x4-a.x43)/2=0.44721360 • x6=(3x5-a.x53)/2=0.44721360 (x6=x5on s’arrête) • a5.x6=2.23606800, la valeur est exacte à 5.10-8 près • En fait a= 2.23606797…

  26. Itérations d’ordre supérieur • On cherche à résoudre f(X)=0 • On suppose que x est proche d’une racine de l’équation • f(x+h)=f(x)+h.f’(x)+(h2/2).f’’(x)+ … (Taylor) • On résoud en h l’équation f(x)+h.f’(x)+ (h2/2).f’’(x)=0 • On trouve h=-(f’(x)/f’’(x))(1-(1-g(x))1/2) avec g(x)=2f(x)f’’(x)/f’(x)2 • f(x) est supposé petit et si f est une fonction suffisamment régulière, g(x) est aussi petit

  27. Itérations d’ordre supérieur (suite) • Si  est petit alors 1-(1-)1/2 /2+2/8 (Taylor) • On applique a g(x) et on trouve (1-(1-g(x))1/2)(g(x)/2).(1+g(x)/4) • D’où h -(f(x)/f’(x)).(1+f(x)f’’(x)/(2.f’(x)2)) • On pose xn+1=xn-(f(xn)/f’(xn)).(1+f(xn)f’’(xn) /(2.f’(xn)2)) • La suite converge cubiquement vers la racine de f(X)=0 • Itération de Householder (1970)

  28. Itérations d’ordre supérieur (fin) • On peut continuer et définir des itérations quartiques, quintiques etc … • Mais toutes les fonctions ne s’y prêtent pas • Il faut que les dérivés successives soient faciles à calculer • Il faut que le surcoût en calcul soit rentable en termes de vitesse de convergence

  29. Exemple: Inverse • Calcul de l’inverse de a • f(x)=1/x-a • On choisit la valeur initiale x0 proche de 1/a • On pose hn=(1-a.xn) • Itération quadratique: xn+1=xn+xnhn (Newton) • Itération cubique: xn+1=xn+xn(hn+hn2) (Householder) • Itération quartique: xn+1=xn+xn(hn+hn2+hn3) • Itération quintique: xn+1=xn+xn(hn+hn2+hn3+hn4) • Etc …

  30. 1.905458103 1.821679118 1.890664087 1.905458103 1.909495007 1.909496839 1.909488296 1.909496839 1.909496839 1.909496839 1.909496839 1.909496839 Exemple: Inverse • a=0.52369817, x0=1.5 quadratique cubique quartique n=1 n=2 n=3 n=4

  31. Calcul des fonctions élémentaires • Calcul des fonctions sin, cos, log, … • On se sert des propriétés de la fonction pour se ramener à un intervalle donné • Pour une précision limitée (une dizaine de chiffres décimaux), l’algorithme CORDIC peut faire le calcul avec uniquement des additions (c’est l’algorithme des calculatrices) • Pour une plus grande précision, on peut utiliser des développement de Taylor en se plaçant au milieu de l’intervalle pour accélérer la convergence • Pour une précision maximale, on utilise des techniques comme le binary splitting ou l’AGM

  32. Exemple: sinus • On a sin(x+2)=sin(x), sin(x+)=-sin(x), sin(-x)=-sin(x), sin(x)=cos(/2-x) • On peut donc se limiter à calculer dans l’intervalle [-/4, /4] • sin(x)=x-x3/6+x5/120+…=(-1)nx2n+1/(2n+1)! • cos(x)=1-x2/2+x4/24+…=(-1)nx2n/(2n)! • Les développements tronqués se calcule par l’algorithme de Horner pour limiter les calculs • Exemple: 1-x2/2+x4/24-x6/720=1-(x2/1.2).(1-(x2/3.4)(1-x2/5.6))

  33. Exemple numérique • Calcul de sin(8) à 5 décimales près • sin(8)=sin(3-8)=cos(8-5/2) • On a 8-5/26<10-5, on peut donc se limiter aux trois premiers termes du développement de cos(x) • On utilise donc cos(x)1-(x2/1.2).(1-(x2/3.4)) • Le calcul par l’algorithme de Horner est stable (les erreurs d’arrondis ne se propage pas) • On trouve numérique sin(8)=0.98935826008034 • En fait sin(8)=0.98935824662338

  34. 0.23876388301 0.066666666667 0.10991119991 0.076190476191 0.12332346666 0.30969097075 0.082783882784 0.14041543333 0.43656365692 0.090231990232 0.71828182846 0.16291649048 0.099111999112 0.19381941508 1.7182818285 Exemple: Calcul de e à 10 décimales • On a e=1+1/2! + 1/3! + … • On a 15! > 1010, on peut donc s’arrêter au 15 ème terme • On pose x15=1/15 et xn=(1+xn-1)/n et on a exp(1)1+x0 • En fait e2.718281828459 • La valeur trouvée est correcte à 10 décimales près

  35. Pour aller plus loin … • Site web: http://numbers.computation.free.fr • Arithmétique en précision arbitraire, Paul Zimmermann, Publication interne INRIA no 4272, septembre 2001 • AGM: Fast multiple-precision evaluation of elementary functions, Richard P. Brent, Journal of the ACM, vol 2 no 23 pages 242-251 (1976) • Calcul de : Pi and the AGM, J.M. et P.B. Borwein sur le site http://www.cecm.sfu.ca

More Related