- lnstitut Fourier Université de Grenoble | gR LE CALCUL NuMERlauE DE LA C_STANTE DIEULER . pa_ Jean-Pierre DEMAILLY - o. UN PEU DIHISTOIRE _ consLanEe d'L¬uler V = ]irrl l + ± ... ± - Iogn , diLc + + n_+m 2 fl parfois constante de _IascheIoni, a fait l'objet d'ulI gralI_ nornbl_e de travaux et d'é\Taluations numériques de_uis le _Ie siècle. Elle garde néanmoins Encore beaucoup de son mystère aujourd'hui. Ainsi, on ne sait toujour_ pas si v est ou non _rratioMel, EII dépit de plusieurs tentatives de dé_on_tration, dont celle de P. Appe_l l_] en 1926 qui avorta par suite d'une el_reur matériclle, et celle de h. Froda rl_l en 1965 qui reUo_e 5'ur un cl¬itère d'irrationnalité elI_re il__ complèternent démontl_é. Signalol_s t_utef_is que certa__s résultats de tra_scendance pour de5 e_r_ssion5' frdi_ant _Itervenir V ont ét_ ob- tenus par K. mahler El8] . Nous nou__ intéresserons ici wrtout à la mise el_ oeuvL`e d'algo- rithmes rapideoIent Lonverg_nts qui, o_tre l'intérêt calculat_ire: laisr sent espérer l'obtention de r¸sultats aritlI_oétiques. I_a prel_ière évaIuation de V est due naturelle.ment à Ieonhard EuIer, qui obtint la valeur o. 577218 en 1735 [12] , _ientôt étendue _ar Mascheroni et Uuelques autres. En 1781, Euler déterlnina la valeur _lus précise o. 577L156649rJ1532 [13] . Il fut suivi notammcnt par C. _.. G_LISS, avec 22 décirnales e_actes, puis par un certá_] num_re d_ rndth_n_aticiens a_IRlais du _Xe si_c]e. tc lecLeur pourra col]sulter J. N'. L. Gl'disher l15] pou3_ ]'hisLnriUue d_taill_ d_s calcl_]s ant___ieur_ á I87O. N.. Sh_nh6 l]9] publird ]lo d_cim_]cs, do_t. lo] ex_c__s, el] 186T-_] ; _cu aUr_b¬, l_ L`_li'_re mdth_m'dLicicn-'dst__>IIol_I_ lngtllu_ Fourlar - B. P. 74 - 38402 Sgl__-Merfln-d.Hbr_a C_d_x L_borB_olr8 da M6_hémg_lque_ Purea _aeotlé _u CNRS - Tél . l_6) 5_-46-00 - - - - - - J.C. AdaIos [ll calcula laborieusement 263 déc_ales, record qui devait tenir depuis sa publication en 1878 jusqu'à l'apparition des pre- miers ordinateurs et les 328 décimales de J.W. Wrench Jr [23] en 1952. Tous ces calc_ls, ainsi que celui ultérieur de D.E. Muth [17] en 19_2 avec 12T1 D , rel>osaient sur le d_velo_l>elI_ent asyn__totique de l +...+ ± - _ogn par la forlnule d'_uler - Mac L3urin. k teml>s n dc calcul prohi_itif des nombres de Bernoulli requis dans cette rn_thode conduisit mra w'. SN7eeney l21] à introduire un nouvel algorithrne plus F. , ' = - '_ -xdx . Syeeney obtint efflcace base sur la formule Y _ogx e o ainsi 3566 D en 1963, et sa méthode fut reprise successivement par Beyer-Waterman [3] , l4] en 19T4 ( T_l4 D , dont 48T9 exactes) et par R.P. Brent [7l , l8] T 20700 D en lg77T. Enfin en 1980, R.P. Brent et E. mc Millan llo] découvrirent un nouvel algorithme plus performant, utilisant les fonctions _e Bcs_el, et calculèrent 301_0 D l9] . Voici un bref aperçu des algoritrmes évoUués plus haLlt, avec analyse corn_arée des temps de calcul. |. FORMULE DIEULER - MAC LAURIN - Cette formule sera utilisée sous la forme suivaote (cf. _urbaki [5] ) : " . = f" ± (f(n) +f(l)) E f(LT f(x)dx + i±l l ' k _ [fl'j-'),n, f(2j-l),,,] + Rk + .', (2j)| - J± où b. est la suite des nombres de Bernoulli, définie par ' '_ bj zj ±= 'o_ 7 '-l j= c d_ sorte que bO = l , b ± - _ , b = _ , b, ± O . . . . l 2 - _e rerte Rk ert doMé par R ___ "k+"(x)dx k - (2k+1Tl B2k+l((x))' l . l ` ± _ m bxm-j ert le rn-iè_e polynôme de Bernoulli et ou B (x) , o j j m .± (h'F l partie fractionnaire de x . si nous posons f(x) = ± , la for- x mule ci-dessus e_tran]e d'aUrès D. Muth [17] : , + l I _ _ '...' ñ - b2 , , b,k , , n B,k+,((x]T = Ioen +± ± +_ _- +...+- --- - l dx. + 2 2n 2 n2 2k n2k , x2k+'L áua_ld n_ +_ , il vient l b2 b2k +_ B2k+l(lxl) V = - +- +...+- - r dk- . 2 2 2k L, x2k+2 Par soustra_tion, nous obtenons donc : v l + ' ' Iogn l b2 b2k '_ B2k+l'_x)T ± _ _...+ñ - - _ ' _ '...+ ,_2k - l 2k+2 d% . n x Ce dévelo_pement asp_totique diverge quand k _ +_ , mais il doMe ce_endant de _nnes valeurs app_ochées de v lorsque n et k sont bien chDisis. En effet, l'i_entité classique B2k+IT_x)T = '(-')k-' . _ 'i"'rnx (2k+1) | r±l ,2n,2k+1 entraAe _B2k+IT{x)T | s 4 (2k+1) ! 2k+1 (2nT et grace à la formule de Rirling, nous en déduisons |'_ B2k+l"x"dx 4 n( k k l 2k+2 ' - - - . n x n n nTTe _e rerte est donc très petit tant que k demeure sensiblement rtérieur à n . An_l se du ten_ s d_ calcul. D_slSIons p_r E cet zlgorlLrme et par E,(d) ]c tcm_s de _ _ ' -d . On aLLri_uc par conv_ntiolI une calcul dc v _ ]a prLclsbn IO unité de temps à chaque opération arithalétique élémentaire _ortant sur l mot de la machine). _a somme de 2 nornbres ayant d déc_nalYs nécessite donc, à des conrtantes près que nous négligerons ici, d uni- tés de tealps ; idem pour les produits ou quotients de nornbres en mul- tipréci,sion par des "petits'' nombres (l o_ot-n_achine). L'évaluation la plus brutale _e l + ± ... ± ' + + reclan_e alors d] 2 n unit_s dc tem_s. n_ns la _ratiUue, on chaisil_'d pour n une Uuiss_l]cc d_ 2 ou 10 , et _ogn ser_ évalué elI d' uniL_s de Lemps à partir des s_ries entières Iog2 ± 2Argth_ , IDg_ ± 2Argth_ B e_onen tielleIoent convergentes. , ± _ peuvent être calculés au D'autre part les non_bres B,k ,k n mo¨en de la formule de récurrence des noalbres de Bernoulli : B l k-_ k-' 2k+1 a2j 2k = 2k+1 _ -j_, 2j L(k-j) n ± n L'éta_e de récurrence exige environ kd unités de tem_s, à condition de stocker en mémoire les résultats partiels 'k?' B'j . , soit au 2J 2(k-J) total k' _ . " d pour l'evaluatlon de la somrne b2 b2k On aboutit donc à . _ +...'_ . . E (d) _ nd +d' +k'd . l Bien entendu, n et k doivent être choisis de manière que le terme d'erreur soit < 10-d . . , ce qul IInpose kIog_ _ dIDglO k . IP choix optimal est obtenu pour n _ k' , d'où kIDgk _ d IDglO , k d_ogIO h togd e' E (d) _ Cd'(IDgd)-' . I SiSnalons qu'l] e%lsLe une forn_ule (p_u colInuY) permctL_nt d_ conLourner le calcul dEs b,k . Cet_< forrnul_ e_Urirnc lc r*st< sous forn_c dc s_rie double rapidement convergente : +_ +_ e ! ± ± ... _ - `°g" + ± ' _ _l 2e+'2kn 2kn+1)...(2kn+e) . Y l + + + 2 o l 2n k l e , pour vérifier cette identité, on part de l'_ltégrale convergente ± '_-'±-± dx où p q >o I(P, q) ._ q ,-x _ _ . O l-x Pour tout entier n z l , un calcul ais_ domle T= ' - "-',dx - d IDg_ I(n,n) Tl +h +...+x O l-x ' de sorte que I(n,n) ± l +± ... ± - _ogn _V Uuand n _ +m . + + 2 n-l IP changement de variable x ± t' dans I(p,q) fournit d'autre part l' ident ité I(P,q) + I(pr,rT = I(pr,qr) . _ppliquons cerLe formule par récurrence avec p ± q , r = 2 . Il vient I(n,nT + I(2n,2) +...+ I(2kn, 2) = I(2kn, 2kn) , d'où à la limite quand k _ +_ : = , É_ I(2kn, 2) v I(n n) + h-l l l _ l 2kn-1 dx = l +- +...+-.- _ogn + r X - . 2 n-l k±l .J O l+x La formule aMoncée s'en déduit maintenant par développement en série de l/(l+xT : )-l +_ e l - l l-x _ ,, (i-x) _ - _ '-_ - e±o ,e+l . En particulier, pour n = l , on obtient la formule simple (E,) o ± _ + _ _ e! k±l _±l 2e"2k,2k+,, ,,k . . . +e) tP term_ g_néral de la série cs, majoré ,ar min,,-e-k-l,,e,-k,e, l_ précision 10-d est donc a_Leil]Le en somn_ónt póur h,e s d IDSIO/IDg2 e s d Lo 10 k LDg'L - _D6e ce qui laisse environ '°g'°d _ogd termes à calculer. _e temps de Iog2 calcul requis par (E,) admet donc l'estimation E,(d) _ Cd'IDgd , asymptotiquement _lférieure à celle de l'algoritrme (E,) . Nous allons voir qu'il existe en fait des algorithmes simples en oTd') , mais cela n'em_&che p_s l'algoritmle IE ) d'&tre salIs doute le plus effic3ce _our l les calculs manuels (d± loo) . z. ALGORITHME DE s EENEY IP point de départ en ert l'égalité L- +_ -'d, - V r'(l) - _ogt e - - , o Uui découle de l identité des définitiDns d'Euler et de Weierstrass de la fonction r : _ ± -l _'_ x -'d, e-yxll " l+x r(l+x) = t e ± e - . ` O n±l " Gr&ce à une iotégration par parties, on obtient v = F(x) - _ogx - R(x) avec x -t _ n-l n F(x) ±T _d, ± E (-') x O t n±l nn! ' -t -x , , , ,,h, k+l +_ e-'dt R(x) ± f+_ ±dt = ± l-j +...+ - . + (-l) (k+l)!F . . x t x x k k+2 x x t _a méthode (s,T la plus simple, utilisée par Sweeney l21] et Beyer- Waterman [3], consiste à choisir un entier x assez grand de manière que R(x) soit négligeable et à calculer la valeur appIoch__ V _ F(x) - IDgx . Comme O < R(x) < _ , ' . . -d esL la preclslon 10 x obtenue pour x_ d IDgIO . - l Etant doMé p z o , soit a l'unique racine réelle >o de p l'équation a T_oga - l) ± p , p p de sorte que aO ± e± 2.718 , a _ 3.591 , a, _ 4.319 , a, _ 4.971 . l La forn_ule de qirling n,on,re _ue n e_ " - - _ - est de l'ordre de e px n ! n pour n ± a h . I_a série F(x) doit dunc &tre son_mi_ ici _usUu'à u n ± a X . l _e calcul de cbaque terale de la série, après factorisation suivant la règle de _.rner rlx) = ± ±-± ±- ± l x l l l 2 2 ... n-l _-ñ II ... ' réclame 3 opérations óritholétiques (l multiplication, l divison, l sous- traction). Une difficulté su_plémentaire se présente ici du fait qu'une partie des chiffres significatifs ert perdue par corn_ensation des termes de signes opposés. Ie terrne de valeur absolue maximale étant de l'ordre de _ x pour n =x , c ' 'x si x z l , e e ".-O _ +_ -2xchvdv e, , TT -2x KOTx) = f e _ _e o _a constante _'Euler peut donc s'écrire sous la forme SO(x) KO(x) V=_-'Dgx-_' KO(h') -4k avEc O < < TTe si xzl . I (_) ° KO(x) Dans l'algoritbme (B) utilisé par Brent, le rerte _ est _ure- ment et sialplement négligé ; la précision l o-d ° . ert donc attelnte pour ± _ d IDglO , et les séries IO(x) , SO(x) doivent étre somm_es x jusUu'g n = a,h . Ie c_lcul de 2 2 2 2 I,(_) = I + _ _ ... _ l + x ... ... ,2 ' ' r ,,L n2 2 (n+I) n_c_ssiLe 2 _p_r'd_ions ariLrm_Liuucs pour ch_Uue tcrmc, et cclui dc - - - -- - - - -- 'O(x) x2 2 , , x2 _ ± H -- ...± - ... - -(...) ... + + + N 2 n2 n+l N 2 l ln+l) en exige 4 . IP temps requis par l'algoritrme TB) ert donc B(d) = a, x _d _oglO x 6 x d _ 12.4d' . KO(x) Como_e dans le cas dc la méthode de S_Te__ey, le rerte _ U_ut étre évalué au moyen d'un déveloPl>ement asym_totiUue. On 'd en cffeL IO(x)KO(h) ± _ _ exU(2_(cosu - chv))du dv | | nO l +_ 2n de = __ e_(-4x r)drl ie, . O O ll-re grâce au changement de variable reie = . ' _ . A développement sln , en série 2n '_ (2k)!' 2k Osr_ l l d6 = E _r l . ' _ fo ,,-rei8| h_o 2 . on déduit alors le dévelo_Uement asym_totique (divergent) IO(x)I_,(x, ~ _ L T2k)!' x k!4,,6x,2k _e terme général dE ce développement passe par un minirnuln poul_ k = 2x , et o_ peut véri_.ier que le rerte corresp_ndant ert de l'ordre de grandeur du terme k = 2x , soit e-4x . - ' 8rr x"' `a `a'e`' a_pr°chee KO (x) l _ T2k)!' _ _ 4x,0,x,2 k=o ,k!,4,,6x,2k ert donc affectée d'une erreur ne dépass_nt Uós e-8x . Cette méthode arnène à chaisir x ± _d LoglO et conduit au temps de calcul B'(d) ± ± ± IDglOd' _ 9._ d' , 4 a + 3 2 optimal _aro_i tous ]es _lgoritm_es _résentés ici. Nous terIninons en doM_nt lc "hit-p_l__d_'' de c__ al6oritm_es, ran6_s U_r ord_¬e d'effic_cité crois__nLc. Wgorithmes 'u'e' '__'eeney Brent : E E s s2 si si s, sj B B_ l 2 l Te_ps de cd c _ogd 4g 6 37 6 26 T 26 o 22.7 16.9 12.4 9._ calcul /d' _ . . . . (_ogd) Si&]alons qu'il eh'iste des algorithn_es théori_uemen_ encore _lus ra_idcs, pern_etLal]L d'é_Taluer _7 à ]o-d _>r_s el] O(d(logd)' _ogIDgd) uni__s de ternps . Ces derniers reposent sul¬ l'utilis_tion d_ l'alRorithme dc mul- tiplicótion rapide de Schó.nhage-_rassen l20l et sur une factorisation par blocs de la série F(x) (resp. ex , IO(x) , SO(x)) ; voir Brent [6] . _ tels algorithmes "rapides" sont toutefois trés difficiles à pI_grammer et ne l'emportent sur les aIgoritrmes "classiques" présentés ici que lorsque d est très grand. 4. DEVELOPPEMENT EN FRACTION CONTINUE DE Y -. Ies valeurs numériUues de Y obtenues par les auteurs mention- nés ci-dessus ont été utilisées pour déterminer les développen_ents en fraction continue de v et e' , qui sont COMUS maintenant jusqu'à l'o rdl'e 29 000 (Brent - _lac _'_illan llllT. La distributiol] rtatistique des réduites successives ne fait ap_aranl_e aucune diUérence significa- tive au niveau de 5 % par rapport à la loi de Gauss - _smin (cf. #intchine [16]) : f ± log, l+± - log, l+± n n n+l doMant la frÉquenc_ des réduites égales à n da__s le développen_ellt de presque tout non_bre réel. On obtient d'auLre part le résultat sui- vant, qui rend e_rén_en_ent lrrl_robable l_ rationdlité de Y ou e' : Th_or_me. si v ou e` .- p/d Uour d_s entiers p,ö posi- tifs, alors Q _ IOl5OOO . | BIBLIOGRAPHIE rl] J.C. ADAMS. - On the value of Euler's conrtant ; Proc. Roy. _c. _ondon, v.27, 1878, p. 88-9_. Voir aussi v.42, 1887, p. 22-25. l2] P. APPELL. - Rr la nature aritrmétique de la constante d'Euler ; Com_tes Rend. Ac. k. Paris, v.l82, 12 avril IJ26, p. 897-899 ; 19 avril 1926, p. 9_9. l3] U7.A. BL¬YEl_ g m.s. _'ATEI_MAN. - ErIor analysis _f a con_l>u- tation of Euler's consLalIt ; mdth. of L'ornt>. , v. 28, 19T4, p. 599-604. l4] W.A. BEYER g M.S. _iATERMAN. - _ciInals and partial quotients of Euler's conrtant and en 2 ; UMT 19, Math. of Comp. , v.28, 1974, p. 66T. Errata : Math. of Comp. , MTE 549, v.32, 19T8, p. 317-318. [5] N. Bou_BAa. - Fonctions d'une variable réelle ; chap. N, é l, n°7 ; Nermann, Paris, 1951. l6] R.P. BRENT. - The complexity of multi_le-precision aritrmetic ; Co3nplexity of Computational p_oblem _lving (R. s. Anderssen and R.P. Brent, Editors), Univ. of Queensland Press, Brisbane, 19T6, p. 126-165. l7] R.P. BRENT. - Computation of the regular continued fraction for Euler's constant ; _lath. of Comp. , v.3l, July 1977, _. T7l-77T. [8] R.P. BRENT. - v and e_(Y) to 20700 D and their regular continued fractions to 20000 partial quotients ; UMT l, __ath. of Comp. , v.32, 1978, p. 311. l9] R.P. B_ENT. - Euler's conrtant and its eqonential to 30,100 decimals ; _Iath. of Comp. , UMT File. llo] R.P. BRENT g E.M. Mc _IILLAN. - _me neNT algorithms for bigh-precision computation of Euler's constant i _lath. of ComU. , v.34, JalIuary 1980, p. 305-312. [ll] _.P. BWENT g E.M. Mc NIILL,_N. - Thc firrt 29,000 partial quotients ln the regul_r cont_Iu_d fractiolI for Euler's cons- tant and its eqonel]tial i m'dth. of ComU. , umT _¬il_. - . [12] L. EULER. - _ progressionibus harmonicis observationes ; Euleri ppera Oalnia, _r. l, v.l4, Teubner, Ieipzig and Berlin, 1925, p. 93-100. [13] L. EULER. - De summis serierum numeIos Bernoullianos invol- ventium i Euleri Opera OrrLnia, _r. l, v.l5, Teubner, _eipzig and Berlin, 1927, p. 91-130. Voir en particulier p. 115. _P calcul ert donné en détail _. 569-583. l]4] A. FRODA. - La constante d'Euler ert irratioMelle ; Atti Accad. Naz. L_]cei R_nd. Cl. ki. _is. Ma_. N_tur. (8) 38 (19_5), p. 338-344. [15l J.W.L. G_AIMER. - Nirtory of Euler's conrtant ; Messenger of Math. , v. l, 1872, p. 25-30. [16] A. Ya. QINTcHINE (A. Ja. HIN_ ). - Continu rractians ; IN ed 3e éditian, traduction anglaise par P. _V'yM, NoordhoU, Groningen, 1963. l17l D.E. auTH. - Euler's constant to 1271 places ; _Iath. of Con_p. , v. 16, 1962, p. 275-281. [18] K. MAHIE_. - Ap_lications of a theore_ by A.B. Shidbvskii ; pIoc. Roy. _c. Iondon, A 305 (19_8), p. 149-173. [19l w. SHAh`KS. - On the numerical value of Euler's conrtant ; PIoc. Roy. bc. Iondon, v.l5, 1867, p. 429-432 ; v. 20, 1871, p. 29-34. [20] A. SCHÖ_HAGE g v. STRAUE_. - kbnelle Multiplikation gI_sser ZaNen ; Computing, v.7, 1971, p. 281-292. [21] D.W. S__'EEhEY. - On the coml>utation of Euler's conrtant ; _Iath. of Cornp. , v. 17, 19_3, p. ITO-178. [22] G.N. NiAT_N. - A Treetise on the Theory of Bessel nnctions ; 2"d ed_ion, Cambridge Univ. Press, IDndon, 1944. E23l J.W. U.RE_CH, Jr. - A neNT calculation of Euler's conkant ; _lath. Tables g other Aids to Com_. , v.6, 1952, _. 255. _r. l24] E. ARTIN. - The Gamma Function ; _blt, _inehart and Winrton, 1964, traduit de : Ertührung in die Theorie der Gamrnafun_dion, Teubner, 1931. [25] R. CAMPBELL. - Ies intégrales eulérieMes et leurs a_plications ; Anod, Paris, 196_. --- (juin 1984)