# -*- coding: utf-8 -*- def f (L1, rmax, lambd): r""" explication (retourne la matrice dont le coefficient (i,j) est ) INPUT: - bla: EXAMPLES:: sage: f (1, 2, 3) [ e^(-3) 4.0*e^(-3)] sage: f (2, 2, pi) [ e^(-pi) (pi + 1)*e^(-pi)] [ 0.5*e^(-2*pi) 0.25*(2*pi + 1)*e^(-2*pi)] sage: f(0, 2, pi) [] """ assert isinstance(L1, (Integer, int)) and L1 >= 0, 'L1 must be a non negative integer' assert isinstance(rmax, (Integer, int)) and rmax >= 2, 'rmax must be greater than 1' #faire la même chose pour lambda réel positif, mais... attention, le type des réels est... reponse = matrix (L1, rmax, symbolic_expression (0)) t0 = exp (-lambd) t2 = 1 for k in range (L1): t1 = (k + 1) * lambd t2 = t0 * t2 for j in range (rmax): t3 = 1 t4 = 1 for i in range (1, j + 1): t4 = t4 * (j - i + 1) t3 = t1 * t3 + t4 reponse[k, j] = SR(t2*t3*RationalField()((k+1)**(-j-1))) return reponse def Z (r, lambd, L1, F): # r designe une sequence d'entiers non nuls. # lambd est la constante lambda choisie dans l'algorithme de R. Crandall pour écrire Zeta (s) sous forme d'un produit moulien. # L1 est l'ordre de troncature de la série donnant Z. # F désigne la matrice dont l'entré (i, j) est f(i, j, lambda) # n1 et i sont des indices de sommations. # l et q sont des variables, t est un vecteur de variables. j = len (r) q = 0 t = [0 for x in range (j+1)] t[-1] = 1 for n1 in range (j, L1 + 1): q = q + 1 # _ for i in range (j-1, 0, -1): # | Boucle descendante, donné par le pseudo-code de R. Crandall t[i] += t[i + 1] / (q - i + j - 1)**(r[i]) # _| t[0] += F[n1-1,r[0]-1] * t[1] return (t[0] / factorial(r[0] - 1)) def pca (x, y): # x et y désignent deux séquences de même longueur return [sum ([x[j] * y[i-j] for j in range(i + 1)]) for i in range (len (x))] # FFT à utiliser pour accélérer l'algorithme! def pppp (x, y): # x et y désignent deux séquences de même longueur return [x[i] * y[i] for i in range (len(x))] def Y (s, q, lambd, L2, B, G, S): # s désigne la séquence d'entiers sur laquelle on va évaluer le multizêta qui nous intérresse. # q désigne un entier. # lambd est la constante lambda choisi dans l'algorithme de R. Crandall pour écrire Zeta (s) sous forme d'un produit moulien. # L2 est l'ordre de troncature de la série donnant Y. # B désigne les "nombres" de Bernouilli (ou presque...) . # G désigne des quotients de factoriels. # S désigne des sommes de termes de s. # N est un indice de sommation. # l désigne la longueur de la séquence s. # X et Y sont des vecteurs k = len (s) X = B # - for m in range (k - 1): # | X = pppp (X,G[m]) # | Accélération de Bayley due en partie à l'utilisation de la transformation de Fourier rapide X = pca (X,B) # - Y = 0 for N in range (L2 + 1): Y += lambd**(N + S[k-1] + q - k) / (N + S[k-1] + q - k) * X[N] return Y / factorial (s[0] - 1) ; def evalue_zeta_avec_liberte_de_choix_de_lambda_L1_et_L2_intermediaire (s, lambd, L1, L2): # s désigne la séquence d'entiers non vide sur laquelle on va évaluer le multizêta qui nous intérresse. # lambd est la constante lambda choisie dans l'algorithme de R. Crandall pour écrire Zeta (s) sous forme d'un produit moulien. # L1 est l'ordre de troncature de la série donnant Z. # L2 est l'ordre de troncature de la série donnant Y. # k et q sont des indices de sommations. # d, z, t, B, G et S sont des variables : d désigne la longueur de la séquence s. # z désigne la valeur que l'on obtiendra du multizêta en cours de calcul. # B désigne les "nombres" de Bernouilli (ou presque...) . # G désigne des quotients de factoriels. d = len (s) ######################################################### # Calcul des nombres de Bernouilli, une fois pour toute # ######################################################### B = [bernoulli(k) / factorial (k) for k in range (L2 + 1)] ######################################################################################### # Calcul de la matrice carré F contenant à la place (i, j) la valeur de f(i, j, lambda) # ######################################################################################### F = f(L1, max (s), lambd) # Calcul de la matrice F ########################################################### # Calcul des quotients de factoriels, une fois pour toute # ########################################################### S = [s[0]] for i in range (d - 1): S.append (S[i] + s[i+1]) G = matrix (d - 1, L2 + 1, lambda m,n: factorial (n + S[m] - m - 2) / factorial (n + S[m + 1] - m - 2)) ###################################### # Calcul de la somme donnant zeta(s) # ###################################### z = Y(s,0,lambd,L2,B,G,S) z = z + Z(s,lambd,L1,F) # termes isolés de zeta(s) for k in range (1, d): # - t = s # | for q in range (s[k]): # | Somme double dans l'expression de zeta(s) z += (-1)**q / factorial (q) * Y(t[0:k],q,lambd,L2,B,G,S) * Z(t[k:d],lambd,L1,F) # | t[k] += - 1 # - return z def evalue_zeta_avec_liberte_de_choix_de_lambda_L1_et_L2 (s, lambd, L1, L2): # s désigne la séquence d'entiers sur laquelle on va évaluer le multizêta qui nous intérresse. # lambd est la constante lambda choisie dans l'algorithme de R. Crandall pour écrire Zeta (s) sous forme d'un produit moulien. # L1 est l'ordre de troncature de la série donnant Z. # L2 est l'ordre de troncature de la série donnant Y. if s == []: return 1 elif s[0] == 1: raise ValueError, 'Attention ! Multizeta divergent...' else: return evalue_zeta_avec_liberte_de_choix_de_lambda_L1_et_L2_intermediaire (s, lambd, L1, L2) def evalue_zeta (s, d): # s désigne la séquence d'entiers sur laquelle on va évaluer le multizêta qui nous intérresse. # d désigne le nombre de décimales exacts souhaité. # epsilon désigne la précision voulu dans le calcul de zeta(s) lambd = pi / 5 epsilon = 10**(- (d + 1)) L1 = integer_ceil ((d + 1) * ln (10) / lambd) L2 = integer_ceil (- (d + 1) * ln (10) / ln (lambd / (2*pi))) return numerical_approx (evalue_zeta_avec_liberte_de_choix_de_lambda_L1_et_L2 (s, lambd, L1, L2), digits = d)