André Eberhard, quelle était la précision des nombres flottants que vous utilisiez (avant 1967) ?¶
Page 280 du livre "Procédures Algol en Analyse Numérique", tome 1: ajustement polynomial aux moindres carrés, André Eberhard implante la méthode naïve où on ajuste les coefficients $a$ de $$p(x) = \sum_{i=0}^m a_i x^i.$$
Il est bien connu que cette méthode est instable, la matrice du système linéaire étant un avatar de la matrice de Hilbert.
C'est un ajustement avec poids: $$\min_{\{a\}} \sum_{i=1}^n g_i(p(x_i) -y_i)^2.$$
On considère les matrices $B_{ij}= x_i^j$ et $ G_{ii}= g_i$ (indicées à partir de $0$).
On résout $$ (B^t G B)\ a = (B^t G)\ y.$$
L'inversion de la matrice est effectuée en utilisant un avatar de la routine GRESOLSYSLIN (Tome 1 des procédures, auteur N. Gastinel).
Cette routine ne fait pas de choix de pivot max., mais juste un test:
IF A[I,I]=0 ALLER A PERMUTER LIGNE;et PERMUTER LIGNE se contente de permuter la ligne I avec la première ligne J rencontrée où A[J,I] est non nul.
Dans le code Algol, quelle était donc la précision des réels?¶
car il y a un seul type REEL (real).
Le problème résolu dans le livre :¶
$x$: 21 points équirépartis sur $[0,1]$
$y = \exp(x)$.
On ajuste un polynôme de degré $m = 4$.
les poids $g$ sont égaux à 1 partout.
L'auteur trouve les coefficients :
$0.10003376\times 10, \ 0.99835205, \ 0.50927734, \ 0.14135742, \ 0.69702148\times 10^{-1}.$
SageMath is a free open-source mathematics software system licensed under the GPL.
La premiere procédure fabrique la matrice $B_{ij}= x_i^j$ :
def matrixB(order,x):
#
n=len(x)
B=matrix(RealNumber,n,order+1)
for i in [0..n-1]:
for j in [0..order]:
B[i,j]=x[i]^j
return B
Deuxième procédure : inversion par la méthode de Gauss (GRESOLSYSLINE):
def invgauss(A,B,n):
pivots=[]
for i in range(0,n):
if A[i,i]==RealNumber(0):
raise ValueError
else:
d=A[i,i]
for j in range(i+1,n):
s=A[j,i]/A[i,i]
for k in range(i,n):
A[j,k]-=s*A[i,k]
for k in range(0,n):
B[j,k]-=s*B[i,k]
for k in range(0,n):
B[n-1,k]=B[n-1,k]/A[n-1,n- 1]
for i in range(n-2,-1,-1):
for k in range(0,n):
s=RealNumber(0)
for j in range(i+1,n):
s+=A[i,j]*B[j,k]
B[i,k]=(B[i,k]-s)/A[i,i]
La troisième résout le problème, comme dans le livre:
def slv(order,x,y,g):
G=diagonal_matrix(RealNumber,g)
B= matrixB(order,x)
M=transpose(B)*G*B
M1=identity_matrix(RealNumber,order+1)
invgauss(M,M1,order+1)
V=transpose(B)*G*y
return M1*V
RealNumber ?¶
GNU MPFR est une bibliothèque de flottants à précision quelconque, qui permet aussi de choisir différents modes d'arrondi. Dans Sage, RealField(n) désigne les flottants de MPFR, avec n bits de précision.
Exemple : RealNumber= RealField(53) donne les double precision habituels.
Quel est le conditionnement de la matrice qu'on inverse ?
BigPrec=350
RealNumber=RealField(BigPrec)
order=4
n=21
x=vector(RealNumber,[(i-1)/(n-1) for i in [1..n]])
cond= lambda m : m.norm(Infinity)*(m^(-1)).norm(Infinity)
B=matrixB(order,x)
cn=cond(transpose(B)*B)
print("conditionnement:",cn.n(15))
On calcule une solution de référence presque exacte :
(avec 350 bits, on aura une erreur de l'ordre de $10^{-100}$)
La solution de référence, tronquée pour l'impression :
print([s.n(53) for s in cref])
Rappel : dans le livre, on a :
$0.10003376\times 10, \ 0.99835205, \ 0.50927734, \ 0.14135742, \ 0.69702148\times 10^{-1}.$
On va calculer les coefficients en utilisant différentes précisions et différents mode d'arrondi (grace à GNU MPFR).¶
On compare les résultats à la solution de référence (norme $\infty$).
fig, ax = plt.subplots()
for t in rndchoices:
ax.plot(precs, DicRef[t])
#ax.plot(precs,[cn*2^(-k) for k in precs])
ax.set(xlabel='Digits', ylabel='Distance',title='')
ax.set_yscale('log')
ax.grid()
A quelle distance de la solution de référence se situe la solution du livre ?
book=[0.10003376*10,0.99835205,0.50927734,0.14135742,0.69702148*10^(-1)]
RR.scientific_notation(True)
v=max([abs(x) for x in vector(Rbig,book)-cref])
print("Distance ",v.n()," de la solution de référence.")
Regardons un peu à quoi cette valeur correspond :
D = dict(zip(precs,DicRef["RNDN"]))
p,w= [(s,D[s]) for s in D.keys() if D[s]>=v][-1]
print("précision (bits) ",p,", distance: ",w.n())
print("avec ",p+1,"bits, distance:",D[p+1].n())
On trouve 27 bits de précision.¶
En précision double actuelle, on aurait (53 bits de précision):
D[53].n(6)
27 bits de précision, c'est assez étrange !¶
Pour le flottant IEEE 754 (actuel, normalisé):
Format simple précision (32 bits) : 23 pour la mantisse.
Format double précision (64 bits) : 52 bits pour la mantisse.
Machines avec des flottants à 27 bits de mantisse ?¶
PDP 6,...,10.
IBM Série 7000.
Finalement :¶
- Jean-François Maître a découvert un article qui donne le pot aux roses : IBM 7044.
IBM série 7000:
Mots de 36 bits (mémoire : 32k mots ?).
- Simple precision : mantisse de 27 bits.
Mais pourquoi ne pas avoir utilisé la double précision ?¶
Double-precision introduite sur l'IBM 7094 (1962): mantisse de 54-bits.
In 1963, IBM introduced lower cost machines with a similar architecture, but fewer instructions and simplified I/O, called the IBM 7040 and 7044.
André Eberhardt ne pouvait pas calculer en double précision.
Point d'argent, point de suisse — (Jean Racine, Les Plaideurs, 1668, acte 1, scène 1).
Instabilité :¶
peut-on vraiment ajuster un polynôme de degré 8 avec 27 bits de précision ?
fig1, a= plt.subplots(rows, cols,figsize=figsize, constrained_layout=True)
for j in range(0,2):
for i in range(0,3):
k=3*j+i+4
if k>9: break
w=map(lambda s: ep(s,D[k]),x)
a[j,i].set_title("deg. "+str(k))
a[j,i].plot(x,[exp(s) for s in x],linestyle="dashed")
a[j,i].plot(x, [ep(s,D[k]) for s in x])
