Outils pour utilisateurs

Outils du site


math_decimal_bernoulli

Calcul des nombres de Bernoulli

Objectif

Pour la théorie des nombres de Bernoulli, voir entre autres ici: http://fr.wikipedia.org/wiki/Nombre_de_Bernoulli

Je m'y suis intéressé parce que ces nombres interviennent dans certains calculs, en particulier dans le calcul de la tangente par série entière.

Code proposé

Je suis parti de l'équation récurrente suivante:

  • sum{k=0}{n}{{{C_{n+1}}^{k}}*{B_k}} = 0
  • avec B_0 = 1
  • et B_i = 0 si i est un nombre impair supérieur à 1

Dans cette formule:

  • {C_{n+1}}^k est le nombre de combinaison de n+1 objets pris k à k (voir analyse combinatoire sur ce site)
  • {B_k} est le nombre de Bernoulli pour la valeur k

Mais en quoi cette équation est-elle “récurrente”?

Pour le savoir, développons pour n=3:

  • {{{C_{4}}^0}*{B_0}} ~+~ {{{C_{4}}^1}*{B_1}} ~+~ {{{C_{4}}^2}*{B_2}} ~+~ {{{C_{4}}^3}*{B_3}} ~=~ 0

Il faudrait donc que l'on puisse calculer B_3.

Mais là, on voit bien le problème: on connait {B_0}=1 par définition, mais pour calculer B_3, il faut calculer avant: B_1 et B_2.

Comment? En appliquant la formule pour n=1 et n=2 et, à chaque fois, en calculant le dernier B.

Voilà le code qui fait cela (Python version 2.6):

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from decimal import *
 
##############################################################################
def combin(n, k):
    """Nombre de combinaisons de n objets pris k a k"""
    if k > n//2:
        k = n-k
    x = 1
    y = 1
    i = n-k+1
    while i <= n:
        x = (x*i)//y
        y += 1
        i += 1
    return x
 
##############################################################################
def bernoulli(n):
    """calcul du nombre de Bernoulli Bn """ 
    prec = getcontext().prec
    getcontext().prec = 100   
    B = [Decimal("1")]
    for m in xrange(1,n+1):
        if m>1 and m%2!=0:
            B.append(Decimal("0"))
        else:
            c = Decimal("0")
            for k in xrange(0,m):
                c += Decimal(str(combin(m+1,k)))*B[k]
            B.append(-c/Decimal(str(combin(m+1,m))))
    getcontext().prec = prec 
    return float(B[n])

Pour que les calculs intermédiaires soient corrects, il est important qu'ils soient faits avec un nombre suffisant de chiffres significatifs. Ici, après essais, j'ai pris 100 chiffres! Au delà, on ne gagne plus rien (j'ai essayé avec 400 chiffres).

Dans ce code, j'ai choisi que le résultat soit retourné en flottants. Mais on pourrait facilement continuer les calculs en mode Decimal() si c'était utile. Cependant, ici, en flottant, le nombre de Bernoulli maxi utilisable sera B_258=1.3352784187354634e+306 puisqu'au delà, on dépasse les capacités de stockage des flottants (double précision du C).

Voilà les résultats:

B( 0 ) = 1.0
B( 1 ) = -0.5
B( 2 ) = 0.16666666666666666
B( 3 ) = 0.0
B( 4 ) = -0.033333333333333333
B( 5 ) = 0.0
B( 6 ) = 0.023809523809523808
B( 7 ) = 0.0
B( 8 ) = -0.033333333333333333
B( 9 ) = 0.0
B( 10 ) = 0.07575757575757576
B( 11 ) = 0.0
B( 12 ) = -0.2531135531135531
B( 13 ) = 0.0
B( 14 ) = 1.1666666666666667
B( 15 ) = 0.0
B( 16 ) = -7.0921568627450977
B( 17 ) = 0.0
B( 18 ) = 54.971177944862156
B( 19 ) = 0.0
B( 20 ) = -529.12424242424242
B( 21 ) = 0.0
B( 22 ) = 6192.123188405797
B( 23 ) = 0.0
B( 24 ) = -86580.253113553117
B( 25 ) = 0.0
B( 26 ) = 1425517.1666666667
B( 27 ) = 0.0
B( 28 ) = -27298231.067816094
B( 29 ) = 0.0
B( 30 ) = 601580873.9006424
B( 31 ) = 0.0
B( 32 ) = -15116315767.092157
B( 33 ) = 0.0
B( 34 ) = 429614643061.16669
B( 35 ) = 0.0
B( 36 ) = -13711655205088.332
B( 37 ) = 0.0
B( 38 ) = 488332318973593.19
B( 39 ) = 0.0
B( 40 ) = -19296579341940068.0
B( 41 ) = 0.0
B( 42 ) = 8.4169304757368256e+17
B( 43 ) = 0.0
B( 44 ) = -4.0338071854059454e+19
B( 45 ) = 0.0
B( 46 ) = 2.1150748638081993e+21
B( 47 ) = 0.0
B( 48 ) = -1.2086626522296526e+23
B( 49 ) = 0.0
B( 50 ) = 7.5008667460769642e+24
B( 51 ) = 0.0
B( 52 ) = -5.0387781014810688e+26
B( 53 ) = 0.0
B( 54 ) = 3.6528776484818122e+28
B( 55 ) = 0.0
B( 56 ) = -2.8498769302450882e+30
B( 57 ) = 0.0
B( 58 ) = 2.3865427499683627e+32
B( 59 ) = 0.0
B( 60 ) = -2.1399949257225335e+34
B( 61 ) = 0.0
B( 62 ) = 2.0500975723478097e+36
B( 63 ) = 0.0
B( 64 ) = -2.0938005911346379e+38
B( 65 ) = 0.0
B( 66 ) = 2.2752696488463515e+40
B( 67 ) = 0.0
B( 68 ) = -2.6257710286239577e+42
B( 69 ) = 0.0
B( 70 ) = 3.2125082102718032e+44
B( 71 ) = 0.0
B( 72 ) = -4.1598278166794712e+46
B( 73 ) = 0.0
B( 74 ) = 5.6920695482035283e+48
B( 75 ) = 0.0
B( 76 ) = -8.2183629419784578e+50
B( 77 ) = 0.0
B( 78 ) = 1.2502904327166994e+53
B( 79 ) = 0.0
B( 80 ) = -2.001558323324837e+55
B( 81 ) = 0.0
B( 82 ) = 3.3674982915364376e+57
B( 83 ) = 0.0
B( 84 ) = -5.947097050313545e+59
B( 85 ) = 0.0
B( 86 ) = 1.1011910323627977e+62
B( 87 ) = 0.0
B( 88 ) = -2.1355259545253502e+64
B( 89 ) = 0.0
B( 90 ) = 4.3328896986641194e+66
B( 91 ) = 0.0
B( 92 ) = -9.1885528241669332e+68
B( 93 ) = 0.0
B( 94 ) = 2.0346896776329074e+71
B( 95 ) = 0.0
B( 96 ) = -4.700383395803573e+73
B( 97 ) = 0.0
B( 98 ) = 1.1318043445484249e+76
B( 99 ) = 0.0
B( 100 ) = -2.8382249570693707e+78
B( 101 ) = 0.0
B( 102 ) = 7.4064248979678853e+80
B( 103 ) = 0.0
B( 104 ) = -2.0096454802756605e+83
B( 105 ) = 0.0
B( 106 ) = 5.6657170050805942e+85
B( 107 ) = 0.0
B( 108 ) = -1.6584511154136216e+88
B( 109 ) = 0.0
B( 110 ) = 5.0368859950492378e+90
B( 111 ) = 0.0
B( 112 ) = -1.5861468237658186e+93
B( 113 ) = 0.0
B( 114 ) = 5.1756743617545625e+95
B( 115 ) = 0.0
B( 116 ) = -1.7488921840217116e+98
B( 117 ) = 0.0
B( 118 ) = 6.1160519994952182e+100
B( 119 ) = 0.0
B( 120 ) = -2.2122776912707833e+103
B( 121 ) = 0.0
B( 122 ) = 8.2722776798770969e+105
B( 123 ) = 0.0
B( 124 ) = -3.1958925111415708e+108
B( 125 ) = 0.0
B( 126 ) = 1.2750082223387793e+111
B( 127 ) = 0.0
B( 128 ) = -5.2500923086774131e+113
B( 129 ) = 0.0
B( 130 ) = 2.2301817894241627e+116
B( 131 ) = 0.0
B( 132 ) = -9.7684521930955207e+118
B( 133 ) = 0.0
B( 134 ) = 4.409836197845295e+121
B( 135 ) = 0.0
B( 136 ) = -2.0508570886464089e+124
B( 137 ) = 0.0
B( 138 ) = 9.8214433279791277e+126
B( 139 ) = 0.0
B( 140 ) = -4.8412600798208881e+129
B( 141 ) = 0.0
B( 142 ) = 2.4553088801480982e+132
B( 143 ) = 0.0
B( 144 ) = -1.2806926804084748e+135
B( 145 ) = 0.0
B( 146 ) = 6.8676167104668579e+137
B( 147 ) = 0.0
B( 148 ) = -3.7846468581969106e+140
B( 149 ) = 0.0
B( 150 ) = 2.1426101250665291e+143
B( 151 ) = 0.0
B( 152 ) = -1.2456727137183695e+146
B( 153 ) = 0.0
B( 154 ) = 7.4345787551000156e+148
B( 155 ) = 0.0
B( 156 ) = -4.5535795304641704e+151
B( 157 ) = 0.0
B( 158 ) = 2.861211281685887e+154
B( 159 ) = 0.0
B( 160 ) = -1.843772355203387e+157
B( 161 ) = 0.0
B( 162 ) = 1.2181154536221047e+160
B( 163 ) = 0.0
B( 164 ) = -8.2482187185314122e+162
B( 165 ) = 0.0
B( 166 ) = 5.7225877937832942e+165
B( 167 ) = 0.0
B( 168 ) = -4.0668530525059105e+168
B( 169 ) = 0.0
B( 170 ) = 2.9596092064642052e+171
B( 171 ) = 0.0
B( 172 ) = -2.2049522565189457e+174
B( 173 ) = 0.0
B( 174 ) = 1.6812597072889599e+177
B( 175 ) = 0.0
B( 176 ) = -1.3116736213556958e+180
B( 177 ) = 0.0
B( 178 ) = 1.0467894009478039e+183
B( 179 ) = 0.0
B( 180 ) = -8.543289357883371e+185
B( 181 ) = 0.0
B( 182 ) = 7.1287821322486547e+188
B( 183 ) = 0.0
B( 184 ) = -6.0802931455535905e+191
B( 185 ) = 0.0
B( 186 ) = 5.2996776424849921e+194
B( 187 ) = 0.0
B( 188 ) = -4.719425916874586e+197
B( 189 ) = 0.0
B( 190 ) = 4.2928413791402983e+200
B( 191 ) = 0.0
B( 192 ) = -3.9876744968232205e+203
B( 193 ) = 0.0
B( 194 ) = 3.781978041935888e+206
B( 195 ) = 0.0
B( 196 ) = -3.6614233683681192e+209
B( 197 ) = 0.0
B( 198 ) = 3.617609027237286e+212
B( 199 ) = 0.0
B( 200 ) = -3.6470772645191356e+215
B( 201 ) = 0.0
B( 202 ) = 3.7508755436454406e+218
B( 203 ) = 0.0
B( 204 ) = -3.934586729643903e+221
B( 205 ) = 0.0
B( 206 ) = 4.2088211148190084e+224
B( 207 ) = 0.0
B( 208 ) = -4.5902296220617918e+227
B( 209 ) = 0.0
B( 210 ) = 5.1031725772629572e+230
B( 211 ) = 0.0
B( 212 ) = -5.7822762303656954e+233
B( 213 ) = 0.0
B( 214 ) = 6.6762482167835884e+236
B( 215 ) = 0.0
B( 216 ) = -7.8535307644450417e+239
B( 217 ) = 0.0
B( 218 ) = 9.4106894067058724e+242
B( 219 ) = 0.0
B( 220 ) = -1.1484933873465185e+246
B( 221 ) = 0.0
B( 222 ) = 1.4272958742848785e+249
B( 223 ) = 0.0
B( 224 ) = -1.8059559586909309e+252
B( 225 ) = 0.0
B( 226 ) = 2.3261535307660807e+255
B( 227 ) = 0.0
B( 228 ) = -3.0495751715499594e+258
B( 229 ) = 0.0
B( 230 ) = 4.0685806076433976e+261
B( 231 ) = 0.0
B( 232 ) = -5.5231031321974362e+264
B( 233 ) = 0.0
B( 234 ) = 7.6277279396434395e+267
B( 235 ) = 0.0
B( 236 ) = -1.0715571119697886e+271
B( 237 ) = 0.0
B( 238 ) = 1.5310200895969188e+274
B( 239 ) = 0.0
B( 240 ) = -2.2244891682179836e+277
B( 241 ) = 0.0
B( 242 ) = 3.286267919069014e+280
B( 243 ) = 0.0
B( 244 ) = -4.9355928955960348e+283
B( 245 ) = 0.0
B( 246 ) = 7.5349571200832511e+286
B( 247 ) = 0.0
B( 248 ) = -1.1691485154584178e+290
B( 249 ) = 0.0
B( 250 ) = 1.8435261467838938e+293
B( 251 ) = 0.0
B( 252 ) = -2.9536826172968082e+296
B( 253 ) = 0.0
B( 254 ) = 4.8079321277501568e+299
B( 255 ) = 0.0
B( 256 ) = -7.9502125045885252e+302
B( 257 ) = 0.0
B( 258 ) = 1.3352784187354634e+306

Et voilà une liste apte à être intégrée dans un code par copier-coller:

bernoulli = [1.0, -0.5, 0.166666666667, 0.0, -0.0333333333333, 0.0, 0.0238095238095, 0.0, -0.0333333333333, 0.0,
0.0757575757576, 0.0, -0.253113553114, 0.0, 1.16666666667, 0.0, -7.09215686275, 0.0, 54.9711779449, 0.0, -529.124242424,
0.0, 6192.12318841, 0.0, -86580.2531136, 0.0, 1425517.16667, 0.0, -27298231.0678, 0.0, 601580873.901, 0.0, -15116315767.1,
0.0, 429614643061.0, 0.0, -1.37116552051e+13, 0.0, 4.88332318974e+14, 0.0, -1.92965793419e+16, 0.0, 8.41693047574e+17,
0.0, -4.03380718541e+19, 0.0, 2.11507486381e+21, 0.0, -1.20866265223e+23, 0.0, 7.50086674608e+24, 0.0, -5.03877810148e+26,
0.0, 3.65287764848e+28, 0.0, -2.84987693025e+30, 0.0, 2.38654274997e+32, 0.0, -2.13999492572e+34, 0.0, 2.05009757235e+36,
0.0, -2.09380059113e+38, 0.0, 2.27526964885e+40, 0.0, -2.62577102862e+42, 0.0, 3.21250821027e+44, 0.0, -4.15982781668e+46,
0.0, 5.6920695482e+48, 0.0, -8.21836294198e+50, 0.0, 1.25029043272e+53, 0.0, -2.00155832332e+55, 0.0, 3.36749829154e+57,
0.0, -5.94709705031e+59, 0.0, 1.10119103236e+62, 0.0, -2.13552595453e+64, 0.0, 4.33288969866e+66, 0.0, -9.18855282417e+68,
0.0, 2.03468967763e+71, 0.0, -4.7003833958e+73, 0.0, 1.13180434455e+76, 0.0, -2.83822495707e+78, 0.0, 7.40642489797e+80,
0.0, -2.00964548028e+83, 0.0, 5.66571700508e+85, 0.0, -1.65845111541e+88, 0.0, 5.03688599505e+90, 0.0, -1.58614682377e+93,
0.0, 5.17567436175e+95, 0.0, -1.74889218402e+98, 0.0, 6.1160519995e+100, 0.0, -2.21227769127e+103, 0.0, 8.27227767988e+105,
0.0, -3.19589251114e+108, 0.0, 1.27500822234e+111, 0.0, -5.25009230868e+113, 0.0, 2.23018178942e+116, 0.0, -9.7684521931e+118,
0.0, 4.40983619785e+121, 0.0, -2.05085708865e+124, 0.0, 9.82144332798e+126, 0.0, -4.84126007982e+129, 0.0, 2.45530888015e+132,
0.0, -1.28069268041e+135, 0.0, 6.86761671047e+137, 0.0, -3.7846468582e+140, 0.0, 2.14261012507e+143, 0.0, -1.24567271372e+146,
0.0, 7.4345787551e+148, 0.0, -4.55357953046e+151, 0.0, 2.86121128169e+154, 0.0, -1.8437723552e+157, 0.0, 1.21811545362e+160,
0.0, -8.24821871853e+162, 0.0, 5.72258779378e+165, 0.0, -4.06685305251e+168, 0.0, 2.95960920646e+171, 0.0, -2.20495225652e+174,
0.0, 1.68125970729e+177, 0.0, -1.31167362136e+180, 0.0, 1.04678940095e+183, 0.0, -8.54328935788e+185, 0.0, 7.12878213225e+188,
0.0, -6.08029314555e+191, 0.0, 5.29967764248e+194, 0.0, -4.71942591687e+197, 0.0, 4.29284137914e+200, 0.0, -3.98767449682e+203,
0.0, 3.78197804194e+206, 0.0, -3.66142336837e+209, 0.0, 3.61760902724e+212, 0.0, -3.64707726452e+215, 0.0, 3.75087554365e+218,
0.0, -3.93458672964e+221, 0.0, 4.20882111482e+224, 0.0, -4.59022962206e+227, 0.0, 5.10317257726e+230, 0.0, -5.78227623037e+233,
0.0, 6.67624821678e+236, 0.0, -7.85353076445e+239, 0.0, 9.41068940671e+242, 0.0, -1.14849338735e+246, 0.0, 1.42729587428e+249,
0.0, -1.80595595869e+252, 0.0, 2.32615353077e+255, 0.0, -3.04957517155e+258, 0.0, 4.06858060764e+261, 0.0, -5.5231031322e+264,
0.0, 7.62772793964e+267, 0.0, -1.07155711197e+271, 0.0, 1.5310200896e+274, 0.0, -2.22448916822e+277, 0.0, 3.28626791907e+280,
0.0, -4.9355928956e+283, 0.0, 7.53495712008e+286, 0.0, -1.16914851546e+290, 0.0, 1.84352614678e+293, 0.0, -2.9536826173e+296,
0.0, 4.80793212775e+299, 0.0, -7.95021250459e+302, 0.0, 1.33527841874e+306, 0.0]


Amusez-vous bien!

math_decimal_bernoulli.txt · Dernière modification: 2010/03/16 14:05 par tyrtamos

Outils de la page