Partout, code compris, sans tenir compte des majuscules, accents, espaces ni symboles. Une formule est aussi calculée : \frac{1}{6}\pi^2, zeta(2) ou 1/6*pi^2 trouvent Bâle par sa valeur, 0,91596 trouve Catalan par ses décimales, et log_b(xy) reconnaît la loi L·1 même avec d'autres lettres. « log » seul se lit ln. / pour chercher, Échap pour effacer.
# mpmath, 25 chiffres significatifs from mpmath import mp, mpf, nsum, inf, ln mp.dps = 25 s = nsum(lambda k: mpf(-1)**k/(k+1), [0, inf]) print(s) # 0.6931471805599453094172321 print(ln(2)) # 0.6931471805599453094172321
# Leibniz 1674 — c'est aussi β(1), la fonction β de Dirichlet en 1 from mpmath import mp, mpf, nsum, inf, pi mp.dps = 25 s = nsum(lambda k: mpf(-1)**k/(2*k+1), [0, inf]) print(s) # 0.7853981633974483096156609 print(pi/4) # 0.7853981633974483096156609
# mélange de deux constantes : π et ln2 réunis from mpmath import mp, mpf, nsum, inf, pi, ln, sqrt mp.dps = 25 s = nsum(lambda k: mpf(-1)**k/(3*k+1), [0, inf]) rhs = (sqrt(3)*pi + 3*ln(2))/9 print(s) # 0.8356488482647210533371035 print(rhs) # 0.8356488482647210533371035
# la famille continue : chaque m livre sa perle from mpmath import mp, mpf, nsum, inf, pi, sqrt, ln mp.dps = 25 s = nsum(lambda k: mpf(-1)**k/(4*k+1), [0, inf]) rhs = (pi + 2*ln(1 + sqrt(2)))/(4*sqrt(2)) print(s) # 0.8669729873399110375739952 print(rhs) # 0.8669729873399110375739952
# Euler 1735 : ζ(2) from mpmath import mp, zeta, pi mp.dps = 25 print(zeta(2)) # 1.644934066848226436472415 print(pi**2/6) # 1.644934066848226436472415
# Euler encore : toutes les ζ(2k) sont des multiples rationnels de π²ᵏ from mpmath import mp, zeta, pi mp.dps = 25 print(zeta(4)) # 1.082323233711138191516004 print(pi**4/90) # 1.082323233711138191516004
# η(2) = (1 − 2¹⁻²)·ζ(2) : alterner divise Bâle par deux from mpmath import mp, mpf, nsum, inf, pi mp.dps = 25 s = nsum(lambda n: mpf(-1)**(n+1)/n**2, [1, inf]) print(s) # 0.8224670334241132182362076 print(pi**2/12) # 0.8224670334241132182362076
# développer −ln(1−x) en série et intégrer terme à terme redonne ζ(2) from mpmath import mp, quad, log, pi mp.dps = 25 I = quad(lambda x: -log(1-x)/x, [0, 1]) print(I) # 1.644934066848226436472415 print(pi**2/6) # 1.644934066848226436472415
retirer à \(\zeta(s)\) sa moitié paire \(\zeta(s)/2^{s}\) ne laisse que les impairs : \(\sum_k(2k+1)^{-s}=(1-2^{-s})\,\zeta(s)\) pour \(\operatorname{Re}s>1\). Avec \(s=2\), c'est \(\tfrac34\cdot\tfrac{\pi^{2}}{6}\) ; avec \(s=4\), \(\tfrac{15}{16}\cdot\tfrac{\pi^{4}}{90}=\tfrac{\pi^{4}}{96}\). Contraste saisissant : les mêmes impairs avec des signes alternés donnent \(\sum(-1)^{k}/(2k+1)^{2}=G\), la constante de Catalan (III·1), dont on ne connaît aucune forme close.
# s = 2 à 25 chiffres, s = 4, puis la règle des impairs pour 50 valeurs de s from mpmath import mp, mpf, nsum, inf, pi, zeta, nstr import random mp.dps = 25 random.seed(1) print(nsum(lambda k: 1/(2*k+1)**2, [0, inf]))# 1.233700550136169827354311 print(pi**2/8) # 1.233700550136169827354311 print(nstr(nsum(lambda k: 1/(2*k+1)**4, [0, inf]) - pi**4/96, 3))# -2.58e-26 # pour un s quelconque, par un autre chemin : Σ (2k+1)^-s = 2^-s · ζ(s, ½) (Hurwitz) print(nstr(max(abs(2**-x * zeta(x, mpf(1)/2) - (1 - 2**-x)*zeta(x))# 2.58e-26 for x in (mpf(random.uniform(1.2, 12)) for _ in range(50))), 3))
la même valeur que II·1, mais construite uniquement avec les nombres premiers : c'est le produit d'Euler (L·47) en \(s=2\). Son inverse, \(\prod_p(1-p^{-2})=6/\pi^{2}\approx0{,}608\), est la probabilité que deux entiers pris au hasard soient premiers entre eux : pour chaque premier \(p\), les deux ne doivent pas être tous deux multiples de \(p\), ce qui arrive avec probabilité \(1-1/p^{2}\) (voir VIII·1 et VIII·2). Le produit converge lentement ; pour les 25 chiffres, le code passe par la fonction zêta des premiers, \(\ln\prod_p=\sum_{k\geq1}P(2k)/k\) avec \(P(s)=\sum_p p^{-s}\).
# produits partiels jusqu'à 10², 10⁴, 10⁶ ; puis le produit complet, par la zêta des premiers N = 10**6 crible = bytearray([1])*(N+1); crible[0] = crible[1] = 0 for q in range(2, int(N**0.5)+1): if crible[q]: crible[q*q::q] = bytearray(len(crible[q*q::q])) premiers = [q for q in range(N+1) if crible[q]] from mpmath import mp, mpf, fprod, exp, nsum, inf, primezeta, pi, nstr mp.dps = 25 partiel = lambda x: fprod(1/(1 - mpf(q)**-2) for q in premiers if q <= x) print(nstr(partiel(100), 12), nstr(partiel(10**4), 12), nstr(partiel(10**6), 12))# 1.64194519662 1.64491792075 1.64493395536 print(exp(nsum(lambda k: primezeta(2*k)/k, [1, inf]))) # tous les premiers# 1.644934066848226436472415 print(pi**2/6) # 1.644934066848226436472415
# β(2) — on ignore encore si G est irrationnel from mpmath import mp, mpf, nsum, inf, catalan mp.dps = 25 s = nsum(lambda k: mpf(-1)**k/(2*k+1)**2, [0, inf]) print(s) # 0.9159655941772190150546035 print(catalan) # 0.9159655941772190150546035
# β impair : forme close en π ; β pair (comme G) : aucune connue from mpmath import mp, mpf, nsum, inf, pi mp.dps = 30 # 5 chiffres de garde s = nsum(lambda k: mpf(-1)**k/(2*k+1)**3, [0, inf]) rhs = pi**3/32 mp.dps = 25 print(s) # 0.9689461462593693804836348 print(rhs) # 0.9689461462593693804836348
# arctan x = Σ(−1)ᵏx²ᵏ⁺¹/(2k+1) : on retombe sur III·1 from mpmath import mp, quad, atan, catalan mp.dps = 25 I = quad(lambda x: atan(x)/x, [0, 1]) print(I) # 0.9159655941772190150546035 print(catalan) # 0.9159655941772190150546035
# Wallis 1656 : le premier produit infini pour π from mpmath import mp, mpf, nprod, inf, pi mp.dps = 25 p = nprod(lambda n: (4*n**2)/(4*n**2 - 1), [1, inf]) print(p) # 1.570796326794896619231322 print(pi/2) # 1.570796326794896619231322
Eugène Catalan, 1873 : le même que les nombres de X·2 et la constante G de III·1. On prend les pairs sur les impairs, par blocs qui doublent de taille, avec des exposants qui se divisent par deux. Même allure que Wallis (IV·1), mais c'est \(e\) qui sort, pas \(\pi\). La convergence est lente, l'erreur n'étant divisée que par 2 à chaque bloc : le code calcule donc chaque bloc exactement grâce à \(\Gamma\), pour aller jusqu'au 90ᵉ.
# 5, 20 puis 90 blocs : la lenteur, puis les 25 chiffres from mpmath import mp, mpf, loggamma, log, exp, e, nstr mp.dps = 70 # chiffres de garde # log du produit a·(a+2)·…·b, calculé exactement avec Γ (comme L·17) lp = lambda a, b: ((b-a)//2 + 1)*log(2) + loggamma(mpf(b)/2 + 1) - loggamma(mpf(a)/2) def P(K): # les K premiers blocs s = log(2) for k in range(1, K+1): s += (lp(2**k+2, 2**(k+1)) - lp(2**k+1, 2**(k+1)-1)) / 2**k return exp(s) print(nstr(P(5), 25)) # 2.689109944174376163286841 print(nstr(P(20), 25)) # 2.718280930017321431871364 print(nstr(P(90), 25)) # 2.718281828459045235360287 print(nstr(e, 25)) # 2.718281828459045235360287
Nicholas Pippenger, 1980. Même architecture que IV·2, mais chaque impair \(m\) est encadré par ses deux voisins pairs, \((m-1)(m+1)/m^{2}\) : l'erreur est alors divisée par 4 à chaque bloc au lieu de 2. Les facteurs de Wallis (IV·1) sont exactement de la forme \(\frac{2n}{2n-1}\cdot\frac{2n}{2n+1}\) ; ici on inverse le rôle des pairs et des impairs, et les exposants \(1/2^{k}\) changent π en \(e\).
# 5, 20 puis 50 blocs from mpmath import mp, mpf, loggamma, log, exp, e, nstr mp.dps = 70 # chiffres de garde # log du produit a·(a+2)·…·b, calculé exactement avec Γ (comme L·17) lp = lambda a, b: ((b-a)//2 + 1)*log(2) + loggamma(mpf(b)/2 + 1) - loggamma(mpf(a)/2) def P(K): s = log(2)/2 for k in range(2, K+1): # impairs m entre 2^(k-1) et 2^k a, b = 2**(k-1)+1, 2**k-1 s += (lp(a-1, b-1) + lp(a+1, b+1) - 2*lp(a, b)) / 2**k return exp(s) print(nstr(P(5), 25)) # 1.359362096221469097802607 print(nstr(P(20), 25)) # 1.359140914229728639590219 print(nstr(P(50), 25)) # 1.359140914229522617680144 print(nstr(e/2, 25)) # 1.359140914229522617680144
# Viète 1593 : le plus ancien produit infini de l'histoire from mpmath import mp, mpf, sqrt, pi mp.dps = 30 # 5 chiffres de garde a, p = mpf(0), mpf(1) for _ in range(60): a = sqrt(2 + a) p *= a/2 rhs = 2/pi mp.dps = 25 print(p) # 0.6366197723675813430755351 print(rhs) # 0.6366197723675813430755351
avec \(n-1\) racines en tout. C'est le demi-périmètre du polygone régulier à \(2^{n}\) côtés inscrit dans un cercle de rayon 1, soit \(2^{n}\sin(\pi/2^{n})\) : \(2\sqrt2\) pour le carré, puis 3,0615 pour l'octogone, 3,1214, 3,1365… Archimède, en poussant jusqu'à 96 côtés (avec des polygones inscrits et circonscrits), encadrait déjà \(3+\tfrac{10}{71}<\pi<3+\tfrac17\). Sœur de Viète (V·1) : mêmes radicaux imbriqués, en produit chez Viète. Un piège de calcul : sur ordinateur, \(2-\sqrt{2+\cdots}\) soustrait deux nombres presque égaux, et la formule telle quelle se dégrade après une quinzaine de doublements, jusqu'à donner 0. La version stable divise au lieu de soustraire : le côté suivant vaut \(s/\sqrt{2+\sqrt{4-s^{2}}}\).
# la formule en double précision, sa ruine, la version stable ; puis 60 doublements à 60 chiffres import math from mpmath import mp, mpf, sqrt, pi, nstr def naif(n): # double précision, formule telle quelle r = 0.0 for _ in range(n-2): r = math.sqrt(2 + r) return 2**(n-1) * math.sqrt(2 - r) def stable(n): # côté des 2ⁿ-gones : s ← s / √(2 + √(4 − s²)) s = math.sqrt(2.0) # le carré, n = 2 for _ in range(n-2): s = s / math.sqrt(2 + math.sqrt(4 - s*s)) return 2**(n-1) * s print([round(naif(n), 10) for n in (2, 3, 4, 10, 15)])# [2.8284271247, 3.0614674589, 3.1214451523, 3.1415877253, 3.1415926548] print(naif(25), naif(28), naif(30)) # l'arrondi dévore tout# 3.142451272494134 3.4641016151377544 0.0 print(stable(30), stable(60)) # 3.141592653589795 3.141592653589795 mp.dps = 60 # ou alors beaucoup de chiffres de garde r = mpf(0) for _ in range(58): r = sqrt(2 + r) print(nstr(2**59 * sqrt(2 - r), 25))# 3.141592653589793238462643 print(nstr(pi, 25)) # 3.141592653589793238462643
# −ln(1−x) en x = 1/2 : converge bien plus vite que I·1 from mpmath import mp, mpf, nsum, inf, ln mp.dps = 25 s = nsum(lambda n: 1/(n*mpf(2)**n), [1, inf]) print(s) # 0.6931471805599453094172321 print(ln(2)) # 0.6931471805599453094172321
# Euler, via la réflexion Li₂(x) + Li₂(1−x) = ζ(2) − ln x·ln(1−x) from mpmath import mp, mpf, nsum, inf, pi, ln mp.dps = 25 s = nsum(lambda n: 1/(n**2*mpf(2)**n), [1, inf]) rhs = pi**2/12 - ln(2)**2/2 print(s) # 0.5822405264650125059026563 print(rhs) # 0.5822405264650125059026563
# le carré de l'intégrale se calcule en coordonnées polaires from mpmath import mp, quad, exp, sqrt, pi, inf mp.dps = 25 I = quad(lambda x: exp(-x**2), [0, inf]) print(I) # 0.8862269254527580136490837 print(sqrt(pi)/2) # 0.8862269254527580136490837
# Serret 1844 : x = tan θ, puis la symétrie θ ↦ π/4 − θ from mpmath import mp, quad, log, pi mp.dps = 25 I = quad(lambda x: log(1+x)/(1+x**2), [0, 1]) print(I) # 0.2721982612879502663125861 print(pi*log(2)/8)# 0.2721982612879502663125861
# pas de 25 chiffres ici : l'erreur relative décroît en O(log n / n). # φ(1)+…+φ(n) ~ 3n²/π², car deux entiers au hasard sont premiers # entre eux avec probabilité 6/π². On regarde donc la convergence. from math import sqrt, pi N = 10**6 phi = list(range(N+1)) for p in range(2, N+1): if phi[p] == p: # p est premier for m in range(p, N+1, p): phi[m] -= phi[m] // p S = 0 for n in range(1, N+1): S += phi[n] if n in (10**3, 10**4, 10**5, 10**6): print(n, n*sqrt(3)/sqrt(S)) print(pi) # 1000 3.140412759431513 # 10000 3.141534213179507 # 100000 3.1415847755808897 # 1000000 3.1415926460191628 (7 chiffres justes) # π 3.141592653589793
# la somme partielle converge lentement (≈ 9 chiffres à 10⁶) ; # l'égalité exacte vient de Σ μ(n)/nˢ = 1/ζ(s), vérifiée à 25 chiffres. # C'est aussi la densité des couples premiers entre eux (voir VIII·1). from math import fsum from mpmath import mp, zeta, pi N = 10**6 mu, prime = [1]*(N+1), [True]*(N+1) for p in range(2, N+1): if prime[p]: for m in range(2*p, N+1, p): prime[m] = False for m in range(p, N+1, p): mu[m] = -mu[m] for m in range(p*p, N+1, p*p): mu[m] = 0 print(fsum(mu[n]/n**2 for n in range(1, N+1)))# 0.6079271020404619 mp.dps = 25 print(1/zeta(2)) # 0.6079271018540266286632768 print(6/pi**2) # 0.6079271018540266286632768
le logarithme du PPCM est la seconde fonction de Chebyshev : \(\ln\operatorname{ppcm}(1,\dots,n)=\psi(n)=\sum_{m\leq n}\Lambda(m)\), car chaque premier \(p\) y figure avec l'exposant \(\lfloor\log_p n\rfloor\), un par puissance \(p^{k}\leq n\) (voir L·30). Par exemple \(\operatorname{ppcm}(1,\dots,10)=2^{3}\cdot3^{2}\cdot5\cdot7=2520\) et \(\psi(10)=\ln2520\). Dire que \(\psi(n)/n\to1\), c'est exactement le théorème des nombres premiers (1896) ; Chebyshev avait déjà encadré ce rapport entre 0,92 et 1,11 vers 1850. Sa sœur, \(\theta(x)=\ln(x\#)\) avec la primorielle \(x\#\) (produit des premiers jusqu'à \(x\)), tend vers 1 de la même façon. C'est aussi pourquoi les dénominateurs de \(H_n\) grandissent comme \(e^{n}\) (L·36).
# ppcm(1..10) ; ψ et θ jusqu'à un million ; puis ppcm(1..1000)^(1/1000) calculé exactement N = 10**6 crible = bytearray([1])*(N+1); crible[0] = crible[1] = 0 for q in range(2, int(N**0.5)+1): if crible[q]: crible[q*q::q] = bytearray(len(crible[q*q::q])) premiers = [q for q in range(N+1) if crible[q]] from math import lcm, log, exp from functools import reduce from mpmath import mp, mpf, e, nstr mp.dps = 25 print(reduce(lcm, range(1, 11)))# 2520 psi = lambda x: sum(log(q) * int(log(x)/log(q) + 1e-12) for q in premiers if q <= x) theta = lambda x: sum(log(q) for q in premiers if q <= x) print([round(exp(psi(x)/x), 5) for x in (10**2, 10**4, 10**6)])# [2.56114, 2.72193, 2.71716] print([round(theta(x)/x, 5) for x in (10**2, 10**4, 10**6)])# [0.83728, 0.9896, 0.99848] print(nstr(mpf(reduce(lcm, range(1, 1001)))**(mpf(1)/1000), 8), nstr(e, 8))# 2.7092746 2.7182818
car \(V_{2k}(1)=\pi^{k}/k!\) : la série de l'exponentielle fait le reste. \(e^{\pi}\) est la constante de Gelfond, dont on sait qu'elle est transcendante.
# mpmath, 25 chiffres significatifs from mpmath import mp, mpf, pi, gamma, nsum, inf, exp mp.dps = 25 V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n s = nsum(lambda k: V(2*k), [0, inf]) print(s) # 23.14069263277926900572909 print(exp(pi)) # 23.14069263277926900572909
c'est \(\Gamma(\tfrac12)\). Avec \(t=x^{2}\), on retombe sur deux fois l'intégrale de Gauss (VII·1). C'est par là que \(\sqrt\pi\) entre dans le volume des boules de dimension impaire.
# quadrature directe gênée par la singularité en 0 : on pose t = eᵘ, # ce qui donne une fonction lisse ; les queues au-delà de [−140, 5] # pèsent moins de 10⁻²⁵. from mpmath import mp, quad, exp, sqrt, pi mp.dps = 25 I = quad(lambda u: exp(u/2 - exp(u)), [-140, 0, 5]) print(I) # 1.772453850905516027298167 print(sqrt(pi)) # 1.772453850905516027298167
chaque terme est la part du cube unité remplie par sa boule inscrite, de rayon \(\tfrac12\) : 78,5 % du carré, 52,4 % du cube, 0,25 % en dimension 10. Comme \(V_{2k}(\tfrac12)=(\pi/4)^{k}/k!\), la somme vaut \(e^{\pi/4}\) ; plus généralement \(\sum_k V_{2k}(r)=e^{\pi r^{2}}\), et IX·1 est le cas \(r=1\).
# mpmath, 25 chiffres significatifs from mpmath import mp, mpf, pi, gamma, nsum, inf, exp mp.dps = 25 V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n s = nsum(lambda k: V(2*k, mpf(1)/2), [0, inf]) print(s) # 2.19328005073801545655977 print(exp(pi/4)) # 2.19328005073801545655977
c'est le cas \(n=3\) de la formule de Dobinski (1877), \(\sum_{k} k^{n}/k! = e\,B_n\), où \(B_n\) est le nombre de Bell : le nombre de façons de répartir \(n\) objets en groupes. Il y en a 5 pour 3 objets : {abc}, {ab|c}, {ac|b}, {bc|a}, {a|b|c}. Le code compte ces répartitions une par une, puis vérifie la formule pour \(n=1\) à 6.
# Bell compté par énumération, puis Dobinski à 25 chiffres def partitions(elems): # toutes les partitions d'un ensemble if not elems: yield [] return x, reste = elems[0], elems[1:] for p in partitions(reste): yield [[x]] + p # x seul dans un bloc for i in range(len(p)): # ou x rejoint un bloc yield p[:i] + [[x] + p[i]] + p[i+1:] from mpmath import mp, nsum, factorial, inf, e mp.dps = 25 B = [sum(1 for _ in partitions(list(range(n)))) for n in range(7)] print(B) # [1, 1, 2, 5, 15, 52, 203] print(all(abs(nsum(lambda k: k**n/factorial(k), [0, inf]) - e*B[n]) < 1e-20# True for n in range(1, 7))) print(nsum(lambda k: k**3/factorial(k), [0, inf]))# 13.59140914229522617680144 print(5*e) # 13.59140914229522617680144
où \(C_n=\binom{2n}{n}/(n+1)\) : 1, 1, 2, 5, 14, 42… Les nombres de Catalan comptent les parenthésages, les arbres binaires, les partitions non croisées, les chemins qui ne passent jamais sous l'horizon. Leurs inverses, sommés, font surgir \(\pi\) et \(\sqrt3\).
# mpmath, 25 chiffres significatifs from mpmath import mp, nsum, binomial, inf, pi, sqrt mp.dps = 25 C = lambda n: binomial(2*n, n)/(n+1) print([int(C(n)) for n in range(8)])# [1, 1, 2, 5, 14, 42, 132, 429] print(nsum(lambda n: 1/C(n), [0, inf]))# 2.806133050770763489152924 print(2 + 4*sqrt(3)*pi/27) # 2.806133050770763489152924
\(\binom{2n}{n}\) compte les chemins de \(n\) pas vers la droite et \(n\) vers le haut : 1, 2, 6, 20, 70… Sœur de X·2, avec le même \(\sqrt3\,\pi/27\), à un facteur 2 près : Catalan n'est qu'un binomial central divisé par \(n+1\).
# mpmath, 25 chiffres significatifs from mpmath import mp, nsum, binomial, inf, pi, sqrt, mpf mp.dps = 25 print(nsum(lambda n: 1/binomial(2*n, n), [0, inf]))# 1.736399858718715077909795 print(mpf(4)/3 + 2*sqrt(3)*pi/27)# 1.736399858718715077909795
où \(a_n\) est le nombre de Fubini : les classements de \(n\) concurrents, ex æquo permis (1, 1, 3, 13, 75, 541…). Il croît comme \(n!/\bigl(2(\ln2)^{n+1}\bigr)\), et la convergence est fulgurante : l'erreur est divisée par environ 9 à chaque cran. Comme Bell (X·1), il a sa formule à la Dobinski : \(a_n=\sum_{k} k^{n}/2^{k+1}\).
# valeurs exactes par récurrence, puis convergence du rapport from math import comb from mpmath import mp, mpf, ln, nsum, inf mp.dps = 25 a = [1] for n in range(1, 61): # on choisit d'abord les premiers ex æquo a.append(sum(comb(n, k)*a[n-k] for k in range(1, n+1))) print(a[:8]) # [1, 1, 3, 13, 75, 541, 4683, 47293] print(nsum(lambda k: k**4/mpf(2)**(k+1), [0, inf]))# 75.0 r = lambda n: mpf(n*a[n-1])/a[n] print(r(5)) # 0.6931608133086876155268022 print(r(10)) # 0.6931471804369557443633155 print(r(20)) # 0.6931471805599453093587796 print(r(60)) # 0.6931471805599453094172321 print(ln(2)) # 0.6931471805599453094172321
autrement dit \(n!\approx\sqrt{2\pi n}\,(n/e)^{n}\). De Moivre trouve la forme vers 1730 ; Stirling identifie la constante \(\sqrt{2\pi}\), grâce au produit de Wallis (IV·1, et L·27 pour le mécanisme). En version logarithmique, \(\ln n! = n\ln n - n + \tfrac12\ln(2\pi n) + \tfrac{1}{12n} - \tfrac{1}{360n^{3}} + \dots\) : les corrections viennent de la formule d'Euler–Maclaurin et sont faites des nombres de Bernoulli de L·25. C'est d'elle que sortent le pic en dimension 5 et la formule d'A121546.
# convergence lente en 1/(12n) ; corrigée par Euler–Maclaurin, elle devient fulgurante from mpmath import mp, mpf, loggamma, log, exp, sqrt, pi mp.dps = 25 r = lambda n: exp(loggamma(n+1) - log(n)/2 - n*log(n) + n) print(r(10)) # 2.52759712035971761410261 print(r(1000)) # 2.506837169017401820105396 print(r(10**6)) # 2.506628483516698758337376 print(r(10**6)*exp(-1/(12*mpf(10)**6))) # 1er terme d'Euler–Maclaurin# 2.506628274631000502183361 print(sqrt(2*pi)) # 2.506628274631000502415765
le cas \(b=e\), \(x=1\) de la série de l'antilogarithme (L·28) : la base des logarithmes naturels, écrite comme une somme. C'est la plus rapide des formules pour \(e\) : 10 termes donnent 6 chiffres, 30 en donnent plus de 25. En alternant les signes, \(\sum(-1)^{n}/n!=1/e\) : c'est la probabilité qu'en distribuant des lettres au hasard, aucune n'arrive à son destinataire (les « dérangements », que le code compte pour 8 lettres).
# 10 puis 30 termes ; puis les dérangements de 8 lettres, comptés un par un from itertools import permutations from math import factorial as fact from mpmath import mp, mpf, factorial, e mp.dps = 30 # chiffres de garde pour les sommes s10 = sum(1/factorial(mpf(n)) for n in range(10)) s30 = sum(1/factorial(mpf(n)) for n in range(30)) mp.dps = 25 print(s10) # 2.718281525573192239858907 print(s30) # 2.718281828459045235360287 print(e) # 2.718281828459045235360287 d = sum(1 for p in permutations(range(8)) if all(p[i] != i for i in range(8))) print(mpf(d)/fact(8)) # 0.3678819444444444444444444 print(1/e) # 0.3678794411714423215955238
tirée de \(\ln(b-1)-\ln b=\ln\bigl(1-\tfrac1b\bigr)\) : en sommant pour \(b=2,\dots,n\), les logarithmes se simplifient en cascade et laissent \(\ln n\) ; en retranchant la série harmonique \(\tfrac12+\dots+\tfrac1n\), il reste \(1-\gamma\). Chaque terme est l'aire d'un petit croissant, entre la courbe \(1/x\) et le rectangle de hauteur \(1/b\) posé sur \([b-1,b]\). \(\gamma\approx0{,}5772\) est la constante d'Euler–Mascheroni, l'écart entre la série harmonique et le logarithme ; on ne sait même pas si elle est irrationnelle. Attention au cas \(b=1\) : \(\ln 0\) n'existe pas, d'où le départ en \(b=2\).
# mpmath, 25 chiffres ; puis la convergence lente de la définition de γ from mpmath import mp, nsum, log, inf, euler, harmonic mp.dps = 25 s = nsum(lambda b: log(b/(b-1)) - 1/b, [2, inf]) print(s) # 0.4227843350984671393934879 print(1 - euler) # 0.4227843350984671393934879 # la même constante, vue lentement : 1 + 1/2 + … + 1/n − ln n → γ print(harmonic(10**6) - log(10**6))# 0.5772161649014495272731789 print(euler) # 0.5772156649015328606065121
cinq constantes, \(e\), \(i\), \(\pi\), 1 et 0, dans une seule égalité. C'est le cas \(x=\pi\) de la formule d'Euler (L·33) : \(e^{i\pi}=\cos\pi+i\sin\pi=-1\). Et comme \(-1=i^{2}\), on a aussi \(\cos(n\pi)=(-1)^{n}=(i^{2})^{n}\) pour tout entier \(n\). On peut aussi la voir comme la série de l'antilogarithme (L·28) avec un exposant imaginaire : \(\sum_k (i\pi)^{k}/k!=-1\).
# e^{iπ} + 1, puis la série de l'exponentielle, puis cos(nπ) = (i²)ⁿ from mpmath import mp, exp, pi, j, cos, nsum, factorial, inf, nstr mp.dps = 25 print(nstr(abs(exp(j*pi) + 1), 3)) # zéro, au bruit d'arrondi près# 1.88e-26 print(nstr(nsum(lambda k: (j*pi)**k/factorial(k), [0, inf]), 20))# (-1.0 - 4.1642941483854355195e-54j) print(all(abs(cos(n*pi) - (j**2)**n) < 1e-20 for n in range(-10, 11)))# True
c'est la moitié du nombre d'or \(\varphi=(1+\sqrt5)/2\), qui gouverne le pentagone régulier. Les dix nombres \(e^{i\pi n/5}\), pour \(n=0\) à 9, sont les racines dixièmes de l'unité, les sommets d'un décagone (L·34) ; le premier vaut \(e^{i\pi/5}=\tfrac{1+\sqrt5}{4}+i\,\tfrac{\sqrt{10-2\sqrt5}}{4}\). Pour les angles \(\pi/n\), de telles formes en radicaux n'existent que pour certains \(n\) : 3, 4, 5, 6, 8, 10, 12, 15, 16, 17… (Gauss).
# mpmath, 25 chiffres ; puis e^{iπ/5} en radicaux from mpmath import mp, cos, sin, exp, pi, j, sqrt, nstr mp.dps = 25 print(cos(pi/5)) # 0.8090169943749474241022934 print((1+sqrt(5))/4) # 0.8090169943749474241022934 print(nstr(abs(exp(j*pi/5) - ((1+sqrt(5))/4 + j*sqrt(10-2*sqrt(5))/4)), 3))# 0.0
un imaginaire élevé à une puissance imaginaire donne un nombre réel : \(i^{i}=e^{i\ln i}=e^{i\cdot i\pi/2}=e^{-\pi/2}\approx0{,}2079\). C'est aussi \(1/\sqrt{e^{\pi}}\), l'inverse de la racine de la constante de Gelfond (IX·1), et donc un nombre transcendant (Gelfond, 1929) : les fractions ne peuvent que l'approcher, comme \(537842140/2587277449\), juste à 19 chiffres. Comme \(i=e^{i(\pi/2+2\pi n)}\) pour tout entier \(n\), \(i^{i}\) a en réalité une infinité de valeurs, toutes réelles : \(e^{-\pi/2-2\pi n}\). Au passage, \(\log_i x=\ln x/\ln i=2\ln x/(i\pi)\).
# la valeur principale, trois autres branches, deux approximations, puis log_i from mpmath import mp, mpf, j, pi, exp, sqrt, log, nstr mp.dps = 25 print(j**j) # (0.2078795763507619085469556 + 0.0j) print(exp(-pi/2)) # 0.2078795763507619085469556 print(1/sqrt(exp(pi))) # 0.2078795763507619085469556 print([nstr(exp(j*j*(pi/2 + 2*pi*n)), 6) for n in (-1, 0, 1)]) # d'autres valeurs, toutes réelles# ['(111.318 + 0.0j)', '(0.20788 + 0.0j)', '(0.000388203 + 0.0j)'] print(nstr(mpf(537842140)/2587277449 - exp(-pi/2), 3))# 1.87e-20 print(nstr(29/(42 + sqrt(9507)) - exp(-pi/2), 3))# -3.02e-9 print(nstr(abs(log(2)/log(j) - 2*log(2)/(j*pi)), 3))# 0.0
car \(\tfrac{1-i}{1+i}=-i\) et \(\ln(-i)=-i\pi/2\). Même idée, plus directe : \(\pi=-i\ln(-1)\), puisque \(e^{i\pi}=-1\) (XII·1). Formule attribuée à Fagnano, au XVIIIᵉ siècle. Le logarithme complexe a une infinité de valeurs, écartées de \(2\pi i\) ; on prend ici la valeur principale, celle dont la partie imaginaire est entre \(-\pi\) et \(\pi\).
# valeur principale du logarithme complexe, 25 chiffres from mpmath import mp, j, log, pi mp.dps = 25 print((1-j)/(1+j)) # (0.0 - 1.0j) print(2*j*log((1-j)/(1+j))) # (3.141592653589793238462643 + 0.0j) print(-j*log(-1)) # (3.141592653589793238462643 + 0.0j) print(pi) # 3.141592653589793238462643
de \(x^{2}=x+1\), on tire \(x=1+1/x\), et en remplaçant \(x\) par lui-même à l'infini, cette fraction continue faite uniquement de 1. Ses réduites sont les rapports de nombres de Fibonacci : 1, 2, 3/2, 5/3, 8/5, 13/8… Ce sont les plus lentes qui soient, car des 1 partout, c'est le pire cas : \(\varphi\) est le nombre le plus difficile à approcher par des fractions. Avec \(x^{2}+x-1=0\), on obtient de même \(1/\varphi=[0;1,1,1,\dots]\) ; mais l'autre racine, \(-\varphi\), échappe à l'itération \(x\leftarrow1/(1+x)\), qui la repousse et file vers \(1/\varphi\). Quant à \(-1/\varphi\), c'est le « conjugué » de \(\varphi\) au sens de Galois (L·49).
# itération, réduites de Fibonacci, puis l'autre racine qui s'échappe from fractions import Fraction from mpmath import mp, mpf, sqrt, nstr mp.dps = 25 x = mpf(1) for _ in range(60): x = 1 + 1/x # φ = 1 + 1/φ, itéré print(x) # 1.618033988749894848204587 print((1 + sqrt(5))/2) # 1.618033988749894848204587 r = [Fraction(1)] for _ in range(7): r.append(1 + 1/r[-1]) print([str(f) for f in r]) # rapports de Fibonacci# ['1', '2', '3/2', '5/3', '8/5', '13/8', '21/13', '34/21'] y = mpf('-1.618') # on part tout près de −φ… for _ in range(60): y = 1/(1 + y) print(nstr(y, 10)) # … et l'itération file vers 1/φ# 0.6180339887
la fraction continue de \(e\) n'est pas périodique (\(e\) n'est pas racine d'une équation du second degré, L·48), mais elle suit un motif parfait par triplets \(1,2k,1\), découvert par Euler (1737). Celle de \(\pi\), à l'inverse, semble n'avoir aucun ordre : \([3;7,15,1,292,1,1,1,2,\dots]\). Couper juste avant le grand 292 donne \(355/113\), juste à 7 chiffres ; c'est la fraction de Zu Chongzhi (Vᵉ siècle), le même que pour le volume de la boule (L·10).
# les termes de e, le motif qui la reconstruit, puis ceux de π et 355/113 from mpmath import mp, mpf, e, pi, floor, nstr mp.dps = 60 # chiffres de garde pour extraire les termes def fc(x, n): # les n premiers termes de la fraction continue t = [] for _ in range(n): a = int(floor(x)); t.append(a); x = 1/(x - a) return t def valeur(t): # [a0; a1, …, ak] évaluée depuis la fin v = mpf(t[-1]) for a in reversed(t[:-1]): v = a + 1/v return v print(fc(e, 16)) # [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1] motif = [2] + [a for k in range(1, 20) for a in (1, 2*k, 1)] print(nstr(valeur(motif) - e, 3)) # le motif redonne e# 4.76e-59 print(fc(pi, 12)) # [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1] print(nstr(mpf(355)/113 - pi, 3))# 2.67e-7
on tire des nombres au hasard entre 0 et 1 et on les additionne jusqu'à dépasser 1 : il en faut en moyenne exactement \(e\). La raison : avoir besoin de plus de \(n\) tirages, c'est que \(U_1+\dots+U_n\leq1\), et la probabilité en vaut \(1/n!\), le volume d'un simplexe (le coin du cube de dimension \(n\)). La moyenne est la somme de ces probabilités, c'est-à-dire la série de \(e\) (XI·1). À ne pas confondre avec \(1+\tfrac12+\dots+\tfrac18=2{,}71786\ldots\), qui ne tombe près de \(e\) que par coïncidence (\(H_8\approx\ln8+\gamma+\tfrac1{16}\)).
# un million de tirages ; P(N > n) contre 1/n! ; la somme exacte ; puis la coïncidence H₈ import random from fractions import Fraction from mpmath import mp, mpf, factorial, nsum, inf, e random.seed(1) def tirages(): total, n = 0.0, 0 while total <= 1: total += random.random(); n += 1 return n T = [tirages() for _ in range(10**6)] print(sum(T) / len(T)) # simulation : un million d'essais# 2.718362 print([round(sum(t > n for t in T) / len(T), 4) for n in range(1, 5)])# [1.0, 0.5004, 0.1664, 0.0417] print([round(1/float(factorial(n)), 4) for n in range(1, 5)]) # P(N > n) = 1/n!# [1.0, 0.5, 0.1667, 0.0417] mp.dps = 25 print(nsum(lambda n: 1/factorial(n), [0, inf]))# 2.718281828459045235360287 print(e) # 2.718281828459045235360287 H8 = sum(Fraction(1, k) for k in range(1, 9)) print(H8, float(e - H8)) # la coïncidence# 761/280 0.0004246856019023782
pour \(a>0,\ a\neq1\) et \(M,N>0\). Sans cette dernière condition, c'est faux : \(M=N=-1\) donne \(\log_a 1=0\), alors que \(\log_a(-1)\) n'existe pas dans \(\mathbb R\).
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal from mpmath import mp, mpf, log, nstr import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): a, M, N = base(), pos(), pos() err = max(err, abs(log(M*N, a) - (log(M, a) + log(N, a)))) print(nstr(err, 3)) # 6.62e-24
pour \(a>0,\ a\neq1\) et \(M,N>0\).
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal from mpmath import mp, mpf, log, nstr import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): a, M, N = base(), pos(), pos() err = max(err, abs(log(M/N, a) - (log(M, a) - log(N, a)))) print(nstr(err, 3)) # 3.31e-24
pour \(a>0,\ a\neq1\), \(M>0\) et \(q\) réel quelconque.
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal from mpmath import mp, mpf, log, nstr import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): a, M = base(), pos() q = mpf(random.uniform(-10, 10)) err = max(err, abs(log(M**q, a) - q*log(M, a))) print(nstr(err, 3)) # 2.65e-23
pour \(a,b>0\), \(a,b\neq1\) et \(N>0\). C'est elle qui permet de tout calculer avec \(\ln\) seul.
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal from mpmath import mp, mpf, log, nstr import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): a, b, N = base(), base(), pos() err = max(err, abs(log(N, a) - log(N, b)/log(a, b))) print(nstr(err, 3)) # 6.62e-24
pour \(N>0\). C'est une notation : sans base écrite, \(\log\) désigne la base 10 (sur la calculatrice, pas toujours ailleurs).
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal from mpmath import mp, mpf, log, nstr, log10 import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): N = pos() err = max(err, abs(log(N, 10) - log10(N))) print(nstr(err, 3)) # 0
pour \(N>0\). Encore une notation, pour la base \(e\).
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal from mpmath import mp, mpf, log, nstr, e, ln import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): N = pos() err = max(err, abs(log(N, e) - ln(N))) print(nstr(err, 3)) # 1.03e-25
pour \(a>0,\ a\neq1\) et \(N\) réel quelconque. Avec \(a=e\) : \(\ln e^{N}=N\).
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal from mpmath import mp, mpf, log, nstr import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): a = base() N = mpf(random.uniform(-20, 20)) err = max(err, abs(log(a**N, a) - N)) print(nstr(err, 3)) # 4.14e-25
pour \(a>0,\ a\neq1\) et \(N>0\) (ici \(N\) doit être positif, contrairement à L·7). Avec \(a=e\) : \(e^{\ln N}=N\).
# 1000 tirages au hasard dans le domaine ; on affiche l'écart maximal (relatif) from mpmath import mp, mpf, log, nstr import random mp.dps = 25 random.seed(1) pos = lambda: mpf(random.uniform(0.01, 100)) base = lambda: random.choice([pos(), mpf(random.uniform(0.01, 0.99))]) err = 0 for _ in range(1000): a, N = base(), pos() err = max(err, abs(a**log(N, a) - N)/N) print(nstr(err, 3)) # 6.04e-26
pour tout entier \(n\geq1\). Plus généralement, \(\varphi(n^{k})=n^{k-1}\varphi(n)\) : élever au carré ne crée aucun nouveau facteur premier.
# vérifiée exhaustivement pour n < 10 000 (la preuve tient en deux lignes) def phi(n): r, m, p = n, n, 2 while p*p <= m: if m % p == 0: while m % p == 0: m //= p r -= r // p p += 1 return r - r//m if m > 1 else r print(all(phi(n*n) == n*phi(n) for n in range(1, 10**4)))# True
pour tout entier \(n\geq0\) et \(r>0\). La fonction \(\Gamma\) prolonge la factorielle, avec \(\Gamma(k+1)=k!\) et \(\Gamma(\tfrac12)=\sqrt\pi\) (voir IX·2) : c'est elle qui traite d'un coup les dimensions paires et impaires.
| n | volume de la boule | aire de la sphère |
|---|---|---|
| 1 | \(2r\) | \(2\) |
| 2 | \(\pi r^{2}\) | \(2\pi r\) |
| 3 | \(\tfrac43\pi r^{3}\) | \(4\pi r^{2}\) |
| 4 | \(\tfrac12\pi^{2} r^{4}\) | \(2\pi^{2} r^{3}\) |
| 5 | \(\tfrac{8}{15}\pi^{2} r^{5}\) | \(\tfrac83\pi^{2} r^{4}\) |
Une surprise : pour \(r=1\), le volume augmente jusqu'à la dimension 5, puis décroît vers 0. En grande dimension, presque tout le volume se loge contre le bord : une coquille d'épaisseur 1 % en contient \(1-0{,}99^{n}\), soit 5 % en dimension 5 et 63 % en dimension 100. Une orange qui ne serait presque plus que de la peau.
# vérification indépendante par découpage en tranches, pour n = 1 à 10 from mpmath import mp, mpf, pi, gamma, quad, sqrt, nstr mp.dps = 25 V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n # une boule de dimension n = un empilement de tranches, # chacune étant une boule de dimension n − 1 err = max(abs(quad(lambda x: V(n-1, sqrt(1-x**2)), [-1, 1]) - V(n)) for n in range(1, 11)) print(nstr(err, 3)) # 1.03e-25 nmax = max(range(30), key=V) print(nmax, V(nmax)) # 5 5.263789013914324596711729
pour \(n\geq1\) et \(r>0\). La peau, c'est ce qu'on ajoute quand \(r\) grandit d'un cran. Ici \(S_n\) borde la boule de \(\mathbb R^{n}\) ; les mathématiciens la notent souvent \(\mathbb S^{n-1}\), d'après sa propre dimension. L'aire de la sphère unité culmine en dimension 7, à environ 33,07.
# dérivée numérique de V comparée à S, 200 tirages (écart relatif) from mpmath import mp, mpf, pi, gamma, diff, nstr import random mp.dps = 25 random.seed(1) V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n S = lambda n, r=1: 2*pi**(mpf(n)/2) / gamma(mpf(n)/2) * r**(n-1) err = 0 for _ in range(200): n, r = random.randint(1, 12), mpf(random.uniform(0.1, 5)) err = max(err, abs(diff(lambda t: V(n, t), r) - S(n, r)) / S(n, r)) print(nstr(err, 3)) # 7.5e-26 nmax = max(range(1, 30), key=S) print(nmax, S(nmax)) # 7 33.07336179231980818717474
pour \(n\geq2\), avec \(V_0=1\) et \(V_1=2\). De quoi tout retrouver sans \(\Gamma\), et voir pourquoi le volume finit par décroître : dès que \(n>2\pi\), le facteur \(2\pi/n\) passe sous 1.
# vérifiée pour n = 2 à 59 from mpmath import mp, mpf, pi, gamma, nstr mp.dps = 25 V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n err = max(abs(V(n) - 2*pi/n * V(n-2)) for n in range(2, 60)) print(nstr(err, 3)) # 5.17e-26
pour un cube de côté \(a\) en dimension \(n\geq1\), de volume \(a^{n}\) : ses \(2n\) faces sont des cubes de dimension \(n-1\). Comme pour la boule (L·11), le bord est la dérivée du volume, à condition de dériver par rapport au rayon de la boule inscrite, \(r=a/2\), car \(V=(2r)^{n}\) (voir L·15). La grande diagonale mesure \(a\sqrt n\) et grandit sans fin.
| n | volume | bord | sommets | arêtes |
|---|---|---|---|---|
| 1 | \(a\) | \(2\) | 2 | 1 |
| 2 | \(a^{2}\) | \(4a\) | 4 | 4 |
| 3 | \(a^{3}\) | \(6a^{2}\) | 8 | 12 |
| 4 | \(a^{4}\) | \(8a^{3}\) | 16 | 32 |
| 5 | \(a^{5}\) | \(10a^{4}\) | 32 | 80 |
# dérivée numérique de V par rapport à r = a/2, 200 tirages (écart relatif) from mpmath import mp, mpf, diff, nstr import random mp.dps = 25 random.seed(1) V = lambda n, r: (2*r)**n # r = a/2 : rayon de la boule inscrite S = lambda n, a: 2*n*a**(n-1) err = 0 for _ in range(200): n, a = random.randint(1, 12), mpf(random.uniform(0.1, 5)) err = max(err, abs(diff(lambda r: V(n, r), a/2) - S(n, a)) / S(n, a)) print(nstr(err, 3)) # 2.53e-26
nombre de faces de dimension \(k\) du cube de dimension \(n\) : une face fixe \(n-k\) coordonnées, chacune à 0 ou à 1, et laisse les \(k\) autres libres. Pour le tesseract : 16 sommets, 32 arêtes, 24 carrés, 8 cubes. En comptant tout, cube compris, \(\sum_k f_k(n)=3^{n}\) : chaque coordonnée a trois états.
# énumération exhaustive de toutes les faces, n = 1 à 8 (la dernière ligne montre n = 8) from itertools import product from math import comb ok = True for n in range(1, 9): f = [0]*(n+1) for face in product('01*', repeat=n): # 0, 1 : fixée ; * : libre f[face.count('*')] += 1 ok &= all(f[k] == 2**(n-k)*comb(n, k) for k in range(n+1)) ok &= sum(f) == 3**n print(ok) # True print(f) # [256, 1024, 1792, 1792, 1120, 448, 112, 16, 1]
pour tout solide de \(\mathbb R^{n}\) dont chaque face touche une même sphère de rayon \(r\), et pour la boule elle-même : on le découpe en pyramides de sommet le centre et de hauteur \(r\), chacune de volume \(r\times\text{base}/n\). Comme \(V\) est proportionnel à \(r^{n}\), cela revient à \(S=dV/dr\), ce qui explique à la fois L·11 et L·13.
# trois familles de solides (boule, cube, hyperoctaèdre), n = 1 à 12 from mpmath import mp, mpf, pi, gamma, sqrt, factorial, nstr mp.dps = 25 err = 0 for n in range(1, 13): # boule de rayon 1 V, S, r = pi**(mpf(n)/2)/gamma(mpf(n)/2+1), 2*pi**(mpf(n)/2)/gamma(mpf(n)/2), 1 err = max(err, abs(V - r*S/n)/V) # cube de côté 2, boule inscrite de rayon 1 V, S, r = mpf(2)**n, 2*n*mpf(2)**(n-1), 1 err = max(err, abs(V - r*S/n)/V) # hyperoctaèdre |x1|+…+|xn| ≤ 1 : 2ⁿ facettes, rayon inscrit 1/√n V, S, r = mpf(2)**n/factorial(n), mpf(2)**n*sqrt(n)/factorial(n-1), 1/sqrt(n) err = max(err, abs(V - r*S/n)/V) print(nstr(err, 3)) # 2.47e-26
dans le cube \([-2,2]^{n}\), on place \(2^{n}\) boules de rayon 1 centrées en \((\pm1,\dots,\pm1)\), et \(\rho_n\) est le rayon de la plus grande boule centrale qui tient entre elles. Dimension 2 : \(\rho\approx0{,}41\). Dimension 4 : \(\rho=1\), aussi grosse que les autres. Dimension 9 : \(\rho=2\), elle touche les faces du cube. Dès la dimension 10, elle en dépasse, alors que toutes les boules des coins restent dedans. En grande dimension, le cube ressemble moins à une boîte qu'à un oursin.
# rayon de la boule centrale ; au-delà de 2, elle sort du cube from math import sqrt # centre → centre d'une boule de coin : √n ; on retire son rayon 1 rho = {n: sqrt(n) - 1 for n in range(1, 21)} print(rho[4]) # 1.0 print(rho[9]) # 2.0 print(min(n for n in rho if rho[n] > 2))# 10
pour tout entier \(n\geq0\) ; \(n=0\) redonne IX·2. C'est par là que \(\sqrt\pi\) entre dans le volume des boules de dimension impaire, et seulement dans celles-là (L·10).
# n = 0 à 39, écart relatif from mpmath import mp, mpf, gamma, factorial, sqrt, pi, nstr mp.dps = 25 err = max(abs(gamma(n + mpf(1)/2) - factorial(2*n)/(4**n*factorial(n))*sqrt(pi)) / gamma(n + mpf(1)/2) for n in range(0, 40)) print(nstr(err, 3)) # 2.28e-26
pour tout \(s\) hors des pôles \(0,-\tfrac12,-1,\dots\), complexe compris. C'est la loi des demi-pas : elle soude les valeurs de \(\Gamma\) décalées de \(\tfrac12\), donc les dimensions paires et impaires, et c'est d'elle que vient le \(-\tfrac12\) du pic des boules (A121546).
# 200 points complexes tirés au hasard, écart relatif from mpmath import mp, mpf, mpc, gamma, sqrt, pi, nstr import random mp.dps = 25 random.seed(1) err = 0 for _ in range(200): s = mpc(random.uniform(0.05, 8), random.uniform(-5, 5)) l = gamma(s)*gamma(s + mpf(1)/2) err = max(err, abs(l - 2**(1-2*s)*sqrt(pi)*gamma(2*s))/abs(l)) print(nstr(err, 3)) # 4.55e-26
pour tout \(s\) non entier, complexe compris. En \(s=\tfrac12\) : \(\Gamma(\tfrac12)^{2}=\pi\), c'est-à-dire IX·2.
# 200 points complexes tirés au hasard, écart relatif from mpmath import mp, mpc, gamma, sin, pi, nstr import random mp.dps = 25 random.seed(1) err = 0 for _ in range(200): s = mpc(random.uniform(-6, 6), random.uniform(-3, 3)) l = gamma(s)*gamma(1-s) err = max(err, abs(l - pi/sin(pi*s))/abs(l)) print(nstr(err, 3)) # 3.16e-25
où \(S(s)=\dfrac{2\pi^{s/2}}{\Gamma(s/2)}\) est l'aire de la sphère unité qui borde la boule de dimension \(s\) (L·11), une dimension désormais complexe. C'est l'équation de Riemann (1859), \(\xi(s)=\xi(1-s)\) avec \(\xi(s)=\tfrac12 s(s-1)\,\pi^{-s/2}\Gamma(s/2)\,\zeta(s)\), réécrite : ζ divisée par l'aire de la sphère est symétrique autour de \(\tfrac12\). Exemple : la sphère de dimension −1 a une « aire » de \(-1/\pi\), et Bâle, \(\zeta(2)=\pi^{2}/6\), donne aussitôt \(\zeta(-1)=-\tfrac1{12}\).
B. Riemann, « Ueber die Anzahl der Primzahlen unter einer gegebenen Grösse », Monatsberichte der Berliner Akademie, 1859 : l'équation fonctionnelle, sous la forme en \(\xi\).
J. Tate, Fourier analysis in number fields and Hecke's zeta-functions, thèse, Princeton, 1950 ; publiée dans J. W. S. Cassels et A. Fröhlich (éd.), Algebraic Number Theory, Academic Press, 1967. Le facteur \(\pi^{-s/2}\Gamma(s/2)\) y devient le facteur local de ζ à la place réelle : l'intégrale de la gaussienne \(e^{-\pi x^{2}}\) contre \(|x|^{s}\), la même gaussienne qui donne l'aire des sphères.
A. Karlsson et M. Pallich, « Volumes of spheres and special values of zeta functions of ℤ and ℤ/nℤ », 2022, arXiv:2209.03590. Lecture voisine mais différente : l'aire de la sphère unité de chaque dimension y est un produit de valeurs spéciales d'une autre fonction zêta, celle de ℤ, qui a elle aussi une symétrie \(s\leftrightarrow1-s\).
À notre connaissance, la réécriture \(\zeta(s)/S(s)=\zeta(1-s)/S(1-s)\) et le calcul de \(\zeta(-1)\) par la « sphère de dimension −1 » ne figurent tels quels dans aucune de ces sources ; ils en découlent en une ligne. À ne pas confondre avec la sphère de Riemann, qui désigne tout autre chose : le plan complexe complété par un point à l'infini.
# 200 points complexes tirés au hasard, puis ζ(−1) retrouvé depuis Bâle from mpmath import mp, mpc, zeta, gamma, pi, nstr import random mp.dps = 25 random.seed(1) S = lambda s: 2*pi**(s/2)/gamma(s/2) # aire de la sphère, L·11 err = 0 for _ in range(200): s = mpc(random.uniform(-8, 9), random.uniform(-20, 20)) l = zeta(s)/S(s) err = max(err, abs(l - zeta(1-s)/S(1-s))/abs(l)) print(nstr(err, 3)) # 1.37e-25 print(S(-1)) # -0.3183098861837906715377675 print(zeta(2)/S(2) * S(-1)) # -0.08333333333333333333333333
densité maximale d'un empilement de boules égales en dimension 8, atteinte par le réseau \(E_8\) : les points de \(\mathbb R^{8}\) à coordonnées toutes entières ou toutes demi-entières, de somme paire. Ses 240 vecteurs les plus courts sont les « racines » de \(E_8\) ; chaque boule en touche donc 240 autres, et ces 240 points sont les sommets d'un solide de dimension 8, le polytope de Gosset \(4_{21}\). Ses symétries forment un groupe fini de 696 729 600 éléments ; le groupe de Lie \(E_8\), lui, est de dimension 248. La densité vaut \(V_8(1/\sqrt2)\), un simple volume de L·10. Que personne ne puisse faire mieux, c'est le théorème de Viazovska : Idem le cite, le code en vérifie les ingrédients.
M. Viazovska, « The sphere packing problem in dimension 8 », Annals of Mathematics 185 (2017), arXiv:1603.04246. Médaille Fields 2022.
# construction du réseau, comptage des voisines, covolume et densité from itertools import product from mpmath import mp, mpf, pi, gamma, sqrt, matrix, det mp.dps = 25 V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n # L·10 # coordonnées doublées w = 2v : toutes paires ou toutes impaires, Σw ≡ 0 mod 4 in_e8 = lambda w: (all(x % 2 == 0 for x in w) or all(x % 2 for x in w)) and sum(w) % 4 == 0 racines = [w for w in product(range(-2, 3), repeat=8) if sum(x*x for x in w) == 8 and in_e8(w)] # |v|² = 2 print(len(racines)) # 240 B = matrix([[2,0,0,0,0,0,0,0], [-1,1,0,0,0,0,0,0], [0,-1,1,0,0,0,0,0], [0,0,-1,1,0,0,0,0], [0,0,0,-1,1,0,0,0], [0,0,0,0,-1,1,0,0], [0,0,0,0,0,-1,1,0], [mpf(1)/2]*8]) # une base de E8 print(abs(det(B))) # 1.0 print(V(8, 1/sqrt(2)) / abs(det(B)))# 0.2536695079010480136365634 print(pi**4/384) # 0.2536695079010480136365634
pour tout \(n\geq1\), où \(\sigma_3(n)\) est la somme des cubes des diviseurs de \(n\) : 240 points à distance \(\sqrt2\), puis 2 160, 6 720, 17 520… La raison est profonde : la série qui compte ces points (la série thêta de \(E_8\)) est une forme modulaire, la série d'Eisenstein \(E_4\). Ce sont ces mêmes formes modulaires qui portent la preuve de Viazovska (L·21).
J.-P. Serre, Cours d'arithmétique, PUF, 1970, chap. VII : séries thêta et formes modulaires, dont celle de \(E_8\).
# dénombrement exact des points du réseau, couche par couche, n = 1 à 5 from collections import Counter def compte(valeurs, qmax): # w = 2v ; état (Σw², Σw mod 4) etats = Counter({(0, 0): 1}) for _ in range(8): nouv = Counter() for (q, s), c in etats.items(): for w in valeurs: if q + w*w <= qmax: nouv[(q + w*w, (s + w) % 4)] += c etats = nouv return etats N = 5 entiers = compte(range(-8, 9, 2), 8*N) # v entier demi = compte(range(-9, 10, 2), 8*N) # v demi-entier sigma3 = lambda n: sum(d**3 for d in range(1, n+1) if n % d == 0) print([entiers[(8*n, 0)] + demi[(8*n, 0)] for n in range(1, N+1)])# [240, 2160, 6720, 17520, 30240] print([240*sigma3(n) for n in range(1, N+1)])# [240, 2160, 6720, 17520, 30240]
densité maximale en dimension 24, atteinte par le réseau de Leech : covolume 1, aucun vecteur de longueur \(\sqrt2\), et 196 560 vecteurs de longueur 2, donc 196 560 voisines pour chaque boule de rayon 1. La densité est alors exactement \(V_{24}(1)=\pi^{12}/12!\). Le code tire ces nombres de la série thêta de Leech, \(\theta=E_{12}-\tfrac{65520}{691}\Delta\), identité de formes modulaires que l'on cite ; l'optimalité est un théorème de 2016, cité lui aussi.
H. Cohn, A. Kumar, S. D. Miller, D. Radchenko, M. Viazovska, « The sphere packing problem in dimension 24 », Annals of Mathematics 185 (2017), arXiv:1603.06518.
# nombre de vecteurs de norme² 2, 4, 6 via la série thêta, puis densité from fractions import Fraction from mpmath import mp, mpf, pi, gamma, factorial mp.dps = 25 V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n # L·10 N = 3 P = [1] + [0]*N # ∏ (1 − qⁿ)²⁴, jusqu'à q^N for n in range(1, N+1): for _ in range(24): for k in range(N, n-1, -1): P[k] -= P[k-n] tau = [0] + P[:N] # Δ = q·∏ : τ(1), τ(2), τ(3) sigma11 = lambda n: sum(d**11 for d in range(1, n+1) if n % d == 0) c = Fraction(65520, 691) print([int(c*(sigma11(n) - tau[n])) for n in range(1, N+1)])# [0, 196560, 16773120] print(V(24)) # 0.001929574309403923047903346 print(pi**12/factorial(12)) # 0.001929574309403923047903346
avec \(T(0,0)=1\). Le même moule, un seul poids qui change : où placer le \(n\)-ième objet ? Soit il ouvre un groupe à lui seul (premier terme), soit il rejoint l'existant, et le poids compte de combien de façons : \(k\) groupes où entrer, \(n-1\) places dans les cycles, \(n+k-1\) places dans les files. Le code ne fait pas confiance à la récurrence : il énumère réellement les objets jusqu'à \(n=6\).
| triangle | poids \(w(n,k)\) | compte | somme des lignes |
|---|---|---|---|
| Pascal | \(1\) | sous-ensembles de taille \(k\) | \(2^{n}\) |
| Stirling 2ᵉ | \(k\) | répartitions en \(k\) groupes | Bell |
| Stirling 1ʳᵉ | \(n-1\) | permutations à \(k\) cycles | \(n!\) |
| Lah | \(n+k-1\) | répartitions en \(k\) files | 1, 3, 13, 73, 501… |
# énumération réelle des sous-ensembles, partitions, permutations et files, n = 0 à 6 def partitions(elems): # toutes les partitions d'un ensemble if not elems: yield [] return x, reste = elems[0], elems[1:] for p in partitions(reste): yield [[x]] + p # x seul dans un bloc for i in range(len(p)): # ou x rejoint un bloc yield p[:i] + [[x] + p[i]] + p[i+1:] from itertools import combinations, permutations, product from math import factorial def triangle(w, N): T = [[1]] for n in range(1, N+1): prev = T[-1] + [0] T.append([(prev[k-1] if k else 0) + w(n, k)*prev[k] for k in range(n+1)]) return T def cycles(p): # nombre de cycles d'une permutation vu, c = set(), 0 for i in range(len(p)): if i not in vu: c += 1 while i not in vu: vu.add(i); i = p[i] return c def compte(n, famille): t = [0]*(n+1) if famille == 'pascal': for k in range(n+1): t[k] = sum(1 for _ in combinations(range(n), k)) if famille == 'stirling2': for P in partitions(list(range(n))): t[len(P)] += 1 if famille == 'stirling1': for p in permutations(range(n)): t[cycles(p)] += 1 if famille == 'lah': # chaque groupe mis en file for P in partitions(list(range(n))): t[len(P)] += sum(1 for _ in product(*[permutations(b) for b in P])) return t poids = {'pascal': lambda n, k: 1, 'stirling2': lambda n, k: k, 'stirling1': lambda n, k: n-1, 'lah': lambda n, k: n+k-1} print(all(compte(n, f) == triangle(w, 6)[n] for f, w in poids.items() for n in range(7)))# True print({f: sum(triangle(w, 6)[6]) for f, w in poids.items()})# {'pascal': 64, 'stirling2': 203, 'stirling1': 720, 'lah': 4051}
pour tout entier \(k\geq1\). Les nombres de Bernoulli \(B_{2k}\) (1/6, −1/30, 1/42…) sont les coefficients des sommes de puissances \(1^{p}+2^{p}+\dots+n^{p}\) (Faulhaber) : c'est le fil « sommes de puissances → π ». \(k=1\) donne Bâle (II·1), \(k=2\) donne \(\pi^{4}/90\) (II·2). Côté négatif, \(\zeta(1-2k)=-B_{2k}/2k\) : c'est le \(-1/12\) de L·20. Pour \(\zeta(3),\zeta(5)\dots\), aucune formule de ce genre n'est connue.
# k = 1 à 15, écart relatif ; puis B₂, B₄, B₆ en fractions from mpmath import mp, zeta, bernoulli, bernfrac, pi, factorial, nstr mp.dps = 25 err = max(abs(zeta(2*k) - (-1)**(k+1)*bernoulli(2*k)*(2*pi)**(2*k)/(2*factorial(2*k))) / zeta(2*k) for k in range(1, 16)) print(nstr(err, 3)) # 1.81e-25 print([bernfrac(2*k) for k in (1, 2, 3)])# [(1, 6), (-1, 30), (1, 42)]
pour tout entier \(k\geq0\), où \(\beta\) est la fonction de Dirichlet de III·1. Les \(|E_{2k}|\) (1, 1, 5, 61, 1385…) comptent les permutations en zigzag, qui montent et descendent en alternance. \(k=0\) donne Leibniz (I·2), \(k=1\) donne \(\pi^{3}/32\) (III·2). C'est le miroir exact de L·25 : pour β, ce sont les valeurs paires qui résistent, à commencer par \(\beta(2)=G\), la constante de Catalan.
# k = 0 à 9, écart relatif ; puis les premiers nombres d'Euler from mpmath import mp, mpf, nsum, inf, eulernum, pi, factorial, nstr mp.dps = 25 beta = lambda s: nsum(lambda j: mpf(-1)**j/(2*j+1)**s, [0, inf]) err = max(abs(beta(2*k+1) - (-1)**k*eulernum(2*k)*pi**(2*k+1)/(4**(k+1)*factorial(2*k))) / beta(2*k+1) for k in range(0, 10)) print(nstr(err, 3)) # 1.29e-25 print([int(eulernum(2*k)) for k in range(6)])# [1, -1, 5, -61, 1385, -50521]
pour \(n\geq1\), où \(W_n=\int_0^{\pi/2}\sin^{n}x\,dx\). En pair et impair : \(W_{2m}=\tfrac{\pi}{2}\binom{2m}{m}4^{-m}\) contient π, \(W_{2m+1}=4^{m}\big/\big((2m+1)\binom{2m}{m}\big)\) n'en contient pas. C'est le fil caché de plusieurs fiches : couper une boule en tranches (Zu, L·10) donne \(V_n=2W_nV_{n-1}\), et deux crans d'un coup donnent la récurrence \(V_n=\tfrac{2\pi}{n}V_{n-2}\) de L·12 ; le rapport \(W_{2m}/W_{2m+1}\to1\) donne le produit de Wallis (IV·1) ; les Steinmetz de dimension \(n\) en sont des multiples, d'où « rationnel en dimension impaire, π en dimension paire » ; et le binomial central de X·3 y apparaît.
# n = 1 à 20 : le produit, le lien avec les boules, puis les deux formes closes from mpmath import mp, mpf, quad, sin, pi, gamma, binomial, nstr mp.dps = 25 W = lambda n: quad(lambda x: sin(x)**n, [0, pi/2]) V = lambda n: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) # L·10 print(nstr(max(abs(W(n)*W(n-1) - pi/(2*n)) for n in range(1, 21)), 3))# 6.46e-27 print(nstr(max(abs(V(n) - 2*W(n)*V(n-1)) for n in range(1, 21)), 3))# 1.03e-25 print(nstr(max(abs(W(2*m) - pi/2*binomial(2*m, m)/4**m) for m in range(10)), 3))# 6.46e-27 print(nstr(max(abs(W(2*m+1) - 4**m/((2*m+1)*binomial(2*m, m))) for m in range(10)), 3))# 1.29e-26
pour \(b>0\) et \(x\) réel : c'est le développement de Maclaurin de l'exponentielle, puisque \(b^{x}=e^{x\ln b}\). Les anciens manuels l'appellent « formule de Mac-Laurin » et l'écrivent avec « log b » : ce log doit être le logarithme naturel. Avec le log décimal (L·5), \(b=10\), \(x=1\) donnerait \(e\) au lieu de 10. Les factorielles écrasent vite les termes, d'autant plus que \(x\ln b\) est petit : d'où l'intérêt de se ramener d'abord à un petit exposant. Avec \(b=e\), \(x=1\), on obtient XI·1.
# 200 tirages au hasard ; nombre de termes selon la taille de x·ln b ; puis le piège du log décimal from mpmath import mp, mpf, ln, log10, factorial, nstr import random mp.dps = 35 # chiffres de garde : si x·ln b < 0, les termes alternent et se compensent random.seed(1) def antilog(b, x, N=200): t = s = mpf(1) for k in range(1, N): t *= x*ln(b)/k s += t return s err = 0 for _ in range(200): b, x = mpf(random.uniform(0.1, 10)), mpf(random.uniform(-5, 5)) err = max(err, abs(antilog(b, x) - b**x) / b**x) print(nstr(err, 3)) # 3.36e-28 def termes(u): # termes nécessaires pour 25 chiffres t, k = mpf(1), 0 while abs(t) > mpf(10)**-25: k += 1 t *= u/k return k print(termes(mpf('0.1')), termes(mpf(10)))# 15 64 print(nstr(sum(log10(10)**k/factorial(k) for k in range(60)), 12))# 2.71828182846
pour tout entier \(n\geq1\). Une preuve en une image : on écrit les \(n\) fractions \(\tfrac1n,\tfrac2n,\dots,\tfrac nn\) et on les simplifie ; chacune prend un dénominateur \(d\) qui divise \(n\), et il y en a exactement \(\varphi(d)\) pour chaque \(d\). Pour \(n=6\) : 1 + 1 + 2 + 2 = 6. Avec L·30, elle donne \(\sum_{d\mid n}\Lambda(d)=\ln\bigl(\sum_{d\mid n}\varphi(d)\bigr)\), les deux côtés valant \(\ln n\).
# vérifiée exactement pour n = 1 à 3000 ; puis le détail pour n = 6 def facteurs(n): # {premier: exposant} f, p = {}, 2 while p*p <= n: while n % p == 0: f[p] = f.get(p, 0) + 1 n //= p p += 1 if n > 1: f[n] = f.get(n, 0) + 1 return f diviseurs = lambda n: [d for d in range(1, n+1) if n % d == 0] def phi(n): r = n for p in facteurs(n): r -= r // p return r print(all(sum(phi(d) for d in diviseurs(n)) == n for n in range(1, 3001)))# True print([phi(d) for d in diviseurs(6)])# [1, 1, 2, 2]
pour tout entier \(n\geq1\), où \(\Lambda(d)=\ln p\) si \(d\) est une puissance d'un nombre premier \(p\), et 0 sinon. Chaque premier « pèse » son logarithme : si \(p^{a}\) divise exactement \(n\), les diviseurs \(p,p^{2},\dots,p^{a}\) apportent \(a\ln p\), et la somme sur tous les premiers redonne \(\ln n\). Combinée à L·29, elle s'écrit \(\sum_{d\mid n}\Lambda(d)=\ln\bigl(\sum_{d\mid n}\varphi(d)\bigr)\). \(\Lambda\) est aussi la fonction dont \(-\zeta'/\zeta\) fait la somme : \(\sum\Lambda(n)/n^{s}=-\zeta'(s)/\zeta(s)\).
# n = 2 à 2000 : la loi, puis la forme combinée avec L·29 def facteurs(n): # {premier: exposant} f, p = {}, 2 while p*p <= n: while n % p == 0: f[p] = f.get(p, 0) + 1 n //= p p += 1 if n > 1: f[n] = f.get(n, 0) + 1 return f diviseurs = lambda n: [d for d in range(1, n+1) if n % d == 0] def phi(n): r = n for p in facteurs(n): r -= r // p return r from mpmath import mp, log, fsum, nstr mp.dps = 25 def Lambda(d): f = facteurs(d) return log(list(f)[0]) if len(f) == 1 else 0 print(nstr(max(abs(fsum(Lambda(d) for d in diviseurs(n)) - log(n)) for n in range(2, 2001)), 3))# 1.03e-25 print(nstr(max(abs(fsum(Lambda(d) for d in diviseurs(n))# 1.03e-25 - log(sum(phi(d) for d in diviseurs(n)))) for n in range(2, 2001)), 3))
pour \(x>0\), où \(\operatorname{Ei}(x)=\int_{-\infty}^{x}\frac{e^{t}}{t}\,dt\) (en valeur principale). La constante d'Euler de XI·2 y apparaît d'elle-même : quand \(x\to0\), \(\operatorname{Ei}(x)-\ln x\) tend vers \(\gamma\). Les factorielles au dénominateur font converger la série aussi vite que celle de \(e\) (XI·1).
# 200 tirages au hasard (écart relatif) ; puis γ retrouvé quand x → 0 from mpmath import mp, mpf, ei, log, euler, nsum, factorial, inf, nstr import random mp.dps = 25 random.seed(1) serie = lambda x: euler + log(x) + nsum(lambda k: x**k/(k*factorial(k)), [1, inf]) err = 0 for _ in range(200): x = mpf(random.uniform(0.01, 30)) err = max(err, abs(ei(x) - serie(x)) / abs(ei(x))) print(nstr(err, 3)) # 2.55e-26 x = mpf(10)**-12 # l'écart restant avec γ vaut environ x print(ei(x) - log(x)) # 0.5772156649025328606065122 print(euler) # 0.5772156649015328606065121
pour \(x>1\), où \(\operatorname{li}(x)=\int_0^{x}\frac{dt}{\ln t}\) (en valeur principale). Avec L·30, pour tout entier \(n\geq2\) : \(\operatorname{li}(n)=\operatorname{Ei}\bigl(\sum_{d\mid n}\Lambda(d)\bigr)\) ; attention, c'est bien \(n\), un entier, des deux côtés. Le logarithme intégral est la meilleure approximation simple du nombre \(\pi(x)\) de nombres premiers jusqu'à \(x\) : Gauss l'avait deviné, et c'est le théorème des nombres premiers (Hadamard et de la Vallée Poussin, 1896). Jusqu'à des hauteurs immenses, li(x) dépasse \(\pi(x)\), mais Littlewood a montré en 1914 que l'ordre s'inverse une infinité de fois.
# 200 tirages ; la forme li(n) = Ei(Σ Λ(d)) pour n = 2 à 500 ; puis π(10⁶) face à li(10⁶) def facteurs(n): # {premier: exposant} f, p = {}, 2 while p*p <= n: while n % p == 0: f[p] = f.get(p, 0) + 1 n //= p p += 1 if n > 1: f[n] = f.get(n, 0) + 1 return f diviseurs = lambda n: [d for d in range(1, n+1) if n % d == 0] def phi(n): r = n for p in facteurs(n): r -= r // p return r from mpmath import mp, mpf, li, ei, log, fsum, nstr import random mp.dps = 25 random.seed(1) err = max(abs(li(x) - ei(log(x))) / li(x) for x in (mpf(random.uniform(1.5, 1e6)) for _ in range(200))) print(nstr(err, 3)) # 9.98e-26 def Lambda(d): f = facteurs(d) return log(list(f)[0]) if len(f) == 1 else 0 print(nstr(max(abs(li(n) - ei(fsum(Lambda(d) for d in diviseurs(n)))) for n in range(2, 501)), 3))# 8.27e-24 N = 10**6 # comptage des premiers par crible crible = bytearray([1])*(N+1); crible[0] = crible[1] = 0 for p in range(2, int(N**0.5)+1): if crible[p]: crible[p*p::p] = bytearray(len(crible[p*p::p])) print(sum(crible)) # 78498 print(nstr(li(N), 8)) # 78627.549
pour tout réel \(x\) (et même tout complexe). Avec \(x=\pi\), c'est XII·1. En l'élevant à la puissance \(m\), on obtient la formule de Moivre, \((\cos x+i\sin x)^{m}=\cos mx+i\sin mx\) : \(e^{4ik\pi/n}\) est le carré de \(e^{2ik\pi/n}\). Avec \(x=\pi/4\), on trouve une racine carrée de \(i\) : \(e^{i\pi/4}=\tfrac{\sqrt2}{2}+i\,\tfrac{\sqrt2}{2}\), dont le carré vaut bien \(e^{i\pi/2}=i\).
# 200 tirages : Euler et Moivre ; puis la racine carrée de i from mpmath import mp, mpf, exp, cos, sin, pi, j, sqrt, nstr import random mp.dps = 25 random.seed(1) err = 0 for _ in range(200): x, m = mpf(random.uniform(-50, 50)), random.randint(-9, 9) err = max(err, abs(exp(j*x) - (cos(x) + j*sin(x))), abs((cos(x) + j*sin(x))**m - (cos(m*x) + j*sin(m*x)))) print(nstr(err, 3)) # 7.78e-26 r = sqrt(2)/2 + j*sqrt(2)/2 print(nstr(abs(exp(j*pi/4) - r), 3), nstr(r**2, 5))# 1.29e-26 (0.0 + 1.0j)
pour tout entier \(n\geq2\). Les \(n\) nombres \(e^{2ik\pi/n}\) sont les sommets d'un polygone régulier inscrit dans le cercle unité ; leur centre de gravité est le centre du cercle, donc leur somme est nulle. Leur produit vaut \((-1)^{n+1}\). Pour \(n=2\), c'est \(1+e^{i\pi}=0\), l'identité d'Euler (XII·1) elle-même.
# n = 2 à 60 : la somme, puis le produit from mpmath import mp, exp, pi, j, fsum, fprod, nstr mp.dps = 25 racines = lambda n: [exp(2*j*k*pi/n) for k in range(n)] print(nstr(max(abs(fsum(racines(n))) for n in range(2, 61)), 3))# 4.12e-25 print(nstr(max(abs(fprod(racines(n)) - (-1)**(n+1)) for n in range(1, 61)), 3))# 1.18e-24
chaque opération répète la précédente : la multiplication répète l'addition, la puissance répète la multiplication, la « tour » \(2\uparrow\uparrow2=2^{2}\) répète la puissance, et ainsi de suite sans fin. Appliquée à 2 et 2, chaque étage répète deux fois 2 à l'étage du dessous, et retombe donc sur 4. Aucun autre nombre ne fait ça : 3 donne 6, 9, 27, puis 7 625 597 484 987. Le seul \(x\) non nul tel que \(x+x=x\times x\) est 2. Côté logarithmes : \(\log_b a+\log_b a=\log_b a\times\log_b a\) exactement quand \(\log_b a\) vaut 0 ou 2, c'est-à-dire \(a=1\) ou \(a=b^{2}\).
# 2 ∘ 2 sur huit étages, 3 ∘ 3 sur quatre, puis les solutions de x + x = x × x def hyper(n, a, b): # étage n : 1 = +, 2 = ×, 3 = puissance, 4 = tour… if n == 1: return a + b if n == 2: return a * b if n == 3: return a ** b if b == 1: return a return hyper(n-1, a, hyper(n, a, b-1)) # l'étage n répète l'étage n − 1 print([hyper(n, 2, 2) for n in range(1, 9)])# [4, 4, 4, 4, 4, 4, 4, 4] print([hyper(n, 3, 3) for n in range(1, 5)])# [6, 9, 27, 7625597484987] print([x for x in range(0, 1001) if x + x == x * x])# [0, 2]
où \(\genfrac{[}{]}{0pt}{}{n+1}{2}\) compte les permutations de \(n+1\) objets à exactement deux cycles (Stirling de première espèce, L·24). Additionner les fractions sur le dénominateur commun \(n!\), sans simplifier, donne ces numérateurs : 1, 3, 11, 50, 274… (A000254). Une fois simplifiée, \(H_n\) = A001008/A002805 : 1, 3/2, 11/6, 25/12, 137/60, 49/20… Le rapport \(n!/\operatorname{PPCM}(1,\dots,n)\) (A025527) mesure ce que la simplification enlève. Quand \(n\) grandit, \(H_n-\ln n\) tend vers \(\gamma\) (XI·2).
# n = 1 à 30 en fractions exactes ; numérateurs sur n!, puis Hₙ simplifié from fractions import Fraction from math import factorial s = [[1]] # Stirling 1re espèce (L·24, poids n − 1) for n in range(1, 32): prev = s[-1] + [0] s.append([(prev[k-1] if k else 0) + (n-1)*prev[k] for k in range(n+1)]) H = [Fraction(0)] for n in range(1, 31): H.append(H[-1] + Fraction(1, n)) print(all(Fraction(s[n+1][2], factorial(n)) == H[n] for n in range(1, 31)))# True print([s[n+1][2] for n in range(1, 8)])# [1, 3, 11, 50, 274, 1764, 13068] print([str(H[n]) for n in range(1, 8)])# ['1', '3/2', '11/6', '25/12', '137/60', '49/20', '363/140']
Theisinger, 1915. Soit \(2^{k}\) la plus grande puissance de 2 qui ne dépasse pas \(n\) : \(\tfrac1{2^{k}}\) est le seul terme dont le dénominateur contient \(2^{k}\). En multipliant toute la somme par \(\operatorname{PPCM}(1,\dots,n)/2\), tous les termes deviennent entiers sauf celui-là, qui devient un demi-entier : la somme ne peut donc pas être entière. En prime, le numérateur de \(H_n\) est toujours impair, et son dénominateur contient exactement \(2^{k}\).
# n = 2 à 3000, fractions exactes : dénominateur divisible par 2^k et pas plus, numérateur impair from fractions import Fraction H, ok = Fraction(0), True for n in range(1, 3001): H += Fraction(1, n) k = n.bit_length() - 1 # 2^k ≤ n < 2^(k+1) d = H.denominator ok &= n == 1 or (d % 2**k == 0 and d % 2**(k+1) != 0 and H.numerator % 2 == 1) print(ok) # True
pour tout premier \(p\geq5\) (Wolstenholme, 1862) : le numérateur de \(1+\tfrac12+\dots+\tfrac1{p-1}\) est divisible par \(p^{2}\). C'est pour cela qu'on appelle aussi ces numérateurs « nombres de Wolstenholme ». Forme équivalente : \(\binom{2p-1}{p-1}\equiv1\pmod{p^{3}}\). Les premiers qui vont un cran plus loin (\(p^{4}\) pour le binomial) sont les « premiers de Wolstenholme » : on n'en connaît que deux, 16 843 et 2 124 679. Le code retrouve le premier.
# p premier de 5 à 400 ; puis recherche des premiers de Wolstenholme jusqu'à 20 000 from fractions import Fraction from math import comb est_premier = lambda p: p > 1 and all(p % q for q in range(2, int(p**0.5)+1)) premiers = [p for p in range(5, 400) if est_premier(p)] H = lambda n: sum(Fraction(1, k) for k in range(1, n+1)) print(all(H(p-1).numerator % p**2 == 0 for p in premiers))# True print(all(comb(2*p-1, p-1) % p**3 == 1 for p in premiers))# True def binom_mod(p, m): # C(2p−1, p−1) modulo m, sans le calculer en entier c = 1 for i in range(1, p): c = c * (p+i) % m * pow(i, -1, m) % m return c print([p for p in range(5, 20000) if est_premier(p) and binom_mod(p, p**4) == 1])# [16843]
pour \(-1
# 100 tirages dans ]−1, 1] ; le module ; puis le piège de la virgule flottante pour A = 10⁻¹⁵
from mpmath import mp, mpf, log, log10, nsum, inf, nstr
import math, random
mp.dps = 25
random.seed(1)
merc = lambda x: nsum(lambda k: mpf(-1)**(k+1) * x**k / k, [1, inf])
print(nstr(max(abs(merc(x) - log(1+x)) for x in (mpf(random.uniform(-0.9, 1)) for _ in range(100))), 3))# 0.0
print(1/log(10)) # le module M = log₁₀ e# 0.4342944819032518276511289
A = 1e-15 # en virgule flottante (double) :
print(math.log10(1 + A)) # naïf : 1 + A est arrondi d'abord# 4.821637332766433e-16
print(A / math.log(10)) # la formule A × M# 4.342944819032518e-16
print(math.log1p(A) / math.log(10)) # la fonction prévue pour ça# 4.3429448190325156e-16
with mp.workdps(50): # la référence : 50 chiffres,
ref = log10(1 + mpf('1e-15')) # sinon 1 + A tronque A là aussi
print(nstr(ref, 17)) # 4.3429448190325161e-16
pour tout entier \(a\) premier avec \(n\) (Euler, 1763). Idée de preuve : multiplier par \(a\) ne fait que mélanger les \(\varphi(n)\) restes premiers avec \(n\) ; le produit de ces restes est donc le même avant et après, ce qui force \(a^{\varphi(n)}\equiv1\). Pour \(n=p\) premier, c'est le petit théorème de Fermat (1640) : \(a^{p-1}\equiv1\pmod p\). La réciproque est fausse : les nombres de Carmichael, comme \(561=3\cdot11\cdot17\), vérifient \(a^{n-1}\equiv1\) pour tout \(a\) premier avec eux sans être premiers. Un contre-exemple à Lehmer (C·1) serait forcément un tel nombre, avec en plus \(\varphi(n)\mid n-1\) ; pour 561, \(\varphi=320\) ne divise pas 560.
# n < 300 et tous les a premiers avec n ; puis les nombres de Carmichael sous 10 000 from math import gcd phi = lambda n: sum(1 for k in range(1, n+1) if gcd(k, n) == 1) print(all(pow(a, phi(n), n) == 1 for n in range(2, 300) for a in range(1, n) if gcd(a, n) == 1))# True est_premier = lambda n: n > 1 and all(n % q for q in range(2, int(n**0.5)+1)) carmichael = [n for n in range(3, 10000, 2) if not est_premier(n) and all(pow(a, n-1, n) == 1 for a in range(2, n) if gcd(a, n) == 1)] print(carmichael) # [561, 1105, 1729, 2465, 2821, 6601, 8911] print(phi(561), 560 % phi(561))# 320 240
vrai dès que \(a\) et \(b\) ne sont pas tous deux négatifs, et en particulier pour \(a,b\geq0\). Pour \(a,b>0\), on en tire par exemple \(\dfrac{\sqrt a\sqrt b}{a}=\dfrac{b}{\sqrt a\sqrt b}=\dfrac{\sqrt b}{\sqrt a}\). Quand \(a\) et \(b\) sont tous deux négatifs, le signe bascule : \(\sqrt a\sqrt b=-\sqrt{ab}\). C'est ce qui démonte la fausse preuve célèbre \(-1=i\cdot i=\sqrt{-1}\sqrt{-1}=\sqrt{(-1)(-1)}=\sqrt1=1\) : l'erreur est au troisième signe égal.
# a, b ≥ 0 ; signes mélangés ; tous deux négatifs (où le signe bascule) ; puis la fausse preuve from mpmath import mp, mpf, sqrt, nstr import random mp.dps = 25 random.seed(1) tir = lambda lo, hi: mpf(random.uniform(lo, hi)) pos = max(abs(sqrt(a)*sqrt(b) - sqrt(a*b)) for a, b in ((tir(0, 9), tir(0, 9)) for _ in range(200))) mix = max(abs(sqrt(a)*sqrt(b) - sqrt(a*b)) for a, b in ((tir(-9, 0), tir(0, 9)) for _ in range(200))) neg = max(abs(sqrt(a)*sqrt(b) + sqrt(a*b)) for a, b in ((tir(-9, 0), tir(-9, 0)) for _ in range(200))) print(nstr(pos, 3), nstr(mix, 3), nstr(neg, 3))# 1.03e-25 2.07e-25 2.07e-25 print(sqrt(-1)*sqrt(-1), sqrt((-1)*(-1)))# (-1.0 + 0.0j) 1.0
on additionne de fines bandes verticales de largeur \(dx\), comme Zu découpe la boule en tranches (L·10). Deux pièges que les aide-mémoire taisent souvent. Premier piège : \(\int_a^b f(x)\,dx\) n'est l'aire sous la courbe que si \(f\geq0\) ; les parties sous l'axe comptent négativement, et \(\int_0^{2\pi}\sin x\,dx=0\) alors que l'aire coloriée vaut 4. Second piège : \(\int_a^b(f-g)\) n'est l'aire entre les courbes que si \(f\) reste au-dessus de \(g\) ; sinon, il faut couper aux points de croisement, ou prendre la valeur absolue. On peut aussi découper en bandes horizontales (\(dy\)) quand les courbes s'écrivent plus simplement en \(x=f(y)\).
# aire signée contre aire vraie, pour sin puis pour sin − cos ; puis un découpage en dy (9/2) from mpmath import mp, quad, sin, cos, pi, fabs, sqrt, nstr mp.dps = 25 print(nstr(quad(sin, [0, 2*pi]), 3)) # aire signée : les deux arches s'annulent# -1.36e-31 print(quad(lambda x: fabs(sin(x)), [0, pi, 2*pi]))# 4.0 # entre sin et cos sur [0, 2π] : ils se croisent en π/4 et 5π/4 print(nstr(quad(lambda x: sin(x) - cos(x), [0, 2*pi]), 3))# -3.76e-26 print(quad(lambda x: fabs(sin(x) - cos(x)), [0, pi/4, 5*pi/4, 2*pi]))# 5.656854249492380195206755 print(4*sqrt(2)) # 5.656854249492380195206755 # en bandes horizontales : entre x = y² et x = y + 2, de y = −1 à y = 2 print(quad(lambda y: (y + 2) - y**2, [-1, 2]))# 4.5
pour deux entiers \(a,b\geq1\), et seulement pour deux. Pour trois, la bonne formule est une inclusion–exclusion : \(\operatorname{ppcm}(a,b,c)=\dfrac{abc\,\operatorname{pgcd}(a,b,c)}{\operatorname{pgcd}(a,b)\,\operatorname{pgcd}(b,c)\,\operatorname{pgcd}(a,c)}\). La version naïve échoue dès trois nombres : pour \(1,2,\dots,n\), le PGCD vaut 1, mais \(n!\) dépasse \(\operatorname{ppcm}(1,\dots,n)\) d'un facteur 1, 1, 1, 2, 2, 12, 12, 48, 144… (A025527). Pour les fractions, on prend le PGCD des numérateurs sur le PPCM des dénominateurs (et inversement pour le PPCM) : \(\operatorname{pgcd}\bigl(\tfrac12,\tfrac13,\dots,\tfrac1n\bigr)=1/\operatorname{ppcm}(2,\dots,n)\), soit 1/A003418(n), et \(\operatorname{pgcd}\bigl(\tfrac1a,\tfrac1b\bigr)\cdot\operatorname{ppcm}\bigl(\tfrac1a,\tfrac1b\bigr)=\tfrac1{ab}\).
# 10 000 paires, 10 000 triplets, l'écart n!/ppcm, puis les fractions from math import gcd, lcm, factorial from fractions import Fraction as Fr from functools import reduce import random random.seed(1) tir = lambda: random.randint(1, 10**6) print(all(gcd(a, b)*lcm(a, b) == a*b for a, b in ((tir(), tir()) for _ in range(10000))))# True print(all(lcm(a, b, c)*gcd(a, b)*gcd(b, c)*gcd(a, c) == a*b*c*gcd(a, b, c)# True for a, b, c in ((tir(), tir(), tir()) for _ in range(10000)))) ppcm = lambda n: reduce(lcm, range(1, n+1)) print([factorial(n) // ppcm(n) for n in range(1, 11)])# [1, 1, 1, 2, 2, 12, 12, 48, 144, 1440] pgcd_fr = lambda x, y: Fr(gcd(x.numerator, y.numerator), lcm(x.denominator, y.denominator)) ppcm_fr = lambda x, y: Fr(lcm(x.numerator, y.numerator), gcd(x.denominator, y.denominator)) print(all(reduce(pgcd_fr, [Fr(1, k) for k in range(2, n+1)]) == Fr(1, ppcm(n)) for n in range(2, 40)))# True print(all(pgcd_fr(Fr(1, a), Fr(1, b))*ppcm_fr(Fr(1, a), Fr(1, b)) == Fr(1, a*b)# True for a in range(1, 60) for b in range(1, 60)))
et \(\sin\theta=\dfrac{e^{i\theta}-e^{-i\theta}}{2i}\), pour tout \(\theta\) : c'est la formule d'Euler (L·33) retournée, en additionnant ou soustrayant \(e^{i\theta}\) et \(e^{-i\theta}\). Elle transforme la trigonométrie en calcul d'exponentielles ; par exemple \(1+\cos(2\pi x)=\tfrac12e^{-2\pi ix}+1+\tfrac12e^{2\pi ix}=2\cos^{2}(\pi x)\). Les angles sont en radians : 1 radian vaut \(180/\pi\) degrés.
# 200 tirages complexes ; puis 1 + cos(2πx) sous ses trois formes from mpmath import mp, mpf, mpc, exp, cos, sin, pi, j, nstr import random mp.dps = 25 random.seed(1) err = 0 for _ in range(200): t = mpc(random.uniform(-20, 20), random.uniform(-3, 3)) # complexe compris err = max(err, abs(cos(t) - (exp(j*t) + exp(-j*t))/2), abs(sin(t) - (exp(j*t) - exp(-j*t))/(2*j))) print(nstr(err, 3)) # 2.13e-25 x = mpf('0.3') print(nstr(abs(1 + cos(2*pi*x) - (exp(-2*pi*j*x)/2 + 1 + exp(2*pi*j*x)/2)), 3),# 0.0 1.29e-26 nstr(abs(1 + cos(2*pi*x) - 2*cos(pi*x)**2), 3))
pour \(|x|\leq1\) seulement (Gregory, 1671 ; Madhava l'avait trouvée vers 1400). En degrés, avec \(x=\tan A\) : \(\tfrac{\pi}{180}A=\tan A-\tfrac13\tan^{3}A+\tfrac15\tan^{5}A-\dots\), valable pour \(-45°\leq A\leq45°\). À \(A=45°\), c'est Leibniz (I·2) ; à 60°, \(\tan A=\sqrt3>1\) et la série diverge. Leibniz est d'une lenteur désespérante (1 000 termes pour 3 chiffres). Machin (1706) a contourné le problème : \(\tfrac{\pi}{4}=4\arctan\tfrac15-\arctan\tfrac1{239}\), où chaque terme apporte environ 1,4 chiffre ; il en a tiré 100 décimales de π à la main.
# 100 tirages dans ]−1, 1[ ; la version en degrés ; Leibniz ; la divergence au-delà de 45° ; puis Machin from mpmath import mp, mpf, atan, tan, pi, sqrt, nstr import random mp.dps = 30 # chiffres de garde random.seed(1) greg = lambda x, N: sum(mpf(-1)**k * x**(2*k+1)/(2*k+1) for k in range(N)) err = max(abs(greg(x, 400) - atan(x)) for x in (mpf(random.uniform(-0.9, 0.9)) for _ in range(100))) A = mpf(30) # un angle en degrés deg = (pi*A/180, greg(tan(A*pi/180), 200)) leib = 4*greg(mpf(1), 1000) machin = 4*(4*greg(mpf(1)/5, 21) - greg(mpf(1)/239, 21)) mp.dps = 25 print(nstr(err, 3)) # 8.87e-31 print(deg[0], deg[1]) # 0.5235987755982988730771072 0.5235987755982988730771072 print(nstr(leib, 10)) # Leibniz, 1 000 termes# 3.140592654 print(nstr(greg(sqrt(3), 20), 6), nstr(greg(sqrt(3), 21), 6)) # 60° : ça diverge# -3.82045e+7 1.09095e+8 print(machin) # Machin, 21 termes# 3.141592653589793238462643 print(pi) # 3.141592653589793238462643
aire du polygone régulier à \(n\geq3\) côtés de longueur \(b\) (en degrés, \(\cot(180°/n)\)) ; son périmètre vaut \(nb\). L'aire est aussi le demi-produit du périmètre par l'apothème \(\tfrac b2\cot\tfrac\pi n\) : c'est L·15 (\(V=rS/n\)) en dimension 2. Inscrit dans un cercle de rayon \(R\), le côté vaut \(b=2R\sin\tfrac\pi n\), l'aire \(\tfrac n2R^{2}\sin\tfrac{2\pi}{n}\), et elle tend vers \(\pi R^{2}\) quand \(n\) grandit, comme le demi-périmètre tend vers \(\pi R\) (V·2).
# 100 polygones tirés au hasard, aire par les sommets ; puis un million de côtés from mpmath import mp, mpf, cos, sin, cot, pi, fsum, nstr import random mp.dps = 25 random.seed(1) def polygone(n, b): # sommets réels, puis aire par la formule du lacet R = b / (2*sin(pi/n)) P = [(R*cos(2*pi*k/n), R*sin(2*pi*k/n)) for k in range(n)] return abs(fsum(P[k][0]*P[(k+1) % n][1] - P[(k+1) % n][0]*P[k][1] for k in range(n))) / 2 err = 0 for _ in range(100): n, b = random.randint(3, 40), mpf(random.uniform(0.1, 10)) A = n*b**2*cot(pi/n)/4 err = max(err, abs(polygone(n, b) - A)/A, abs(A - (n*b)*(b/2*cot(pi/n))/2)/A) # ½ × périmètre × apothème print(nstr(err, 3)) # 4.77e-26 n = 10**6 print(nstr(n/2*sin(2*pi/n), 15), nstr(pi, 15)) # rayon 1 : l'aire tend vers π# 3.14159265356912 3.14159265358979
pour \(\operatorname{Re}s>1\) (Euler, 1737). Le mécanisme est un crible d'Ératosthène écrit en analyse : multiplier \(\zeta(s)\) par \((1-2^{-s})\) retire tous les termes pairs (II·5), puis \((1-3^{-s})\) retire les multiples de 3, et ainsi de suite ; quand tous les premiers ont criblé, il ne reste que le terme 1. En \(s=1\), la série harmonique diverge, donc le produit aussi : il y a une infinité de nombres premiers (preuve d'Euler). Plus finement, Mertens (1874) : \(\prod_{p\leq x}\frac{1}{1-1/p}\approx e^{\gamma}\ln x\), où revient la constante d'Euler (XI·2). Ses cousins : \(1/\zeta(s)=\sum\mu(n)/n^{s}\) (VIII·2) et \(-\zeta'/\zeta=\sum\Lambda(n)/n^{s}\) (L·30).
# un s complexe ; le crible après 2, 3, 5, 7 ; puis Mertens N = 10**6 crible = bytearray([1])*(N+1); crible[0] = crible[1] = 0 for q in range(2, int(N**0.5)+1): if crible[q]: crible[q*q::q] = bytearray(len(crible[q*q::q])) premiers = [q for q in range(N+1) if crible[q]] from mpmath import mp, mpf, mpc, zeta, fprod, log, euler, exp, nstr mp.dps = 20 s = mpc(3, 2) # un s complexe, Re s > 1 print(nstr(abs(fprod(1/(1 - mpf(q)**-s) for q in premiers[:20000]) - zeta(s)), 3))# 5.47e-13 c = [1]*201 # le crible sur les coefficients de ζ for q in (2, 3, 5, 7): for n in range(200, 0, -1): if n % q == 0: c[n] -= c[n // q] # multiplier par (1 − q⁻ˢ) print([n for n in range(1, 201) if c[n]][:16])# [1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67] mertens = lambda x: fprod(1/(1 - mpf(1)/q) for q in premiers if q <= x) / log(x) print([nstr(mertens(x), 6) for x in (10**2, 10**4, 10**6)]) # Mertens, en s = 1# ['1.80479', '1.78327', '1.78114'] print(nstr(exp(euler), 6)) # 1.78107
pour tout irrationnel réel \(x\). Euler a prouvé un sens (périodique implique quadratique, 1737), Lagrange l'autre (1770). Exemples : \(\sqrt2=[1;2,2,2,\dots]\), \(\sqrt3=[1;1,2,1,2,\dots]\), \(\sqrt7=[2;1,1,1,4,\dots]\), et \(\varphi=[1;1,1,\dots]\) (XIII·1). Pour \(\sqrt n\), la période finit toujours par \(2a_0\). À l'inverse, \(e\) (XIII·2) et \(\pi\) n'ont pas de période : ils ne sont pas quadratiques.
# quelques racines ; les 283 non-carrés jusqu'à 300 reconstruits ; puis la période qui se ferme sur 2a₀ from math import isqrt from mpmath import mp, mpf, sqrt, nstr mp.dps = 30 def fc_racine(n, nb): # algorithme exact en entiers pour √n a0 = isqrt(n); m, d, a, t = 0, 1, a0, [a0] while len(t) < nb: m = d*a - m; d = (n - m*m)//d; a = (a0 + m)//d; t.append(a) return t def valeur(t): v = mpf(t[-1]) for a in reversed(t[:-1]): v = a + 1/v return v print({n: fc_racine(n, 9) for n in (2, 3, 7, 13)})# {2: [1, 2, 2, 2, 2, 2, 2, 2, 2], 3: [1, 1, 2, 1, 2, 1, 2, 1, 2], 7: [2, 1, 1, 1, 4, 1, 1, 1, 4], 13: [3, 1, 1, 1, 1, 6, 1, 1, 1]} non_carres = [n for n in range(2, 301) if isqrt(n)**2 != n] print(nstr(max(abs(valeur(fc_racine(n, 60)) - sqrt(n)) for n in non_carres), 3))# 1.58e-30 def periode_finit_par_2a0(n): t = fc_racine(n, 2*n + 5); a0 = t[0] k = t.index(2*a0) # la première fois qu'apparaît 2a₀… return t[1:k+1] == t[k+1:2*k+1] # … ferme bien une période qui se répète print(all(periode_finit_par_2a0(n) for n in non_carres))# True
si \(x=[\,\overline{a_0;a_1,\dots,a_{n-1}}\,]\) est purement périodique, où \(\bar x\) est l'autre racine de son équation du second degré (son « conjugué »). C'est le premier article publié par Galois, en 1829, à 17 ans. Il prouve aussi qu'une fraction continue est purement périodique exactement quand \(x>1\) et \(-1<\bar x<0\). Pour \(\varphi=[\,\overline{1}\,]\), la période se lit pareil à l'envers : \(-1/\bar\varphi=\varphi\), c'est-à-dire \(\bar\varphi=-1/\varphi\approx-0{,}618\).
É. Galois, « Démonstration d'un théorème sur les fractions continues périodiques », Annales de mathématiques pures et appliquées (Gergonne), 1829.
# 200 périodes tirées au hasard ; puis l'exemple (1, 2, 3, 4) from mpmath import mp, mpf, sqrt, nstr import random mp.dps = 40 random.seed(1) def racines(per): # x = [per; x] ⇒ q x² + (q' − p) x − p' = 0 p, pp, q, qq = 1, 0, 0, 1 for a in per: p, pp, q, qq = a*p + pp, p, a*q + qq, q A, B, C = q, qq - p, -pp D = sqrt(B*B - 4*A*C) return (-B + D)/(2*A), (-B - D)/(2*A) # x, et son conjugué def valeur(t): v = mpf(t[-1]) for a in reversed(t[:-1]): v = a + 1/v return v ok = True for _ in range(200): per = [random.randint(1, 9) for _ in range(random.randint(1, 5))] x, xb = racines(per) ok &= -1 < xb < 0 and abs(-1/xb - valeur(per[::-1]*(200//len(per)))) < mpf(10)**-25 print(ok) # True x, xb = racines([1, 2, 3, 4]) print(nstr(xb, 12), nstr(-1/xb, 12), nstr(valeur([4, 3, 2, 1]*40), 12))# -0.232666399786 4.29799919936 4.29799919936
où les \(B_j^{+}\) sont les nombres de Bernoulli avec la convention \(B_1^{+}=+\tfrac12\) (les autres sont inchangés : 1, 1/6, 0, −1/30…). Attention au piège : avec la convention \(B_1=-\tfrac12\), la même formule s'arrête à \(n-1\) (pour \(p=1\), \(n=3\), elle donne 3 au lieu de 6). On peut aussi garder \(B_1=-\tfrac12\) et écrire \((-1)^{j}B_j\), ou passer par les polynômes de Bernoulli : \(\sum_{k=0}^{n}k^{p}=\bigl(B_{p+1}(n+1)-B_{p+1}(0)\bigr)/(p+1)\). Faulhaber (1631) calculait ces sommes à la main jusqu'à \(p=17\) ; Jacques Bernoulli (1713) se vantait d'obtenir \(1^{10}+\dots+1000^{10}\) en moins d'un quart d'heure. Cas \(p=3\), le théorème de Nicomaque : \(1^{3}+2^{3}+\dots+n^{3}=(1+2+\dots+n)^{2}\), une pile de cubes qui forme un carré.
# p = 0 à 11, n = 1 à 29 en fractions exactes ; le piège de B₁ ; Bernoulli 1713 ; Nicomaque from fractions import Fraction as Fr from math import comb B = [Fr(1)] # Bernoulli, convention B₁ = −½ for m in range(1, 16): B.append(-sum(comb(m+1, j)*B[j] for j in range(m)) / (m+1)) Bplus = B[:]; Bplus[1] = Fr(1, 2) faul = lambda n, p, b: sum(comb(p+1, j)*b[j]*n**(p+1-j) for j in range(p+1)) / Fr(p+1) somme = lambda n, p: sum(k**p for k in range(1, n+1)) print(all(faul(n, p, Bplus) == somme(n, p) for p in range(12) for n in range(1, 30)))# True print(faul(3, 1, B), faul(3, 1, Bplus), somme(3, 1)) # le piège de la convention# 3 6 6 print(faul(1000, 10, Bplus)) # le calcul de Jacques Bernoulli# 91409924241424243424241924242500 print(all(somme(n, 3) == somme(n, 1)**2 for n in range(1, 500))) # Nicomaque# True
pour \(|x|<2\pi\) : la fonction explose en \(x=\pm2\pi i\), et c'est de là que vient le \((2\pi)^{2k}\) de L·25. Pourquoi les \(B_n\) impairs sont nuls (sauf \(B_1\)) : en notant \(B(x)\) cette fonction, \(B(-x)=\dfrac{xe^{x}}{e^{x}-1}=B(x)+x\). Dans \(B(-x)-B(x)\), les termes pairs s'annulent et les impairs se doublent, alors qu'un seul terme doit survivre, \(x\) : donc \(B_1=-\tfrac12\) et \(B_3=B_5=\dots=0\). Liens avec ζ : \(B_n=-n\,\zeta(1-n)\) pour \(n\geq2\) (pour \(n=1\), on obtient \(+\tfrac12\), l'autre convention), et \(\zeta(2k)\) en découle (L·25).
# les Bₙ tirés de la série elle-même ; B(−x) = B(x) + x ; le rayon 2π ; puis Bₙ = −n ζ(1−n) from fractions import Fraction as Fr from math import factorial from mpmath import mp, mpf, exp, bernoulli, zeta, nstr import random mp.dps = 30 N = 16 a = [Fr(1, factorial(k+1)) for k in range(N)] # (eˣ − 1)/x = Σ xᵏ/(k+1)! c = [Fr(1)] for n in range(1, N): c.append(-sum(a[k]*c[n-k] for k in range(1, n+1))) # on inverse la série print([str(c[n]*factorial(n)) for n in range(9)]) # Bₙ = n! × coefficient# ['1', '-1/2', '1/6', '0', '-1/30', '0', '1/42', '0', '-1/30'] random.seed(1) Bx = lambda x: x/(exp(x) - 1) print(nstr(max(abs(Bx(-x) - Bx(x) - x) for x in (mpf(random.uniform(-5, 5)) for _ in range(100))), 3))# 2.47e-30 serie = lambda x, M: sum(bernoulli(n)*x**n/factorial(n) for n in range(M)) print(nstr(serie(mpf(5), 200) - Bx(mpf(5)), 3), nstr(serie(mpf(7), 200), 3)) # rayon 2π ≈ 6,28# 1.76e-20 2.16e+9 print(all(abs(-n*zeta(1-n) - bernoulli(n)) < mpf(10)**-20 for n in range(2, 30)), -1*zeta(0))# True 0.5
pour tout \(b\geq1\) : la somme des racines \(b\)-ièmes primitives de l'unité vaut la fonction de Möbius (VIII·2). Idée : toutes les racines \(b\)-ièmes se compensent (L·34), et on les regroupe selon leur ordre exact. En sommant sur toutes les fractions réduites \(a/b\) de dénominateur \(b\leq n\) (il y en a A002088(n), numérateurs A038566, dénominateurs A038567), on obtient la fonction de Mertens \(M(n)=\mu(1)+\dots+\mu(n)\) : 1, 0, −1, −1, −2, −1, −2… (A002321). L'hypothèse de Riemann (C·2) équivaut à ce que \(M(n)\) ne dépasse pas \(n^{1/2+\varepsilon}\). Mertens avait conjecturé \(|M(n)|<\sqrt n\) pour tout \(n>1\) : c'est faux, Odlyzko et te Riele l'ont réfuté en 1985, sans qu'on connaisse de contre-exemple explicite.
# b = 1 à 300 ; la somme de Farey contre M(n) ; puis |M(n)| < √n jusqu'à 10⁵ (vrai si loin, faux un jour) import cmath from math import gcd, pi, isqrt N = 10**5 mu = [1]*(N+1); premier = [True]*(N+1) # Möbius par crible for q in range(2, N+1): if premier[q]: for m in range(q, N+1, q): if m > q: premier[m] = False mu[m] = -mu[m] for m in range(q*q, N+1, q*q): mu[m] = 0 prim = lambda b: sum(cmath.exp(2j*pi*a/b) for a in range(1, b+1) if gcd(a, b) == 1) print(all(abs(prim(b) - mu[b]) < 1e-9 for b in range(1, 301)))# True farey = lambda n: sum(prim(b) for b in range(1, n+1)) M = [0] for k in range(1, N+1): M.append(M[-1] + mu[k]) print([round(farey(n).real) for n in range(1, 13)])# [1, 0, -1, -1, -2, -1, -2, -2, -2, -1, -2, -2] print(M[1:13]) # [1, 0, -1, -1, -2, -1, -2, -2, -2, -1, -2, -2] print(max(abs(M[n]) / n**0.5 for n in range(2, N+1)) < 1)# True
où \(d\) compte les diviseurs (A056793 est \(d(\operatorname{ppcm}(1,\dots,n))\)). Le PPCM ne change qu'en passant une puissance de premier : un nouveau premier double le nombre de diviseurs, et \(p^{k}\) fait passer l'exposant de \(p\) de \(k-1\) à \(k\), d'où le facteur \((k+1)/k\). Ce facteur vaut au plus \(\tfrac32\), atteint exactement aux carrés de premiers : 4, 9, 25, 49… Les \(n+1\) où rien ne bouge sont les nombres qui ne sont pas des puissances de premiers (A024619) : 6, 10, 12, 14, 15…
# n + 1 jusqu'à 300, PPCM factorisés pour de vrai ; puis les sauts non entiers from math import lcm from fractions import Fraction est_premier = lambda m: m > 1 and all(m % q for q in range(2, int(m**0.5)+1)) def d(m): # diviseurs, par vraie factorisation de m c, q = 1, 2 while m > 1: e = 0 while m % q == 0: m //= q; e += 1 c *= e + 1; q += 1 return c def puissance(m): # (p, k) si m = p^k, sinon None for q in range(2, m+1): if m % q == 0: k = 0 while m % q == 0: m //= q; k += 1 return (q, k) if m == 1 else None L = [1] for n in range(1, 301): L.append(lcm(L[-1], n)) ok = True for n in range(1, 300): r, pk = Fraction(d(L[n+1]), d(L[n])), puissance(n+1) ok &= r == (1 if pk is None else 2 if pk[1] == 1 else Fraction(pk[1]+1, pk[1])) print(ok) # True print([(n+1, str(Fraction(d(L[n+1]), d(L[n])))) for n in range(1, 60)# [(4, '3/2'), (8, '4/3'), (9, '3/2'), (16, '5/4'), (25, '3/2'), (27, '4/3'), (32, '6/5'), (49, '3/2')] if d(L[n+1]) != d(L[n]) and not est_premier(n+1)])
pour tout polygone simple (sans trou, sans croisement) dont les sommets sont sur une grille de pas 1 : \(B\) compte les points de la grille sur le bord, \(I\) ceux à l'intérieur (Pick, 1899). Exemple : un polygone à 12 sommets avec \(B=20\) et \(I=8\) a pour aire \(10+8-1=17\). Avec \(h\) trous, la formule devient \(\mathcal A=I+\tfrac B2-1+h\). En dimension 3, aucun équivalent n'existe : le tétraèdre de Reeve (1957), de sommets \((0,0,0),(1,0,0),(0,1,0),(1,1,r)\), n'a aucun point de la grille hormis ses quatre sommets, quel que soit \(r\), alors que son volume \(r/6\) grandit sans fin. Encore la dimension 2 qui se laisse faire, et la 3 qui se dérobe.
# ton polygone recompté ; 300 triangles au hasard ; puis les tétraèdres de Reeve (points, volume) from math import gcd from fractions import Fraction as Fr import random def pick(P): # aire (lacet), bord et intérieur comptés pour de vrai n = len(P) A2 = abs(sum(P[i][0]*P[(i+1)%n][1] - P[(i+1)%n][0]*P[i][1] for i in range(n))) B = sum(gcd(abs(P[(i+1)%n][0]-P[i][0]), abs(P[(i+1)%n][1]-P[i][1])) for i in range(n)) def bord(x, y): return any((b[0]-a[0])*(y-a[1]) == (b[1]-a[1])*(x-a[0]) and min(a[0],b[0]) <= x <= max(a[0],b[0]) and min(a[1],b[1]) <= y <= max(a[1],b[1]) for a, b in zip(P, P[1:]+P[:1])) def dedans(x, y): c = False for (x1, y1), (x2, y2) in zip(P, P[1:]+P[:1]): if (y1 > y) != (y2 > y) and x < x1 + Fr(y-y1)*(x2-x1)/(y2-y1): c = not c return c xs, ys = [q[0] for q in P], [q[1] for q in P] I = sum(1 for x in range(min(xs), max(xs)+1) for y in range(min(ys), max(ys)+1) if not bord(x, y) and dedans(x, y)) return Fr(A2, 2), B, I exemple = [(1,2),(3,2),(4,4),(6,2),(6,7),(5,6),(3,6),(2,7),(1,7),(1,5),(3,4),(1,3)] A, B, I = pick(exemple) print(A, B, I, Fr(B, 2) + I - 1)# 17 20 8 17 random.seed(1) ok = True for _ in range(300): # triangles de la grille tirés au hasard T = [(random.randint(-9, 9), random.randint(-9, 9)) for _ in range(3)] if (T[1][0]-T[0][0])*(T[2][1]-T[0][1]) == (T[2][0]-T[0][0])*(T[1][1]-T[0][1]): continue A, B, I = pick(T); ok &= A == Fr(B, 2) + I - 1 print(ok) # True def points_reeve(r): # points de la grille dans le tétraèdre de Reeve S = [(0,0,0), (1,0,0), (0,1,0), (1,1,r)] def dedans(p): # coordonnées barycentriques exactes : p = S0 + u(S1−S0) + v(S2−S0) + w(S3−S0) x, y, z = p w = Fr(z, r); u = x - w; v = y - w return min(u, v, w, 1 - u - v - w) >= 0 return sum(dedans((x, y, z)) for x in (0, 1) for y in (0, 1) for z in range(r+1)) print([(r, points_reeve(r), str(Fr(r, 6))) for r in (1, 2, 5, 20)])# [(1, 4, '1/6'), (2, 4, '1/3'), (5, 4, '5/6'), (20, 4, '10/3')]
avec \(A=\sqrt{a^{2}+b^{2}}\) et \(\theta=\operatorname{atan2}(b,a)\), c'est-à-dire \(\tan\theta=b/a\) (partie imaginaire sur partie réelle) dans le bon quadrant. Piège : \(\arctan(b/a)\) confond \(1+i\) et \(-1-i\), dont les angles diffèrent de π ; c'est la fonction Math.atan2 des langages de programmation qui tranche. En polaire, multiplier revient à multiplier les modules et additionner les angles (L·33). En physique, c'est le phaseur : un point qui tourne sur un cercle, \(A\,e^{i(\omega t+\varphi)}\), dont l'ombre sur l'axe réel est l'oscillation \(A\cos(\omega t+\varphi)\).
# 500 tirages avec atan2, et combien de fois arctan(b/a) se trompe ; le produit ; puis le phaseur from mpmath import mp, mpf, mpc, sqrt, atan2, atan, exp, cos, pi, j, nstr import random mp.dps = 25 random.seed(1) err, piege = 0, 0 for _ in range(500): a, b = mpf(random.uniform(-9, 9)), mpf(random.uniform(-9, 9)) A, t = sqrt(a*a + b*b), atan2(b, a) err = max(err, abs(A*exp(j*t) - (a + j*b))) piege += abs(A*exp(j*atan(b/a)) - (a + j*b)) > 1e-10 # arctan se trompe quand a < 0 print(nstr(err, 3), piege) # 3.27e-25 234 z1, z2 = mpc(1, 2), mpc(-3, 1) # produit : modules multipliés, angles ajoutés print(nstr(abs(abs(z1*z2) - abs(z1)*abs(z2)), 3),# 0.0 1.29e-26 nstr(abs(exp(j*(atan2(z1.imag, z1.real) + atan2(z2.imag, z2.real))) - z1*z2/abs(z1*z2)), 3)) A, w, phi, t = mpf(2), mpf(3), mpf('0.4'), mpf('1.7') # le phaseur print(nstr(abs((A*exp(j*(w*t + phi))).real - A*cos(w*t + phi)), 3))# 0.0
pour tout \(\lambda>0\) et tout entier \(k\geq0\) (Poisson, 1837). L'image du tissu : un long tissu présente en moyenne \(\lambda\) déchirures par mètre, indépendantes les unes des autres. On coupe le mètre en \(n\) morceaux si petits que chacun a au plus une déchirure, avec probabilité \(\lambda/n\) : c'est un pile ou face très déséquilibré, répété \(n\) fois, et la binomiale tend vers la loi de Poisson. Les probabilités somment à 1 grâce à la série de \(e\) (XI·1), et la moyenne vaut \(\lambda\). Même loi, ailleurs : le nombre de points fixes d'une permutation tirée au hasard suit Poisson(1), d'où la probabilité \(1/e\) des dérangements (XI·1). Et l'exemple historique de Bortkiewicz (1898) : les morts par ruade de cheval dans la cavalerie prussienne.
# binomiale contre Poisson pour n = 10, 100, 10⁴ ; somme et moyenne ; puis les points fixes de 8! permutations from itertools import permutations from math import comb, factorial as fact, exp from mpmath import mp, mpf, nsum, inf, factorial, e, nstr mp.dps = 25 lam = 2.5 poisson = lambda k, l=lam: exp(-l) * l**k / fact(k) binom = lambda k, n: comb(n, k) * (lam/n)**k * (1 - lam/n)**(n-k) print([f'{max(abs(binom(k, n) - poisson(k)) for k in range(8)):.1e}' for n in (10, 100, 10**4)])# ['3.7e-02', '3.0e-03', '2.9e-05'] L = mpf(lam) print(nstr(nsum(lambda k: e**-L * L**k / factorial(k), [0, inf]), 20),# 1.0 2.5 nstr(nsum(lambda k: k * e**-L * L**k / factorial(k), [0, inf]), 20)) n = 8 # points fixes des 8! permutations fixes = [sum(1 for i in range(n) if p[i] == i) for p in permutations(range(n))] print([round(fixes.count(k)/fact(n), 4) for k in range(4)])# [0.3679, 0.3679, 0.184, 0.0611] print([round(poisson(k, 1), 4) for k in range(4)])# [0.3679, 0.3679, 0.1839, 0.0613]
pour tout entier \(n\geq0\), avec \(\binom nk=C_n^{k}=\dfrac{n!}{k!\,(n-k)!}\) (Pascal, L·24). Avec \(y\) changé en \(-y\) : \((x-y)^{n}=\sum\binom nk x^{n-k}(-y)^{k}\). Avec \(x=y=1\), la ligne \(n\) du triangle somme à \(2^{n}\) ; avec \(x=1,\ y=-1\), la somme alternée est nulle. Newton (1665) l'a étendu à tout exposant \(\alpha\), réel ou complexe, pour \(|x|<1\) : \((1+x)^{\alpha}=\sum_{k\geq0}\binom\alpha k x^{k}\), où \(\binom\alpha k=\alpha(\alpha-1)\cdots(\alpha-k+1)/k!\) ; la série ne s'arrête alors plus. Exemple : \(1/\sqrt{1-4x}=\sum_k\binom{2k}{k}x^{k}\), où reviennent les binomiaux centraux de X·3.
# 500 cas exacts ; les sommes de lignes ; 100 exposants réels ; puis les binomiaux centraux from math import comb from mpmath import mp, mpf, binomial, nsum, inf, sqrt, nstr import random mp.dps = 25 random.seed(1) ok = True for _ in range(500): # entiers exacts, n jusqu'à 30 x, y, n = random.randint(-50, 50), random.randint(-50, 50), random.randint(0, 30) ok &= (x + y)**n == sum(comb(n, k) * x**(n-k) * y**k for k in range(n+1)) print(ok) # True print([sum(comb(n, k) for k in range(n+1)) for n in range(8)])# [1, 2, 4, 8, 16, 32, 64, 128] err = 0 for _ in range(100): # exposant quelconque (Newton 1665) a, x = mpf(random.uniform(-5, 5)), mpf(random.uniform(-0.8, 0.8)) err = max(err, abs(nsum(lambda k: binomial(a, k) * x**k, [0, inf]) - (1 + x)**a)) print(nstr(err, 3)) # 0 x = mpf('0.1') print(nstr(nsum(lambda k: binomial(2*k, k) * x**k, [0, inf]) - 1/sqrt(1 - 4*x), 3))# -2.58e-26
pour \(a,b\) entiers et \(f\) de classe \(C^{k+1}\) sur \([a,b]\) (Euler 1732, Maclaurin 1742), avec la convention \(B_1=-\tfrac12\) et le reste \(R=\dfrac{(-1)^{k}}{(k+1)!}\displaystyle\int_a^b\tilde B_{k+1}(t)\,f^{(k+1)}(t)\,dt\), où \(\tilde B_{k+1}(t)=B_{k+1}(t-\lfloor t\rfloor)\) est le polynôme de Bernoulli répété périodiquement (L·59). C'est le pont entre sommes et intégrales ; les premiers termes sont \(\tfrac12(f(b)-f(a))+\tfrac1{12}(f'(b)-f'(a))-\dots\) Trois usages dans Idem : Stirling (X·5), l'écart harmonique \(H_n-\ln n\approx\gamma+\tfrac1{2n}-\tfrac1{12n^{2}}+\tfrac1{120n^{4}}\) (XI·2), et Faulhaber (L·50), où la formule devient exacte car les dérivées d'un polynôme finissent par s'annuler. En général, la série infinie diverge : on s'en sert comme d'un développement asymptotique, arrêté au bon moment.
# la formule avec son reste, pour k = 1 à 4 ; puis Hₙ − ln n − γ à n = 10 contre son développement from mpmath import mp, mpf, fsum, quad, diff, bernoulli, bernpoly, factorial, floor, linspace from mpmath import harmonic, log, euler, nstr mp.dps = 30 f = lambda t: 1/(1 + t**2) a, b = 0, 6 def ecart(k): # somme − intégrale − (termes + reste) somme = fsum(f(n) for n in range(a+1, b+1)) - quad(f, [a, b]) termes = fsum((-1)**(r+1) * bernoulli(r+1) / factorial(r+1) * (diff(f, b, r) - diff(f, a, r)) for r in range(k+1)) R = (-1)**k / factorial(k+1) * quad(lambda t: bernpoly(k+1, t - floor(t)) * diff(f, t, k+1), linspace(a, b, b - a + 1)) return somme - termes - R # à la précision près des dérivées numériques (≈ 10⁻¹⁷) print([nstr(ecart(k), 3) for k in (1, 2, 3, 4)])# ['1.95e-17', '1.95e-17', '1.95e-17', '1.95e-17'] n = mpf(10) # l'écart harmonique, développé print(nstr(harmonic(n) - log(n) - euler, 15))# 0.0491674960726754 print(nstr(1/(2*n) - 1/(12*n**2) + 1/(120*n**4), 15))# 0.0491675
pour \(0\leq x\leq1\) et \(n\geq2\) : le \(k\)-ième coefficient de Fourier de \(B_n\), répété périodiquement, vaut \(c_k=-\dfrac{n!}{(2i\pi k)^{n}}\). C'est la preuve de la formule d'Euler pour \(\zeta(2p)\) (L·25) : en \(x=0\) et \(n=2p\), les termes \(k\) et \(-k\) s'additionnent, et il vient \(B_{2p}=(-1)^{p+1}\dfrac{2\,(2p)!}{(2\pi)^{2p}}\zeta(2p)\), soit \(\displaystyle\sum_{n\geq1}\frac1{n^{2p}}=\frac{|B_{2p}|\,2^{2p-1}\pi^{2p}}{(2p)!}\). Pour \(n\) impair, la même série en \(x=0\) donne 0, ce qui ne dit rien de \(\zeta(3),\zeta(5),\dots\), et c'est pourquoi on ne leur connaît aucune forme close.
# coefficients de Fourier pour n = 2 à 6 ; la série qui reconstruit B₄ ; puis ζ(2p) retrouvé from mpmath import mp, mpf, quad, bernpoly, bernoulli, exp, pi, j, factorial, zeta, fsum, nstr mp.dps = 25 coef = lambda n, k: quad(lambda x: bernpoly(n, x) * exp(-2*j*pi*k*x), [0, 1]) print(nstr(max(abs(coef(n, k) + factorial(n)/(2*j*pi*k)**n) for n in range(2, 7) for k in (1, 2, 3, -2)), 3))# 8.08e-28 x, n = mpf('0.3'), 4 # la série reconstruit B₄(0,3) serie = -factorial(n)/(2*j*pi)**n * fsum(exp(2*j*pi*k*x)/mpf(k)**n for k in range(-4000, 4001) if k) print(nstr(serie.real, 12), nstr(bernpoly(n, x), 12))# 0.0107666666667 0.0107666666667 print([nstr(abs(bernoulli(2*p)) * 2**(2*p-1) * pi**(2*p) / factorial(2*p) - zeta(2*p), 3) for p in (1, 2, 3, 10)])# ['5.17e-26', '2.58e-26', '5.17e-26', '1.03e-25']
avec \(B_n(x)=\sum_{k=0}^{n}\binom nk B_k\,x^{\,n-k}\) (attention à l'exposant : \(n-k\), et non \(n-1\)), et \(B_n(0)=B_n\) (L·51). Trois propriétés font tout le travail. La différence : \(B_n(x+1)-B_n(x)=n\,x^{n-1}\), qui télescope les sommes de puissances, d'où Faulhaber : \(\sum_{k=0}^{m}k^{p}=\bigl(B_{p+1}(m+1)-B_{p+1}(0)\bigr)/(p+1)\) (L·50). Le miroir : \(B_n(1-x)=(-1)^{n}B_n(x)\), d'où les \(B_n\) impairs nuls. La dérivée : \(B_n'=n\,B_{n-1}\), avec \(\int_0^1B_n=0\) pour \(n\geq1\). Répétés périodiquement, ils ont la série de Fourier de L·59 et servent de noyau au reste d'Euler–Maclaurin (L·58).
| n | \(B_n(x)\) |
|---|---|
| 0 | \(1\) |
| 1 | \(x-\tfrac12\) |
| 2 | \(x^{2}-x+\tfrac16\) |
| 3 | \(x^{3}-\tfrac32x^{2}+\tfrac12x\) |
| 4 | \(x^{4}-2x^{3}+x^{2}-\tfrac1{30}\) |
| 5 | \(x^{5}-\tfrac52x^{4}+\tfrac53x^{3}-\tfrac16x\) |
| 6 | \(x^{6}-3x^{5}+\tfrac52x^{4}-\tfrac12x^{2}+\tfrac1{42}\) |
# en fractions exactes : la différence, le miroir, l'intégrale nulle ; puis la fonction génératrice from fractions import Fraction as Fr from math import comb B = [Fr(1)] # nombres de Bernoulli, B₁ = −½ for m in range(1, 101): B.append(-sum(comb(m+1, j)*B[j] for j in range(m)) / (m+1)) from mpmath import mp, mpf, bernpoly, exp, factorial, nstr Bx = lambda n, x: sum(comb(n, k) * B[k] * x**(n-k) for k in range(n+1)) xs = [Fr(i, 7) for i in range(-10, 30)] # 40 points exacts : assez pour tout degré ≤ 39 print(all(Bx(n, x+1) - Bx(n, x) == n * x**(n-1) for n in range(1, 30) for x in xs))# True print(all(Bx(n, 1-x) == (-1)**n * Bx(n, x) for n in range(30) for x in xs))# True integrale = lambda n: sum(comb(n, k) * B[k] / (n-k+1) for k in range(n+1)) # ∫₀¹ Bₙ print(all(integrale(n) == 0 for n in range(1, 30)))# True mp.dps = 25 t, x = mpf('1.3'), mpf('0.4') print(nstr(sum(bernpoly(n, x) * t**n / factorial(n) for n in range(80)) - t*exp(x*t)/(exp(t) - 1), 3))# 0.0
pour tout \(n\geq1\) (von Staudt et Clausen, indépendamment, 1840). Conséquence : le dénominateur de \(B_{2n}\) est exactement le produit des premiers \(p\) tels que \(p-1\) divise \(2n\), toujours multiple de 6 ; par exemple \(B_{12}=-691/2730\), et \(2730=2\cdot3\cdot5\cdot7\cdot13\). Le théorème ne vaut que pour les indices pairs : pour \(n\) impair \(\geq3\), \(B_n=0\) et il resterait \(\tfrac12\). Le numérateur 691 est célèbre : c'est le même 691 que dans la série thêta du réseau de Leech (L·23), et Ramanujan a prouvé \(\tau(n)\equiv\sigma_{11}(n)\pmod{691}\), où \(\tau\) sont les coefficients de \(\Delta\).
# indices pairs jusqu'à 100 : l'entier, puis les dénominateurs ; B₁₂ ; l'indice 3 ; puis la congruence de Ramanujan from fractions import Fraction as Fr from math import comb B = [Fr(1)] # nombres de Bernoulli, B₁ = −½ for m in range(1, 101): B.append(-sum(comb(m+1, j)*B[j] for j in range(m)) / (m+1)) est_premier = lambda q: q > 1 and all(q % r for r in range(2, int(q**0.5)+1)) somme = lambda n: B[n] + sum(Fr(1, q) for q in range(2, n+2) if est_premier(q) and n % (q-1) == 0) print(all(somme(n).denominator == 1 for n in range(2, 101, 2)))# True from math import prod print(all(B[n].denominator == prod(q for q in range(2, n+2) if est_premier(q) and n % (q-1) == 0)# True for n in range(2, 101, 2))) print(B[12], somme(3)) # 691/2730, puis l'échec aux indices impairs# -691/2730 1/2 N = 40 # Ramanujan : τ(n) ≡ σ₁₁(n) mod 691 P = [1] + [0]*N for n in range(1, N+1): for _ in range(24): for k in range(N, n-1, -1): P[k] -= P[k-n] tau = lambda n: P[n-1] sigma11 = lambda n: sum(d**11 for d in range(1, n+1) if n % d == 0) print(all((tau(n) - sigma11(n)) % 691 == 0 for n in range(1, N+1)))# True
pour tout premier \(p\) et tout entier \(n\) : c'est le petit théorème de Fermat (L·40) déguisé, puisque \(n^{p-1}\equiv1\) si \(p\nmid n\), et \(\equiv0\) sinon. Avec \(p=5\) : \((4n^{4}+1)\bmod5\) donne 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1… L'hypothèse « \(p\) premier » est indispensable, et même caractéristique : pour tout composé \(m\), un diviseur propre \(d\) de \(m\) fait échouer la formule, puisque \(d^{m-1}\) ne peut pas valoir 1 modulo \(m\). Avec 9, on obtient 1, 0, 6, 1, 3, 3, 1, 6, 0… La formule est donc aussi un test de primalité caché, très lent.
# p = 5 ; l'échec avec 9 ; puis les m jusqu'à 300 où la formule marche : exactement les premiers est_premier = lambda q: q > 1 and all(q % r for r in range(2, int(q**0.5)+1)) indic = lambda m, n: ((m-1) * pow(n, m-1, m) + 1) % m print([indic(5, n) for n in range(11)])# [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1] print([indic(9, n) for n in range(9)])# [1, 0, 6, 1, 3, 3, 1, 6, 0] marche = [m for m in range(2, 300) if all(indic(m, n) == (1 if n % m == 0 else 0) for n in range(3*m))] print(marche == [m for m in range(2, 300) if est_premier(m)])# True
pour tous entiers \(n,m\geq1\) : la somme des produits de \(m\) entiers consécutifs est le produit de \(m+1\) consécutifs, divisé par \(m+1\). Avec \(m=1\), \(1+2+\dots+n=n(n+1)/2\) ; avec \(m=2\), \(1\cdot2+2\cdot3+\dots+n(n+1)=\dfrac{n(n+1)(n+2)}{3}=\dfrac n3\,(n^{2}+3n+2)\). Attention à ne pas écrire \(+1\) à la fin : pour \(n=1\), on obtiendrait 5/3 au lieu de 2. Derrière la formule se cache l'identité « de la crosse de hockey » du triangle de Pascal, \(\sum_{k}\binom{k+m-1}{m}=\binom{n+m}{m+1}\), puisque \(k(k+1)\cdots(k+m-1)=m!\binom{k+m-1}{m}\).
# m = 1 à 6 et n = 1 à 50 en entiers exacts ; le piège du « +1 » ; puis la crosse de hockey from math import comb, prod from fractions import Fraction as Fr print(all(sum(prod(range(k, k+m)) for k in range(1, n+1)) * (m+1) == prod(range(n, n+m+1))# True for m in range(1, 7) for n in range(1, 51))) print([sum(k*(k+1) for k in range(1, n+1)) for n in range(1, 6)])# [2, 8, 20, 40, 70] print([str(Fr(n, 3)*(n*n + 3*n + 1)) for n in range(1, 6)]) # avec « +1 » : faux# ['5/3', '22/3', '19', '116/3', '205/3'] print(all(sum(comb(k+m-1, m) for k in range(1, n+1)) == comb(n+m, m+1)# True for m in range(1, 7) for n in range(1, 51))) # crosse de hockey
pour tout entier \(m\geq2\) (énoncé par Waring en 1770, qui l'attribue à son élève Wilson, et démontré par Lagrange en 1771). Pour \(m\) composé, \((m-1)!\) est même divisible par \(m\), sauf pour \(m=4\) où il reste 2. Une conséquence : la somme \(1+2+\dots+n=n(n+1)/2\) divise \(n!\) exactement quand \(n+1\) n'est pas premier (pour \(n\geq2\)), car \(n!\big/\tfrac{n(n+1)}{2}=\dfrac{2\,(n-1)!}{n+1}\). Ce rapport est pair, sauf pour \(n=1\) et \(n=3\), où il vaut 1. Comme l'indicatrice de Fermat (L·62), c'est un test de primalité exact mais très lent.
# Wilson jusqu'à 2000 ; les restes pour les composés ; puis ta divisibilité et ses rapports impairs from math import factorial est_premier = lambda m: m > 1 and all(m % q for q in range(2, int(m**0.5)+1)) def fact_mod(m): # (m−1)! modulo m, sans calculer l'entier r = 1 for k in range(2, m): r = r * k % m return r print(all(((fact_mod(m) + 1) % m == 0) == est_premier(m) for m in range(2, 2001)))# True print({m: fact_mod(m) for m in range(4, 26) if not est_premier(m)})# {4: 2, 6: 0, 8: 0, 9: 0, 10: 0, 12: 0, 14: 0, 15: 0, 16: 0, 18: 0, 20: 0, 21: 0, 22: 0, 24: 0, 25: 0} T = lambda n: n*(n+1)//2 print(all((factorial(n) % T(n) == 0) == (not est_premier(n+1)) for n in range(2, 501)))# True print([(n, factorial(n)//T(n)) for n in range(1, 501) if factorial(n) % T(n) == 0 and factorial(n)//T(n) % 2])# [(1, 1), (3, 1)]
nombre maximal de régions obtenues en coupant l'espace de dimension \(d\) par \(n\) hyperplans en position générale. En dimension 2, \(1+n+\binom n2=\tfrac{n(n+1)}2+1\) : 2, 4, 7, 11, 16, 22, 29, 37… (la suite « du traiteur paresseux », qui compte aussi les morceaux d'un disque coupé par \(n\) droites). En dimension 3, avec \(n\) plans : 1, 2, 4, 8, 15, 26, 42, 64… (les nombres du « gâteau »). La preuve tient dans la récurrence \(R_d(n)=R_d(n-1)+R_{d-1}(n-1)\) : le \(n\)-ième hyperplan est lui-même découpé en \(R_{d-1}(n-1)\) morceaux, et chacun coupe une région en deux. Tant que \(n\leq d\), on trouve \(2^{n}\), puis la croissance ralentit. Le piège voisin : en reliant \(n\) points d'un cercle par toutes leurs cordes, on obtient 1, 2, 4, 8, 16, puis 31 et non 32 (problème de Moser, \(1+\binom n2+\binom n4\)). Un motif n'est pas une preuve.
# formule contre récurrence pour d ≤ 6 et n ≤ 30 ; les suites en dimensions 2 et 3 ; puis le piège de Moser from math import comb from functools import lru_cache R = lambda d, n: sum(comb(n, k) for k in range(d+1)) @lru_cache(None) def rec(d, n): # la récurrence de la preuve if n == 0 or d == 0: return 1 return rec(d, n-1) + rec(d-1, n-1) print(all(R(d, n) == rec(d, n) for d in range(7) for n in range(31)))# True print([R(2, n) for n in range(1, 9)])# [2, 4, 7, 11, 16, 22, 29, 37] print([R(3, n) for n in range(0, 9)])# [1, 2, 4, 8, 15, 26, 42, 64, 93] print([1 + comb(n, 2) + comb(n, 4) for n in range(1, 9)]) # Moser : 31, pas 32# [1, 2, 4, 8, 16, 31, 57, 99] print([2**(n-1) for n in range(1, 9)])# [1, 2, 4, 8, 16, 32, 64, 128]
pour tout \(a\) premier avec \(q\), quand \(N\to\infty\). Dirichlet (1837) a prouvé que chaque progression \(a,\ a+q,\ a+2q,\dots\) contient une infinité de premiers, en inventant pour cela ses fonctions \(L\) ; de la Vallée Poussin (1896) a montré qu'elles se partagent les premiers à parts égales, \(1/\varphi(q)\) chacune. Exemple avec \(q=9\) : les six suites qui partent de 1, 2, 4, 5, 7, 8 (et avancent de 9 en 9) reçoivent chacune un sixième des premiers. Attention au sens : c'est \(N/\ln N\), pas \(N\ln N\). Finesse célèbre, le biais de Chebyshev (1853) : les classes qui ne sont pas des carrés modulo \(q\) mènent presque toujours la course. Modulo 4, la classe de 3 devance celle de 1 jusqu'à 26 861, et Rubinstein et Sarnak (1994) ont estimé, sous des hypothèses de type Riemann, qu'elle mène environ 99,6 % du « temps » logarithmique.
# les six classes modulo 9 jusqu'à 10⁷ ; leur part ; le biais de Chebyshev ; puis la course modulo 4 import math N = 10**7 crible = bytearray([1])*(N+1); crible[0] = crible[1] = 0 for q in range(2, int(N**0.5)+1): if crible[q]: crible[q*q::q] = bytearray(len(crible[q*q::q])) compte = {r: 0 for r in (1, 2, 4, 5, 7, 8)} for n in range(2, N+1): if crible[n] and n % 9 in compte: compte[n % 9] += 1 print(compte) # {1: 110772, 2: 110836, 4: 110743, 5: 110760, 7: 110679, 8: 110788} print(sum(compte.values()) // 6, round(N / math.log(N) / 6))# 110763 103403 print(compte[2] + compte[5] + compte[8], compte[1] + compte[4] + compte[7]) # biais : non-carrés devant# 332384 332194 c1 = c3 = 0 for n in range(3, N+1): # la course modulo 4 if crible[n]: if n % 4 == 1: c1 += 1 else: c3 += 1 if c1 > c3: break print(n) # 26861
où \(\Phi_d\) a pour racines les racines \(d\)-ièmes primitives de l'unité (L·52), et pour degré \(\varphi(d)\). Les racines \(n\)-ièmes de \(-1\), \((-1)^{1/n}=e^{i\pi/n}\), sont des racines de \(\Phi_{2n}\) : \(\Phi_4=x^{2}+1\), \(\Phi_6=x^{2}-x+1\), \(\Phi_8=x^{4}+1\), \(\Phi_{10}=x^{4}-x^{3}+x^{2}-x+1\), \(\Phi_{12}=x^{4}-x^{2}+1\), \(\Phi_{14}=x^{6}-x^{5}+\dots+1\). Le coefficient juste sous le terme de tête vaut \(-\mu(n)\), puisque la somme des racines est \(\mu(n)\) (L·52). Tous les coefficients valent 0, 1 ou −1… jusqu'à \(\Phi_{105}\), qui contient un −2 (Migotti, 1883) ; ensuite, ils deviennent aussi grands qu'on veut. Encore un motif qui trompe (voir Moser, L·65).
# tous les Φₙ jusqu'à 210 par division exacte ; tes exemples ; le degré φ(n) et le coefficient −μ(n) ; puis le −2 de Φ₁₀₅ from math import gcd def divise(a, b): # division exacte de polynômes entiers (degrés décroissants) a, q = a[:], [] while len(a) >= len(b): c = a[0] // b[0]; q.append(c) for i in range(len(b)): a[i] -= c * b[i] a.pop(0) return q def produit(a, b): r = [0]*(len(a) + len(b) - 1) for i, x in enumerate(a): for j, y in enumerate(b): r[i+j] += x*y return r Phi = {} for n in range(1, 211): P = [1] + [0]*(n-1) + [-1] # xⁿ − 1 for d in range(1, n): if n % d == 0: P = divise(P, Phi[d]) Phi[n] = P ok = True for n in range(1, 211): P = [1] for d in range(1, n+1): if n % d == 0: P = produit(P, Phi[d]) ok &= P == [1] + [0]*(n-1) + [-1] print(ok) # True print({n: Phi[n] for n in (6, 10, 12, 14)})# {6: [1, -1, 1], 10: [1, -1, 1, -1, 1], 12: [1, 0, -1, 0, 1], 14: [1, -1, 1, -1, 1, -1, 1]} phi = lambda n: sum(1 for k in range(1, n+1) if gcd(k, n) == 1) def mu(n): r, m, q = 1, n, 2 while q*q <= m: if m % q == 0: m //= q if m % q == 0: return 0 r = -r q += 1 return -r if m > 1 else r print(all(len(Phi[n]) - 1 == phi(n) for n in range(1, 211)),# True True all(Phi[n][1] == -mu(n) for n in range(2, 211))) print(max(max(abs(c) for c in Phi[n]) for n in range(1, 105)), min(Phi[105]))# 1 -2
où \(C(n)\) est le nombre minimal de comparaisons qu'il faut, dans le pire des cas, pour trier \(n\) objets en les comparant deux à deux. Argument : il y a \(n!\) ordres possibles, et chaque comparaison, qui répond oui ou non, divise au mieux les possibilités par deux ; il en faut donc au moins \(\log_2 n!\). Stirling (X·5) donne \(\log_2 n!\approx n\log_2 n-n\log_2 e\), d'où le « \(n\ln n\) » de l'informatique, à un facteur constant près : c'est \(\log_2\) qui compte. L'algorithme de Ford et Johnson (1959) atteint la borne jusqu'à 11 objets ; pour 12 objets, la borne vaut 29 comparaisons, mais il en faut réellement 30 (Wells, 1965, par recherche sur ordinateur).
# la borne et Ford–Johnson (29 contre 30 pour n = 12) ; le pire cas du tri fusion sur toutes les permutations ; puis Stirling from math import factorial, log2, ceil, e from itertools import permutations borne = lambda n: ceil(log2(factorial(n))) FJ = lambda n: sum(ceil(log2(3*k/4)) for k in range(1, n+1)) # Ford–Johnson, pire cas print([borne(n) for n in range(1, 13)])# [0, 1, 3, 5, 7, 10, 13, 16, 19, 22, 26, 29] print([FJ(n) for n in range(1, 13)])# [0, 1, 3, 5, 7, 10, 13, 16, 19, 22, 26, 30] def tri_fusion(L, c): # tri fusion qui compte ses comparaisons if len(L) <= 1: return L a, b = tri_fusion(L[:len(L)//2], c), tri_fusion(L[len(L)//2:], c) r = [] while a and b: c[0] += 1 r.append(a.pop(0) if a[0] < b[0] else b.pop(0)) return r + a + b def pire_fusion(n): pire = 0 for P in permutations(range(n)): c = [0]; tri_fusion(list(P), c); pire = max(pire, c[0]) return pire print([pire_fusion(n) for n in range(1, 9)], [borne(n) for n in range(1, 9)])# [0, 1, 3, 5, 8, 11, 14, 17] [0, 1, 3, 5, 7, 10, 13, 16] n = 10**6 print(round(log2(factorial(20)), 3), round(20*log2(20) - 20*log2(e) + 0.5*log2(2*3.141592653589793*20), 3))# 61.077 61.071
pour tous réels \(x_i\geq0\), avec égalité seulement quand ils sont tous égaux. Pour deux nombres, c'est \((\sqrt x-\sqrt y)^{2}\geq0\) ; Cauchy (1821) l'a étendue à \(n\) nombres par une récurrence qui double puis redescend, et on peut aussi la tirer de la concavité du logarithme. Elle s'inscrit dans une chaîne : moyenne harmonique \(\leq\) géométrique \(\leq\) arithmétique \(\leq\) quadratique. Sens géométrique : à périmètre égal, le carré est le rectangle d'aire maximale (et le cube, la boîte de volume maximal), cousin en polygones de l'inégalité isopérimétrique \(4\pi A\leq P^{2}\), où c'est le cercle qui gagne.
# 20 000 listes au hasard pour la chaîne H ≤ G ≤ A ≤ Q ; le cas d'égalité ; puis les rectangles de périmètre 4 import random, math random.seed(1) ok_chaine, ecart_min = True, 1.0 for _ in range(20000): n = random.randint(2, 20) x = [random.uniform(0.001, 100) for _ in range(n)] H = n / sum(1/t for t in x) G = math.exp(sum(math.log(t) for t in x) / n) A = sum(x) / n Q = math.sqrt(sum(t*t for t in x) / n) ok_chaine &= H <= G*(1 + 1e-12) and G <= A*(1 + 1e-12) and A <= Q*(1 + 1e-12) print(ok_chaine) # True x = [7.5]*9 # égalité : tous égaux print(sum(x)/9, math.prod(x)**(1/9))# 7.5 7.499999999999999 aires = [(a, round(a*(2 - a), 4)) for a in (0.2, 0.6, 1.0, 1.4, 1.8)] # périmètre 4 print(aires) # [(0.2, 0.36), (0.6, 0.84), (1.0, 1.0), (1.4, 0.84), (1.8, 0.36)]
pour \(x>1\), avec la fonction \(R\) de Riemann, \(R(x)=\sum_{n\geq1}\dfrac{\mu(n)}{n}\operatorname{li}\bigl(x^{1/n}\bigr)\). La somme porte sur les zéros non triviaux \(\rho\) de \(\zeta\), groupés par paires conjuguées ; les deux derniers termes sont la contribution des zéros triviaux. \(\pi_0\) compte les premiers jusqu'à \(x\), en prenant la moyenne aux sauts. Riemann (1859) l'a écrite, von Mangoldt (1895) l'a démontrée. Chaque zéro ajoute une vague qui redessine l'escalier des premiers. À lui seul, \(R\) approche déjà \(\pi\) bien mieux que \(\operatorname{li}\) (L·32) : pour \(x=10^{6}\), \(R\) donne 78 527,4 et \(\operatorname{li}\) 78 627,5, pour 78 498 premiers. L'hypothèse de Riemann (C·2) revient à dire que toutes ces vagues restent de taille \(\sqrt x\).
# π, R et li à un million ; puis la formule en x = 100 avec 0, 10, 50, 100 et 200 zéros : elle converge vers 25 en oscillant from mpmath import mp, mpf, riemannr, li, ei, log, re, quad, inf, zetazero, nstr import math mp.dps = 20 N = 10**6 crible = bytearray([1])*(N+1); crible[0] = crible[1] = 0 for q in range(2, int(N**0.5)+1): if crible[q]: crible[q*q::q] = bytearray(len(crible[q*q::q])) print(sum(crible), nstr(riemannr(N), 9), nstr(li(N), 9))# 78498 78527.3994 78627.5492 def mu(n): r, m, q = 1, n, 2 while q*q <= m: if m % q == 0: m //= q if m % q == 0: return 0 r = -r q += 1 return -r if m > 1 else r # même formule, regroupée sous sa forme finie : π₀(x) = Σ_{n ≤ log₂ x} μ(n)/n · J(x^(1/n)), avec # J(y) = li(y) − Σ_ρ li(y^ρ) − ln 2 + ∫_y^∞ dt / (t(t²−1) ln t), et li(y^ρ) = Ei(ρ ln y) Z = [zetazero(k) for k in range(1, 201)] def J(y, K): return (li(y) - sum(2*re(ei(r*log(y))) for r in Z[:K]) - log(2) + quad(lambda t: 1/(t*(t*t - 1)*log(t)), [y, inf])) x = mpf(100) # π(100) = 25, et 100 n'est pas premier pi0 = lambda K: sum(mu(n)/mpf(n) * J(x**(mpf(1)/n), K) for n in range(1, int(math.log2(100))+1) if mu(n)) print([nstr(pi0(K), 6) for K in (0, 10, 50, 100, 200)])# ['25.6885', '25.2813', '25.0775', '24.9147', '24.9256']
pour \(x>1\) qui n'est pas une puissance de premier (von Mangoldt, 1895), où \(\psi(x)=\sum_{n\leq x}\Lambda(n)=\ln\operatorname{ppcm}(1,\dots,\lfloor x\rfloor)\) (VIII·3) et où \(\rho\) parcourt les zéros non triviaux de \(\zeta\), par paires conjuguées. Le graphe de \(\psi\) est un escalier qui épouse la droite \(y=x\), et les zéros y ajoutent des vagues. Il saute de \(\ln p\) à chaque puissance de premier, donc aussi en 4, 8, 9, 16, 25, 27… ; c'est \(\theta(x)=\sum_{p\leq x}\ln p\) qui ne saute qu'aux premiers, tandis que \(\pi(x)\) saute de 1 et suit la courbe \(x/\ln x\). D'où vient la formule : une série de Dirichlet est une transformée intégrale de son escalier, \(\displaystyle-\frac{\zeta'(s)}{\zeta(s)}=\sum\frac{\Lambda(n)}{n^{s}}=s\int_1^{\infty}\psi(x)\,x^{-s-1}\,dx\) (une intégrale de Stieltjes, qui devient une transformée de Laplace avec \(x=e^{u}\)) ; on l'inverse par la formule de Perron, et les pôles de \(-\zeta'/\zeta\), c'est-à-dire les zéros de \(\zeta\), donnent les termes de la somme.
# les sauts de ψ ; la formule en x = 100,5 avec 0 à 300 zéros ; puis série de Dirichlet = intégrale de l'escalier from math import lcm, log from functools import reduce from mpmath import mp, mpf, zetazero, zeta, pi, re, nstr mp.dps = 20 def Lam(n): # von Mangoldt (L·30) for q in range(2, n+1): if n % q == 0: while n % q == 0: n //= q return log(q) if n == 1 else 0.0 return 0.0 print([(n, round(Lam(n), 3)) for n in range(2, 30) if Lam(n)]) # les sauts : puissances de premiers# [(2, 0.693), (3, 1.099), (4, 0.693), (5, 1.609), (7, 1.946), (8, 0.693), (9, 1.099), (11, 2.398), (13, 2.565), (16, 0.693), (17, 2.833), (19, 2.944), (23, 3.135), (25, 1.609), (27, 1.099), (29, 3.367)] x = mpf('100.5') psi = log(reduce(lcm, range(1, 101))) Z = [zetazero(k) for k in range(1, 301)] f = lambda K: x - sum(2*re(x**r / r) for r in Z[:K]) - log(2*pi) - log(1 - x**-2)/2 print(round(psi, 6), [nstr(f(K), 7) for K in (0, 10, 100, 300)])# 94.045311 ['98.66217', '95.77622', '95.00926', '93.78253'] N, s = 20000, 2 # le pont Dirichlet ↔ intégrale (sommation d'Abel) L = [0.0] + [Lam(n) for n in range(1, N+1)] ps = [0.0]*(N+1) for n in range(1, N+1): ps[n] = ps[n-1] + L[n] integrale = sum(ps[k] * (k**-s - (k+1)**-s) / s for k in range(1, N)) # ∫₁ᴺ ψ(x) x^(−s−1) dx, exacte print(round(s*integrale + ps[N]*N**-s, 12), round(sum(L[n]*n**-s for n in range(1, N+1)), 12))# 0.569910938076 0.569910938076 print(nstr(-zeta(s, derivative=1)/zeta(s), 10))# 0.5699609931
en longueurs de carte, pour \(n\) cartes empilées au bord d'une table, une par étage, chacune décalée au maximum sans que la pile tombe. En partant du haut, la \(k\)-ième carte dépasse de \(\tfrac1{2k}\) celle du dessous : le centre de gravité des \(k\) cartes du haut tombe alors pile sur le bord de la carte suivante, et celui de toute la pile sur le bord de la table. Avec 100 cartes, \(\tfrac12H_{100}\approx2{,}59\) longueurs dans le vide. Comme \(H_n\) tend vers l'infini (L·37 et XI·2), on peut déborder aussi loin qu'on veut, mais chaque longueur supplémentaire coûte environ \(e^{2}\approx7{,}4\) fois plus de cartes, puisque \(\tfrac12H_n\approx\tfrac12(\ln n+\gamma)\). Au passage, \(H_{10n}-H_n\) tend vers \(\ln10\approx2{,}3026\), l'écart entre deux puissances de 10. En s'autorisant plusieurs cartes par étage, on fait bien mieux : un débordement qui grandit comme \(n^{1/3}\) (Paterson et Zwick, 2006).
| débordement | cartes nécessaires |
|---|---|
| 1 longueur | 4 |
| 2 longueurs | 31 |
| 3 longueurs | 227 |
| 4 longueurs | 1 674 |
| 5 longueurs | 12 367 |
| 6 longueurs | 91 380 |
# l'équilibre exact de chaque étage (en fractions) ; les 100 cartes ; le tableau et son rapport e² ; puis tes valeurs de Hₙ et l'écart ln 10 from fractions import Fraction as Fr from mpmath import mp, mpf, harmonic, log, e, nstr def pile(n): # bords droits, du haut vers le bas ; table en x = 0 R = [sum(Fr(1, 2*k) for k in range(j, n+1)) for j in range(1, n+1)] + [Fr(0)] centres = [r - Fr(1, 2) for r in R[:n]] equilibre = all(sum(centres[:j]) / j == R[j] for j in range(1, n+1)) # centre de gravité sur le bord return R[0], equilibre print(all(pile(n)[1] for n in range(1, 60)))# True print(pile(100)[0] == sum(Fr(1, 2*k) for k in range(1, 101)), float(pile(100)[0]))# True 2.5936887588198103 H, n, besoin, cible = 0.0, 0, [], 1 while len(besoin) < 6: n += 1; H += 1/n if H/2 >= cible: besoin.append(n); cible += 1 print(besoin, [round(besoin[i+1]/besoin[i], 2) for i in range(5)], round(float(e**2), 2))# [4, 31, 227, 1674, 12367, 91380] [7.75, 7.32, 7.37, 7.39, 7.39] 7.39 mp.dps = 20 print([nstr(harmonic(10**k), 11) for k in range(1, 10)])# ['2.928968254', '5.1873775176', '7.4854708606', '9.787606036', '12.09014613', '14.392726723', '16.695311366', '18.997896414', '21.300481502'] print(nstr(harmonic(10**9) - harmonic(10**8), 10), nstr(log(10), 10))# 2.302585088 2.302585093
pour \(n\geq2\) ; autrement dit, \(n\equiv1 \pmod{\varphi(n)}\) suffirait-il à rendre \(n\) premier ? La réciproque est vraie, car \(\varphi(p)=p-1\). On ne connaît aucun contre-exemple, et on sait qu'il en faudrait un énorme, avec beaucoup de facteurs premiers.
# une recherche n'est pas une preuve : aucun contre-exemple sous 10⁶ N = 10**6 phi = list(range(N+1)) for p in range(2, N+1): if phi[p] == p: for m in range(p, N+1, p): phi[m] -= phi[m] // p contre = [n for n in range(2, N+1) if (n-1) % phi[n] == 0 and phi[n] != n-1] print(contre) # []
tous les zéros non triviaux de ζ seraient sur la droite critique. C'est le centre de la symétrie de L·20 : l'équation fonctionnelle place les zéros en miroir autour de \(\operatorname{Re}s=\tfrac12\), l'hypothèse dit qu'ils sont tous sur le miroir lui-même. Vérifiée par le calcul pour des milliers de milliards de zéros, jamais démontrée ; l'un des problèmes du millénaire.
# méthode de Turing : on compte les zéros dans toute la bande critique, # puis ceux qui sont sur la droite ; s'il y en a autant, aucun n'est ailleurs. # Une vérification finie n'est pas une preuve. from mpmath import mp, nzeros, zetazero mp.dps = 25 T = 100 # zéros dans la bande 0 < Re s < 1, 0 < Im s ≤ T print(nzeros(T)) # 29 # zéros trouvés sur la droite Re s = ½ k = 0 while zetazero(k+1).imag <= T: k += 1 print(k) # 29 print(zetazero(1)) # (0.5 + 14.13472514173469379045725j)
on sait l'empilement optimal en dimensions 1, 2, 3, 8 et 24 seulement ; la dimension 4 est la première inconnue. Le meilleur candidat est le réseau \(D_4\) (points entiers de somme paire) : parmi les empilements en réseau, il est prouvé optimal depuis Korkine et Zolotarev (1877), et chaque boule y touche 24 voisines, le maximum possible en dimension 4 (Musin, 2008). Ces 24 voisines sont les sommets du 24-cellules, le solide de Platon propre à la dimension 4. Reste à exclure tout empilement irrégulier plus dense : c'est ouvert.
| dim. | densité optimale | démontrée par |
|---|---|---|
| 1 | \(1\) | évident |
| 2 | \(\pi/(2\sqrt3)\approx0{,}9069\) | Thue, puis Fejes Tóth (1940) |
| 3 | \(\pi/(3\sqrt2)\approx0{,}7405\) | Kepler (1611) ; Hales (2005), preuve formelle en 2017 |
| 4 | \(\pi^{2}/16\approx0{,}6169\) ? | ouvert |
| 8 | \(\pi^{4}/384\approx0{,}2537\) | Viazovska (2016), L·21 |
| 24 | \(\pi^{12}/12!\approx0{,}0019\) | Cohn et al. (2016), L·23 |
A. Korkine et G. Zolotareff, « Sur les formes quadratiques positives », Mathematische Annalen 11 (1877).
O. R. Musin, « The kissing number in four dimensions », Annals of Mathematics 168 (2008).
# voisines et densité de D4 ; l'optimalité parmi tous les empilements reste à prouver from itertools import product from mpmath import mp, mpf, pi, gamma, sqrt mp.dps = 25 V = lambda n, r=1: pi**(mpf(n)/2) / gamma(mpf(n)/2 + 1) * r**n # L·10 # D4 : vecteurs entiers de somme paire, covolume 2, distance minimale √2 voisins = [v for v in product(range(-1, 2), repeat=4) if sum(x*x for x in v) == 2 and sum(v) % 2 == 0] print(len(voisins)) # 24 print(V(4, 1/sqrt(2)) / 2) # 0.6168502750680849136771557 print(pi**2/16) # 0.6168502750680849136771557
pour tout premier \(p\), l'ensemble des \(n\) où \(p\) divise le numérateur de \(H_n\) serait fini. Pour \(p=2\) il est vide (L·37) ; pour \(p=3\), c'est \(\{2,7,22\}\) ; pour \(p=5\), \(\{4,20,24\}\) ; pour \(p=7\), il s'allonge déjà. Ces ensembles commandent la suite A110566, \(\operatorname{PPCM}(1,\dots,n)\) divisé par le dénominateur de \(H_n\) : 1, 1, 1, 1, 1, 3, 3, 3… Un premier \(p\) la divise exactement quand \(\lfloor n/p^{k}\rfloor\) tombe dans l'ensemble de \(p\), où \(p^{k}\) est la plus grande puissance de \(p\) qui ne dépasse pas \(n\). Pour \(n=6\) : \(\lfloor6/3\rfloor=2\), qui est dans l'ensemble de 3 ; de fait, \(\tfrac13+\tfrac16=\tfrac12\), et le 3 disparaît du dénominateur.
A. Eswarathasan et E. Levine, « p-integral harmonic sums », 1991.
# les ensembles pour p = 2, 3, 5, 7 jusqu'à 600 ; A110566 ; puis la règle qui la relie aux ensembles from fractions import Fraction from math import lcm N = 600 H, L = [Fraction(0)], [1] for n in range(1, N+1): H.append(H[-1] + Fraction(1, n)); L.append(lcm(L[-1], n)) J = lambda p: [n for n in range(1, N+1) if H[n].numerator % p == 0] print(J(2)) # [] print(J(3)) # [2, 7, 22] print(J(5)) # [4, 20, 24] print(J(7)) # [6, 42, 48, 295, 299, 337, 341] print([L[n] // H[n].denominator for n in range(1, 25)]) # A110566# [1, 1, 1, 1, 1, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 15, 45, 45, 45, 15] premiers = [p for p in range(2, N+1) if all(p % q for q in range(2, int(p**0.5)+1))] Js = {p: set(J(p)) for p in premiers} def prevu(n): # premiers p avec ⌊n/p^k⌋ dans J(p) out = set() for p in premiers: if p > n: break k = 1 while p**(k+1) <= n: k += 1 if n // p**k in Js[p]: out.add(p) return out print(all({p for p in premiers if (L[n] // H[n].denominator) % p == 0} == prevu(n)# True for n in range(1, N+1)))
quand \(x\to\infty\), pour \(\lambda\) fixé. L'idée vient du modèle aléatoire de Cramér (1936) : faire comme si chaque entier \(n\) était premier avec probabilité \(1/\ln n\), comme des déchirures dans un tissu (L·56). Gallagher (1976) a démontré cette loi de Poisson à condition que la conjecture des \(k\)-uplets premiers de Hardy et Littlewood soit vraie, ce qui reste ouvert. Deux précautions : la densité des premiers près de \(x\) est \(1/\ln x\), avec le logarithme naturel (avec \(\log_{10}\), on en prédirait 2,3 fois trop) ; et il faut des intervalles de longueur proportionnelle à \(\ln x\), car un intervalle de longueur 1 ne contient qu'un seul entier. Le modèle de Cramér a ses limites : Maier (1985) a montré qu'il se trompe sur des intervalles un peu plus longs. Et à portée de calcul, les premiers sont plus réguliers que le hasard, puisqu'ils évitent les multiples de 2, 3, 5…
H. Cramér, « On the order of magnitude of the difference between consecutive prime numbers », Acta Arithmetica 2 (1936).
P. X. Gallagher, « On the distribution of primes in short intervals », Mathematika 23 (1976).
# premiers près d'un million contre les deux densités ; puis les répartitions observées près de 10⁷ contre Poisson import math from collections import Counter N = 2 * 10**7 crible = bytearray([1])*(N+1); crible[0] = crible[1] = 0 for q in range(2, int(N**0.5)+1): if crible[q]: crible[q*q::q] = bytearray(len(crible[q*q::q])) x, w = 10**6, 10**4 # la bonne densité : 1/ln x print(sum(crible[x:x+w]), round(w/math.log(x)), round(w/math.log10(x)))# 753 724 1667 x0 = 10**7 def compare(h): # (observé, Poisson) pour k = 0 … 5 C = Counter(sum(crible[a:a+h]) for a in range(x0, x0 + 5*10**6, h)) tot, lam = sum(C.values()), h / math.log(x0) return [(round(C[k]/tot, 3), round(math.exp(-lam)*lam**k/math.factorial(k), 3)) for k in range(6)] print(compare(16)) # λ ≈ 1# [(0.322, 0.371), (0.428, 0.368), (0.202, 0.183), (0.043, 0.06), (0.004, 0.015), (0.0, 0.003)] print(compare(48)) # λ ≈ 3# [(0.027, 0.051), (0.125, 0.152), (0.243, 0.226), (0.275, 0.224), (0.195, 0.167), (0.096, 0.099)]