/* Pour utiliser l'algorithme, sauvegarder ce fichier en le nommant algoCFR (attention, ne pas lui ajouter d'extension !) Puis suivre les indications donnees dans le fichier http://www.math.u-psud.fr/~fischler/modeemploi.pdf ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- Maintenant : en-tete commun a tous les fichiers ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- */ /* dcpelsimples_revolution : prynt*/ /* briques : prwnt */ prynt(aaa1,aaa2,aaa3,aaa4,aaa5,aaa6,aaa7,aaa11,aaa12,aaa13,aaa14,aaa15,aaa16,aaa17) = 0;\ prwnt(aaa1,aaa2,aaa3,aaa4,aaa5,aaa6,aaa7,aaa11,aaa12,aaa13,aaa14,aaa15,aaa16,aaa17) = \ 0; /* prwnt(aaa1,aaa2,aaa3,aaa4,aaa5,aaa6,aaa7,aaa11,aaa12,aaa13,aaa14,aaa15,aaa16,aaa17) = \ print(aaa1,aaa2,aaa3,aaa4,aaa5,aaa6,aaa7,aaa11,aaa12,aaa13,aaa14,aaa15,aaa16,aaa17);\ prynt(aaa1,aaa2,aaa3,aaa4,aaa5,aaa6,aaa7,aaa11,aaa12,aaa13,aaa14,aaa15,aaa16,aaa17) = \ print(aaa1,aaa2,aaa3,aaa4,aaa5,aaa6,aaa7,aaa11,aaa12,aaa13,aaa14,aaa15,aaa16,aaa17); */ /* ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- Maintenant : fichier initialisation_poids19plusclair ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- */ /* Ici on se limite au poids 19 (en fait 99), et a partir du poids 9 on passe a notation moins agreable */ /*Mais des le poids 2 on separe les differents MZV de meme poids*/ /* On suppose aussi prof <= 4 pour ce fichier (passage de 3 a 4 le 1er fev. 2005, non teste) (et <= 6, pour zetalapos; et prof <= 9 pour le reste des fichiers */ /*Il y aura des pb si les MZV a DV log de poids >= 7 apparaissent. Mais probablement ce ne sera pas le cas ; et si c'est le cas on pourra toujours implementer de maniere automatique la regularisation shuffle*/ /* De plus attention, on passera peut-etre a cote de certaines annulations, dues a des relations lineaires en poids >= 7*/ /*Si on doit passer de wtron=9 a qqch de plus grand, utiliser Pari (ou le shell...) pour produire un fichier texte contenant le tableau ldzptr, puis le lire par la commande read....*/ /* tableau se refere a (cf. manuscrit p. 47 et 87, avec correction p. 51) [1, zeta(2), zeta(3), zeta(4), zeta(5), zeta(4,1), zeta(6), zeta(5,1), le_reste] qui comporte 9 cases ; ici zeta(.) est le MZV habituel, et le_reste est en fait tout... */ wmax = 99; wtron = 9; /*mais en fait on ne tronque pas en poids 9, on passe juste a une notation moins agreable*/ ldzpun = vector(wtron); ldzpun[1] = z_1; ldzpun[2] = z_2; ldzpun[3] = z_3; ldzpun[4] = z_4; ldzpun[5] = z_5; ldzpun[6] = z_6; ldzpun[7] = z_7; ldzpun[8] = z_8; ldzpun[9] = z_9; ldzpde = matrix(wtron,wtron); /* profondeur 2*/ /*On ne definit ldzpde[a,b] = z_{ab} que pour 0 < a+b <= wtron, les autres sont tronqu\'es*/ /*avec a,b >= 1 car les DV log sont consideres comme CV*/ ldzpde[1,1] = z_11; ldzpde[1,2] = z_12; ldzpde[1,3] = z_13; ldzpde[1,4] = z_14; ldzpde[1,5] = z_15; ldzpde[1,6] = z_16; ldzpde[1,7] = z_17; ldzpde[1,8] = z_18; ldzpde[2,1] = z_21; ldzpde[2,2] = z_22; ldzpde[2,3] = z_23; ldzpde[2,4] = z_24; ldzpde[2,5] = z_25; ldzpde[2,6] = z_26; ldzpde[2,7] = z_27; ldzpde[3,1] = z_31; ldzpde[3,2] = z_32; ldzpde[3,3] = z_33; ldzpde[3,4] = z_34; ldzpde[3,5] = z_35; ldzpde[3,6] = z_36; ldzpde[4,1] = z_41; ldzpde[4,2] = z_42; ldzpde[4,3] = z_43; ldzpde[4,4] = z_44; ldzpde[4,5] = z_45; ldzpde[5,1] = z_51; ldzpde[5,2] = z_52; ldzpde[5,3] = z_53; ldzpde[5,4] = z_54; ldzpde[6,1] = z_61; ldzpde[6,2] = z_62; ldzpde[6,3] = z_63; ldzpde[7,1] = z_71; ldzpde[7,2] = z_72; ldzpde[8,1] = z_81; ldzptr = matrix(wtron,wtron); /* profondeur 3*/ /*On ne definit ldzpde[a,b][c] = z_{abc} que pour 0 < a+b+c <= wtron, les autres sont tronqu\'es*/ /*avec a,b,c >= 1 car les DV log sont consideres comme CV*/ /*On suppose donc a+b <= 8*/ for(iptrun = 1, wtron,\ for(iptrde = 1, wtron,\ ldzptr[iptrun, iptrde] = vector(wtron);\ );\ ); ldzptr[1,1][1] = z_111; ldzptr[1,1][2] = z_112; ldzptr[1,1][3] = z_113; ldzptr[1,1][4] = z_114; ldzptr[1,1][5] = z_115; ldzptr[1,1][6] = z_116; ldzptr[1,1][7] = z_117; ldzptr[1,2][1] = z_121; ldzptr[1,2][2] = z_122; ldzptr[1,2][3] = z_123; ldzptr[1,2][4] = z_124; ldzptr[1,2][5] = z_125; ldzptr[1,2][6] = z_126; ldzptr[1,3][1] = z_131; ldzptr[1,3][2] = z_132; ldzptr[1,3][3] = z_133; ldzptr[1,3][4] = z_134; ldzptr[1,3][5] = z_135; ldzptr[1,4][1] = z_141; ldzptr[1,4][2] = z_142; ldzptr[1,4][3] = z_143; ldzptr[1,4][4] = z_144; ldzptr[1,5][1] = z_151; ldzptr[1,5][2] = z_152; ldzptr[1,5][3] = z_153; ldzptr[1,6][1] = z_161; ldzptr[1,6][2] = z_162; ldzptr[1,7][1] = z_171; ldzptr[2,1][1] = z_211; ldzptr[2,1][2] = z_212; ldzptr[2,1][3] = z_213; ldzptr[2,1][4] = z_214; ldzptr[2,1][5] = z_215; ldzptr[2,1][6] = z_216; ldzptr[2,2][1] = z_221; ldzptr[2,2][2] = z_222; ldzptr[2,2][3] = z_223; ldzptr[2,2][4] = z_224; ldzptr[2,2][5] = z_225; ldzptr[2,3][1] = z_231; ldzptr[2,3][2] = z_232; ldzptr[2,3][3] = z_233; ldzptr[2,3][4] = z_234; ldzptr[2,4][1] = z_241; ldzptr[2,4][2] = z_242; ldzptr[2,4][3] = z_243; ldzptr[2,5][1] = z_251; ldzptr[2,5][2] = z_252; ldzptr[2,6][1] = z_261; ldzptr[3,1][1] = z_311; ldzptr[3,1][2] = z_312; ldzptr[3,1][3] = z_313; ldzptr[3,1][4] = z_314; ldzptr[3,1][5] = z_315; ldzptr[3,2][1] = z_321; ldzptr[3,2][2] = z_322; ldzptr[3,2][3] = z_323; ldzptr[3,2][4] = z_324; ldzptr[3,3][1] = z_331; ldzptr[3,3][2] = z_332; ldzptr[3,3][3] = z_333; ldzptr[3,4][1] = z_341; ldzptr[3,4][2] = z_342; ldzptr[3,5][1] = z_351; ldzptr[4,1][1] = z_411; ldzptr[4,1][2] = z_412; ldzptr[4,1][3] = z_413; ldzptr[4,1][4] = z_414; ldzptr[4,2][1] = z_421; ldzptr[4,2][2] = z_422; ldzptr[4,2][3] = z_423; ldzptr[4,3][1] = z_431; ldzptr[4,3][2] = z_432; ldzptr[4,4][1] = z_441; ldzptr[5,1][1] = z_511; ldzptr[5,1][2] = z_512; ldzptr[5,1][3] = z_513; ldzptr[5,2][1] = z_521; ldzptr[5,2][2] = z_522; ldzptr[5,3][1] = z_531; ldzptr[6,1][1] = z_611; ldzptr[6,1][2] = z_612; ldzptr[6,2][1] = z_621; ldzptr[7,1][1] = z_711; /*ldzptr et ses potes ne sont implementes qu'en prof <= 3; en prof >= 4, on passe directement a la notation avec valeurzegros, meme en petit poids*/ ze; /*Le but est que ze apparaisse apres les z_{ijk} dans les listes de valeurs zeta*/ zf; zg; /*Attention, en cas de prof 3 ou 4 il peut y avoir une ambiguite ; mais la formulation ancienne provoquait un degree overflow (cf. eventuellement doc de Pari pour trouver une autre solution) Mais le pb est seulement si s[2] ou s[3] est >= 10 : reste rare a priori.... */ valeurzegros(N,s) = \ local(resuvaleurzegros);\ if(N == 1, resuvaleurzegros = ze^(s[1]));\ if(N == 2, resuvaleurzegros = ze^(100*s[1] + s[2]));\ if(N == 3, resuvaleurzegros = zf^(100*s[1] + 10*s[2] + s[3]));\ if(N == 4, resuvaleurzegros = zg^(1000*s[1] + 100*s[2] + 10*s[3] + s[4]));\ resuvaleurzegros valeurconstante(valeurlamb) = [valeurlamb,0,0, 0,0,0, 0,0,0]; valeurgrospoids(valeurgp) = [0,0,0, 0,0,0, 0,0,valeurgp]; xvar_1 = 'xvar_1;\ xvar_2 = 'xvar_2;\ xvar_3 = 'xvar_3;\ xvar_4 = 'xvar_4;\ xvar_5 = 'xvar_5;\ xvar_6 = 'xvar_6;\ xvar_7 = 'xvar_7;\ xvar_8 = 'xvar_8;\ xvar_9 = 'xvar_9;\ ldvar = [xvar_1,xvar_2,xvar_3,xvar_4,xvar_5,xvar_6,xvar_7,xvar_8,xvar_9]; /*et ensuite il faut tout specialiser a 1, a la fin, pour donner un sens a l'ensemble */ xfige_1 = 'xfige_1;\ xfige_2 = 'xfige_2;\ xfige_3 = 'xfige_3;\ xfige_4 = 'xfige_4;\ xfige_5 = 'xfige_5;\ xfige_6 = 'xfige_6;\ xfige_7 = 'xfige_7;\ xfige_8 = 'xfige_8;\ xfige_9 = 'xfige_9;\ varfige = [xfige_1,xfige_2,xfige_3,xfige_4,xfige_5,xfige_6,\ xfige_7,xfige_8,xfige_9]; /*un zetalarge est la valeur en 1 d'un polylog large*/ /*un zetastrict est un MZV usuel*/ /* On utilise que zeta(3) = zeta(2,1) et les autres relations du meme style */ /*Attention, xlistecv ne doit apparaitre nulle part */ /*Utilisation des relations lineaires entre zetas CV en poids <= 6: manuscrit p. 87 et 88 pour le poids 6*/ /* Pour le poids >= 7 et <= wtron, on considere les zetas a DV log comme etant CV, i.e. dans ce cas zetastrictpositif renvoie sur zetastrictconvergent directement. */ /*zetastrictsimplifie utilise les relations entre MZV en poids <= 6 ; ecrit un tel MZV comme un tableau de coord sur les MZV principaux.*/ zetastrictsimplifie(N,s,xlistesimp) = \ local(resulasimp, poidssimp);\ if(N == 0, resulasimp = [1,0,0, 0,0,0, 0,0,0]);\ if(N == 1, \ if(s[1] == 2, resulasimp = [0,1,0, 0,0,0, 0,0,0]);\ if(s[1] == 3, resulasimp = [0,0,1, 0,0,0, 0,0,0]);\ if(s[1] == 4, resulasimp = [0,0,0, 1,0,0, 0,0,0]);\ if(s[1] == 5, resulasimp = [0,0,0, 0,1,0, 0,0,0]);\ if(s[1] == 6, resulasimp = [0,0,0, 0,0,0, 1,0,0]);\ );\ if((N == 2) && (s[1] == 2),\ if(s[2] == 1, resulasimp = zetastrictsimplifie(1,[3],xlistesimp));\ if(s[2] == 2, resulasimp = (3/4)*zetastrictsimplifie(1,[4],xlistesimp));\ if(s[2] == 3, resulasimp = (1/2)*zetastrictsimplifie(1,[5],xlistesimp) + \ 2*zetastrictsimplifie(2,[4,1],xlistesimp));\ if(s[2] == 4, resulasimp = (7/12)*zetastrictsimplifie(1,[6],xlistesimp) + \ 2*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if((N == 2) && (s[1] == 3),\ if(s[2] == 1, resulasimp = (1/4)*zetastrictsimplifie(1,[4],xlistesimp));\ if(s[2] == 2, resulasimp = (1/2)*zetastrictsimplifie(1,[5],xlistesimp) - \ 3*zetastrictsimplifie(2,[4,1],xlistesimp));\ if(s[2] == 3, resulasimp = (1/4)*zetastrictsimplifie(1,[6],xlistesimp) + \ (-1)*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if((N == 2) && (s[1] == 4),\ if(s[2] == 1, resulasimp = [0,0,0, 0,0,1, 0,0,0]);\ if(s[2] == 2, resulasimp = (1/6)*zetastrictsimplifie(1,[6],xlistesimp) + \ (-2)*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if((N == 2) && (s[1] == 5),\ if(s[2] == 1, resulasimp = [0,0,0, 0,0,0, 0,1,0]);\ );\ if((N == 3) && (s[1] == 2),\ if(s[2] == 1,\ if(s[3] == 1, resulasimp = zetastrictsimplifie(1,[4],xlistesimp));\ if(s[3] == 2, resulasimp = (1/2)*zetastrictsimplifie(1,[5],xlistesimp) + \ 2*zetastrictsimplifie(2,[4,1],xlistesimp));\ if(s[3] == 3, resulasimp = (11/16)*zetastrictsimplifie(1,[6],xlistesimp) + \ (-2)*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if(s[2] == 2,\ if(s[3] == 1, resulasimp = \ (1/2)*zetastrictsimplifie(1,[5],xlistesimp) - \ 3*zetastrictsimplifie(2,[4,1],xlistesimp));\ if(s[3] == 2, resulasimp = (3/16)*zetastrictsimplifie(1,[6],xlistesimp));\ );\ if(s[2] == 3,\ if(s[3] == 1, resulasimp = (-1/24)*zetastrictsimplifie(1,[6],xlistesimp)+\ 3*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ );\ if((N == 3) && (s[1] == 3),\ if((s[2] == 1) & (s[3] == 1), \ resulasimp = zetastrictsimplifie(2,[4,1],xlistesimp));\ if((s[2] == 1) & (s[3] == 2), \ resulasimp = (-1/24)*zetastrictsimplifie(1,[6],xlistesimp) + \ 3*zetastrictsimplifie(2,[5,1],xlistesimp));\ if((s[2] == 2) & (s[3] == 1),\ resulasimp = (13/48)*zetastrictsimplifie(1,[6],xlistesimp) + \ (-6)*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if((N == 3) && (s[1] == 4),\ if((s[2] == 1) & (s[3] == 1), \ resulasimp = (-1/16)*zetastrictsimplifie(1,[6],xlistesimp) + \ 2*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if((N == 4) && (s[1] == 2),\ if((s[2] == 1) & (s[3] == 1) & (s[4] == 1), \ resulasimp = zetastrictsimplifie(1,[5],xlistesimp));\ if((s[2] == 1) & (s[3] == 1) & (s[4] == 2), \ resulasimp = (7/12)*zetastrictsimplifie(1,[6],xlistesimp) + \ 2*zetastrictsimplifie(2,[5,1],xlistesimp));\ if((s[2] == 1) & (s[3] == 2) & (s[4] == 1), \ resulasimp = (1/4)*zetastrictsimplifie(1,[6],xlistesimp) + \ (-1)*zetastrictsimplifie(2,[5,1],xlistesimp));\ if((s[2] == 2) & (s[3] == 1) & (s[4] == 1), \ resulasimp = (1/6)*zetastrictsimplifie(1,[6],xlistesimp) + \ (-2)*zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if((N == 4) && (s[1] == 3),\ if((s[2] == 1) & (s[3] == 1) & (s[4] == 1), \ resulasimp = zetastrictsimplifie(2,[5,1],xlistesimp));\ );\ if((N == 5) && (s[1] == 2),\ if((s[2] == 1) & (s[3] == 1) & (s[4] == 1) & (s[5] == 1), \ resulasimp = zetastrictsimplifie(1,[6],xlistesimp));\ );\ poidssimp = sum(indlacv=1,N,s[indlacv]);\ if(poidssimp > 6, resulasimp = horschamp; print("on demande grosse simplif. initpoids6regsimp"); 0/0);\ resulasimp zetastrictconvergent(N,s,xlistecv) = \ local(resulacv, poidscv);\ poidscv = sum(indlacv=1,N,s[indlacv]);\ if(poidscv > wmax, resulacv = horschamp; print("Probleme initpoids6reg"); 0/0);\ if((poidscv <= wmax) && (poidscv > wtron), resulacv = valeurgrospoids(valeurzegros(N,s)));\ if((poidscv <= wtron) && (poidscv > 6), \ if(N == 1, resulacv = valeurgrospoids(ldzpun[s[1]]));\ if(N == 2, resulacv = valeurgrospoids(ldzpde[s[1],s[2]]));\ if(N == 3, resulacv = valeurgrospoids(ldzptr[s[1],s[2]][s[3]]));\ if(N == 4, resulacv = valeurgrospoids(valeurzegros(N,s)));\ if(N >= 5, print("Non implemente en prof >= 5 init"); 0/0;);\ );\ if(poidscv <= 6,\ if(N == 1, resulacv = zetastrictsimplifie(N,s,xlistecv) + valeurgrospoids(ldzpun[s[1]]));\ if(N == 2, resulacv = zetastrictsimplifie(N,s,xlistecv) + valeurgrospoids(ldzpde[s[1],s[2]]));\ if(N == 3, resulacv = zetastrictsimplifie(N,s,xlistecv) + valeurgrospoids(ldzptr[s[1],s[2]][s[3]]));\ if(N == 4, resulacv = zetastrictsimplifie(N,s,xlistecv) + valeurgrospoids(valeurzegros(N,s)));\ if(N >= 5, print("Non implemente en prof >= 5 init"); 0/0;);\ );\ resulacv /* On regularise les DV comme indique pages 18 et 50 des notes manuscrites, et pages 79 \`a 81 pour le poids 5 et pages 82 \`a 85 pour le poids 6 */ /*Probablement xlistepos ne sert a rien ici */ /*Utilisation des formules pour regulariser des zetas DV positifs*/ zetastrictpositif(N,s,xlistepos) = \ local(resulapos, poidspos);\ poidspos = sum(ipo = 1,N,s[ipo]);\ if(poidspos > wmax,\ print("On utilise le poids ",sum(ipo = 1,N,s[ipo]), \ "alors que ce n'est pas impl\'ement\'e ; pb dans initp6reg"); 0/0; \ );\ if(poidspos <= 4,\ if((N == 0) || (s[1] >= 2), \ resulapos = zetastrictconvergent(N,s,xlistepos);\ , \ if(N == 1, resulapos = valeurconstante(0));\ if(N == 2, \ if(s[2] == 1, resulapos = valeurconstante(0));\ if(s[2] == 2, resulapos = -2*zetastrictconvergent(1,[3],xlistepos));\ if(s[2] == 3, resulapos = -zetastrictconvergent(2,[2,2],xlistepos) -\ 2*zetastrictconvergent(2,[3,1],xlistepos));\ );\ if((N == 3) & (s[2] == 1), \ if(s[3] == 1, resulapos = valeurconstante(0));\ if(s[3] == 2, resulapos = 3*zetastrictconvergent(3,[2,1,1],xlistepos));\ );\ if((N == 3) & (s[2] == 2), \ if(s[3] == 1, resulapos = -3*zetastrictconvergent(3,[2,1,1],xlistepos));\ );\ if((N == 4) & (s[2] == 1) & (s[3] == 1) & (s[4]==1), \ resulapos = valeurconstante(0));\ if(sum(indlapos=1,N,s[indlapos]) > 4, resulapos = horschamp ; print("Pb initp6reg_a"); 0/0);\ );\ ,\ if(poidspos == 5,\ if(s[1] == 1,\ if(s[2] == 4,\ resulapos = -zetastrictconvergent(2,[2,3],xlistepos)\ -zetastrictconvergent(2,[3,2],xlistepos)\ -2*zetastrictconvergent(2,[4,1],xlistepos);\ );\ if(s[2] == 3,\ resulapos = -zetastrictconvergent(3,[2,2,1],xlistepos)\ -3*zetastrictconvergent(3,[3,1,1],xlistepos);\ );\ if(s[2] == 2,\ if(s[3] == 2,\ resulapos = -2*zetastrictconvergent(3,[2,1,2],xlistepos)\ -2*zetastrictconvergent(3,[2,2,1],xlistepos);\ );\ if(s[3] == 1,\ resulapos = -4*zetastrictconvergent(4,[2,1,1,1],xlistepos);\ );\ );\ if(s[2] == 1,\ if(s[3] == 3,\ resulapos = zetastrictconvergent(3,[2,1,2],xlistepos)\ +2*zetastrictconvergent(3,[2,2,1],xlistepos)\ +3*zetastrictconvergent(3,[3,1,1],xlistepos);\ );\ if(s[3] == 2,\ resulapos = 6*zetastrictconvergent(4,[2,1,1,1],xlistepos);\ );\ if(s[3] == 1,\ if(s[4] == 2,\ resulapos = -4*zetastrictconvergent(4,[2,1,1,1],xlistepos);\ );\ if(s[4] == 1,\ resulapos = valeurconstante(0);\ );\ );\ );\ );\ if(s[1] >= 2,\ resulapos = zetastrictconvergent(N,s,xlistepos);\ );\ );\ if(poidspos == 6,\ if(s[1] == 1,\ if(s[2] == 5,\ resulapos = -zetastrictconvergent(2,[2,4],xlistepos)\ -zetastrictconvergent(2,[3,3],xlistepos)\ -zetastrictconvergent(2,[4,2],xlistepos)\ -2*zetastrictconvergent(2,[5,1],xlistepos);\ );\ if(s[2] == 4,\ resulapos = -zetastrictconvergent(3,[2,3,1],xlistepos)\ -zetastrictconvergent(3,[3,2,1],xlistepos)\ -3*zetastrictconvergent(3,[4,1,1],xlistepos);\ );\ if(s[2] == 3,\ if(s[3] == 2,\ resulapos = -zetastrictconvergent(3,[2,2,2],xlistepos)\ -2*zetastrictconvergent(3,[3,1,2],xlistepos)\ -2*zetastrictconvergent(3,[3,2,1],xlistepos);\ );\ if(s[3] == 1,\ resulapos = -zetastrictconvergent(4,[2,2,1,1],xlistepos)\ -4*zetastrictconvergent(4,[3,1,1,1],xlistepos);\ );\ );\ if(s[2] == 2,\ if(s[3] == 3,\ resulapos =-2*zetastrictconvergent(3,[2,1,3],xlistepos)\ -zetastrictconvergent(3,[2,2,2],xlistepos)\ -2*zetastrictconvergent(3,[2,3,1],xlistepos);\ );\ if(s[3] == 2,\ resulapos = -2*zetastrictconvergent(4,[2,1,2,1],xlistepos)\ -3*zetastrictconvergent(4,[2,2,1,1],xlistepos);\ );\ if(s[3] == 1,\ if(s[4] == 2,\ resulapos = -3*zetastrictconvergent(4,[2,1,1,2],xlistepos)\ -2*zetastrictconvergent(4,[2,1,2,1],xlistepos);\ );\ if(s[4] == 1,\ resulapos = -5*zetastrictconvergent(5,[2,1,1,1,1],xlistepos);\ );\ );\ );\ if(s[2] == 1,\ if(s[3] == 4,\ resulapos = zetastrictconvergent(3,[2,1,3],xlistepos)\ +zetastrictconvergent(3,[2,3,1],xlistepos)\ +zetastrictconvergent(3,[3,1,2],xlistepos)\ +zetastrictconvergent(3,[3,2,1],xlistepos)\ +zetastrictconvergent(3,[2,2,2],xlistepos)\ -zetastrictconvergent(3,[1,4,1],xlistepos);\ );\ if(s[3] == 3,\ resulapos = zetastrictconvergent(4,[2,1,2,1],xlistepos)\ +3*zetastrictconvergent(4,[2,2,1,1],xlistepos)\ +6*zetastrictconvergent(4,[3,1,1,1],xlistepos);\ );\ if(s[3] == 2,\ if(s[4] == 2,\ resulapos = 3*zetastrictconvergent(4,[2,1,1,2],xlistepos)\ +4*zetastrictconvergent(4,[2,1,2,1],xlistepos)\ +3*zetastrictconvergent(4,[2,2,1,1],xlistepos);\ );\ if(s[4] == 1,\ resulapos = 10*zetastrictconvergent(5,[2,1,1,1,1],xlistepos);\ );\ );\ if(s[3] == 1,\ if(s[4] == 3,\ resulapos = -zetastrictconvergent(4,[2,1,1,2],xlistepos)\ -2*zetastrictconvergent(4,[2,1,2,1],xlistepos)\ -3*zetastrictconvergent(4,[2,2,1,1],xlistepos)\ -4*zetastrictconvergent(4,[3,1,1,1],xlistepos);\ );\ if(s[4] == 2,\ resulapos = -10*zetastrictconvergent(5,[2,1,1,1,1],xlistepos);\ );\ if((s[4] == 1) & (s[5] == 2),\ resulapos = 5*zetastrictconvergent(5,[2,1,1,1,1],xlistepos);\ );\ if((s[4] == 1) & (s[5] == 1),\ resulapos = valeurconstante(0);\ );\ );\ );\ );\ if(s[1] >= 2,\ resulapos = zetastrictconvergent(N,s,xlistepos);\ );\ );\ if((poidspos > 6) && (poidspos <= wmax), resulapos = zetastrictconvergent(N,s,xlistepos));\ );\ resulapos /* Calcul des zetalargepositif en fonction des zetastrictpositif*/ /*Il s'agit toujours de valeurs regularisees*/ /* par def, un zetalargepositif est tq les s_i soient tous >= 1 */ /*Le meme algo s'appliquerait pour dcp des polylogs larges positifs*/ zetalargepositif(N,s,xlistepos) = \ local(nvx,nvN,nvs,resullap);\ resullap = 0;\ \ if(N == 1, resullap = zetastrictpositif(1,s,xlistepos));\ if(N == 2,\ for(ilap_1 = 0, 1,\ nvx = xlistepos;\ nvs = s;\ nvN = 1;\ if(ilap_1 == 0,\ nvN = nvN + 1;\ );\ if(ilap_1 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[2];\ nvs[nvN] = nvs[nvN]+s[2];\ );\ resullap = resullap + zetastrictpositif(nvN,nvs,nvx);\ );\ );\ \ if(N == 3,\ for(ilap_1 = 0, 1,\ for(ilap_2 = 0, 1,\ nvx = xlistepos;\ nvs = s;\ nvN = 1;\ if(ilap_1 == 0,\ nvN = nvN + 1;\ );\ if(ilap_1 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[2];\ nvs[nvN] = nvs[nvN]+s[2];\ );\ if(ilap_2 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[3];\ nvs[nvN] = s[3];\ );\ if(ilap_2 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[3];\ nvs[nvN] = nvs[nvN]+s[3];\ );\ resullap = resullap + zetastrictpositif(nvN,nvs,nvx);\ );\ );\ );\ \ if(N == 4,\ for(ilap_1 = 0, 1,\ for(ilap_2 = 0, 1,\ for(ilap_3 = 0, 1,\ nvx = xlistepos;\ nvs = s;\ nvN = 1;\ if(ilap_1 == 0,\ nvN = nvN + 1;\ );\ if(ilap_1 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[2];\ nvs[nvN] = nvs[nvN]+s[2];\ );\ if(ilap_2 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[3];\ nvs[nvN] = s[3];\ );\ if(ilap_2 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[3];\ nvs[nvN] = nvs[nvN]+s[3];\ );\ if(ilap_3 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[4];\ nvs[nvN] = s[4];\ );\ if(ilap_3 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[4];\ nvs[nvN] = nvs[nvN]+s[4];\ );\ resullap = resullap + zetastrictpositif(nvN,nvs,nvx);\ );\ );\ );\ );\ \ if(N == 5,\ for(ilap_1 = 0, 1,\ for(ilap_2 = 0, 1,\ for(ilap_3 = 0, 1,\ for(ilap_4 = 0, 1,\ nvx = xlistepos;\ nvs = s;\ nvN = 1;\ if(ilap_1 == 0,\ nvN = nvN + 1;\ );\ if(ilap_1 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[2];\ nvs[nvN] = nvs[nvN]+s[2];\ );\ if(ilap_2 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[3];\ nvs[nvN] = s[3];\ );\ if(ilap_2 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[3];\ nvs[nvN] = nvs[nvN]+s[3];\ );\ if(ilap_3 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[4];\ nvs[nvN] = s[4];\ );\ if(ilap_3 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[4];\ nvs[nvN] = nvs[nvN]+s[4];\ );\ if(ilap_4 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[5];\ nvs[nvN] = s[5];\ );\ if(ilap_4 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[5];\ nvs[nvN] = nvs[nvN]+s[5];\ );\ resullap = resullap + zetastrictpositif(nvN,nvs,nvx);\ );\ );\ );\ );\ );\ \ \ if(N == 6,\ for(ilap_1 = 0, 1,\ for(ilap_2 = 0, 1,\ for(ilap_3 = 0, 1,\ for(ilap_4 = 0, 1,\ for(ilap_5 = 0, 1,\ nvx = xlistepos;\ nvs = s;\ nvN = 1;\ if(ilap_1 == 0,\ nvN = nvN + 1;\ );\ if(ilap_1 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[2];\ nvs[nvN] = nvs[nvN]+s[2];\ );\ if(ilap_2 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[3];\ nvs[nvN] = s[3];\ );\ if(ilap_2 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[3];\ nvs[nvN] = nvs[nvN]+s[3];\ );\ if(ilap_3 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[4];\ nvs[nvN] = s[4];\ );\ if(ilap_3 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[4];\ nvs[nvN] = nvs[nvN]+s[4];\ );\ if(ilap_4 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[5];\ nvs[nvN] = s[5];\ );\ if(ilap_4 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[5];\ nvs[nvN] = nvs[nvN]+s[5];\ );\ if(ilap_5 == 0,\ nvN = nvN + 1;\ nvx[nvN] = xlistepos[6];\ nvs[nvN] = s[6];\ );\ if(ilap_5 == 1,\ nvx[nvN] = nvx[nvN]*xlistepos[6];\ nvs[nvN] = nvs[nvN]+s[6];\ );\ resullap = resullap + zetastrictpositif(nvN,nvs,nvx);\ );\ );\ );\ );\ );\ );\ \ resullap lapositif(N,s,xlapos) = \ zetalargepositif(N,s,xlapos) /* ------------------------------------------------------------- Calcul de la vraie valeur (numerique) d'un tableau dont les cases correspondent aux coeffs de [1, zeta(2), zeta(3), zeta(4), zeta(5), zeta(4,1), zeta(6), zeta(5,1)] ------------------------------------------------------------- */ tableaudesvraiszetas = vector(8); tableaudesvraiszetas[1] = 1; tableaudesvraiszetas[2] = zeta(2); tableaudesvraiszetas[3] = zeta(3); tableaudesvraiszetas[4] = zeta(4); tableaudesvraiszetas[5] = zeta(5); tableaudesvraiszetas[6] = 0.09655115998944373446564553142894276403201; /*Valeur de zeta(4,1) donnee par EZ-face*/ tableaudesvraiszetas[7] = zeta(6); tableaudesvraiszetas[8] = 0.04053689727151973782904590793969648233449; /*Valeur de zeta(5,1) donnee par EZ-face*/ vraievaleurduntab(tableauvv) = \ resulvvtab = 0;\ for(ivvtab = 1, length(valeurconstante(0)),\ resulvvtab = resulvvtab + tableauvv[ivvtab] * tableaudesvraiszetas[ivvtab];\ );\ xvar_1 = 1;\ xvar_2 = 1;\ xvar_3 = 1;\ xvar_4 = 1;\ xvar_5 = 1;\ xvar_6 = 1;\ xvar_7 = 1;\ xvar_8 = 1;\ xvar_9 = 1;\ resulvvtab = eval(resulvvtab);\ xvar_1 = 'xvar_1;\ xvar_2 = 'xvar_2;\ xvar_3 = 'xvar_3;\ xvar_4 = 'xvar_4;\ xvar_5 = 'xvar_5;\ xvar_6 = 'xvar_6;\ xvar_7 = 'xvar_7;\ xvar_8 = 'xvar_8;\ xvar_9 = 'xvar_9;\ resulvvtab /* ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- Maintenant : fichier polylogs ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- */ Smax = 20; /*borne quoi exactement ? */ /*Ce programme calcule les la a partir des lapositif */ /*Ce programme ne d\'epend pas de wmax */ /*lanegatifunevar est La_{s[1]}(x_1) avec s[1] <= 0 */ /*donc ici s <= 0 ; ne prend pas de variable N car c'est 1*/ /*ici s est un entier, et PAS un tableau de taille N=1 */ /*attention lanegatifunevar renvoie une constante, et pas un tableau donc parfois necessaire appeler valeurconstante(lanegatifunevar) */ lanegatifunevar(s,xneg) = \ -subst(taba[1,-s+1],xlibre,xneg); /*Polynome P_s(k,x) de TR, p. 13 ; avec s >= 0 et k >= 0 ; par def P_s(k,x) = sum_{l=1} ^k l^s x^(-l) On a P_s(k,x) = x^(-k) A_{x,s}(k) - a_{0,x,s} pour x \neq 0,1 o\`u A_{x,s}(k) = sum_{t=0} ^s a_{t,x,s} k^t est un polynome en k de degre s, dont les coeffs sont des fractions rationnelles en x avec 1 pour seul pole. En x=1 l'evaluation de P_s(k,1) se ferait avec compensation des poles, c'est un probleme qu'on n'a pas besoin de traiter dans l'algo. */ /*Implementation : taba est une matrice carree de taille Smax+1 dont les coeffs sont taba[tbi,sbi] = a_{tbi-1,x,sbi-1} pour 1 \leq tbi \leq sbi qui sont des fracn ratnelles en xlibre */ /*taba est une variable globale, ou quasi... */ calculdetaba() = \ taba = matrix(Smax+1,Smax+1);\ taba[1,1] = 1/(1-xlibre) ;\ for(sbi = 1, Smax,\ taba[1,sbi+1] = - xlibre * deriv(taba[1,sbi],xlibre);\ for(tbi = 2, sbi,\ taba[tbi,sbi+1] = taba[tbi-1,sbi] - xlibre * deriv(taba[tbi,sbi],xlibre);\ );\ taba[sbi+1,sbi+1] = taba[sbi,sbi];\ );\ calculdetaba(); /*les la sont des fonctions de x_1,x_2,x_3 ; elles sont constantes ssi tous les s_i sont >= 1 */ la(N,s,xliste) =\ local(j, resula,nvxliste,nvsliste,decalagesuites);\ if(N <= 1, \ if((N ==0) || (s[1] >= 1), \ resula = lapositif(N,s,xliste);\ ,\ resula = valeurconstante(lanegatifunevar(s[1],xliste[1]));\ );\ , \ j = N; \ while((j >= 1) && (s[j] >= 1), j = j-1);\ if(j == 0, resula = lapositif(N,s,xliste));\ if(j == N, \ print("On utilise utilise un La <= 0 : pas normal ! (cas j=N)polylogs"); 0/0; \ nvxliste = xliste;\ nvxliste[N-1] = xliste[N-1] * xliste[N];\ decalagesuites = vector(length(s));\ decalagesuites[N-1] = 1;\ resula = - subst(taba[1,-s[N]+1],xlibre,xliste[N])*la(N-1,s,xliste)+ \ sum(t=0, -s[N], \ subst(taba[t+1,-s[N]+1],xlibre,xliste[N]) * \ la(N-1,s - t*decalagesuites,nvxliste)\ );\ );\ if((j >= 2) && (j <= N-1),\ print("On utilise utilise un La <= 0 : pas normal ! (cas 2 <= j <= N-1)polylogs"); 0/0; \ nvxliste = xliste;\ nvxliste[j-1] = xliste[j-1] * xliste[j];\ for(jnv = j, N-1, nvxliste[jnv] = xliste[jnv+1]);\ nvsliste = vector(N-1);\ for(jnv = 1, j-1, nvsliste[jnv] = s[jnv]);\ for(jnv = j, N-1, nvsliste[jnv] = s[jnv+1]);\ decalagesuites = vector(N-1);\ decalagesuites[j-1] = 1;\ resula = sum(t=0, -s[j], \ subst(taba[t+1,-s[j]+1],xlibre,xliste[j]) * \ la(N-1,nvsliste - t*decalagesuites,nvxliste)\ );\ nvxliste = xliste;\ nvxliste[j] = xliste[j] * xliste[j+1];\ for(jnv = j+1, N-1, nvxliste[jnv] = xliste[jnv+1]);\ /*nvsliste est la meme que juste au-dessus*/ \ decalagesuites = vector(N-1);\ decalagesuites[j] = 1;\ resula = resula - sum(tpr=0, -s[j], \ sum(tsec=tpr,-s[j],(-1)^(tsec-tpr)*binomial(tsec,tpr)*xliste[j]*\ subst(taba[tsec+1,-s[j]+1],xlibre,xliste[j])) * \ la(N-1,nvsliste - tpr*decalagesuites,nvxliste)\ );\ );\ if(j == 1,\ print("On utilise utilise un La <= 0 : pas normal !polylogs",s); 0/0; \ nvxliste = vector(N-1);\ for(jnv = 1, N-1, nvxliste[jnv] = xliste[jnv+1]);\ nvsliste = vector(N-1);\ for(jnv = 1, N-1, nvsliste[jnv] = s[jnv+1]);\ resula = (lanegatifunevar(s[1],xliste[1]) + \ subst(taba[1,-s[1]+1],xlibre,xliste[1])) * \ la(N-1,nvsliste,nvxliste);\ nvxliste = xliste;\ nvxliste[1] = xliste[1] * xliste[2];\ for(jnv = 2, N-1, nvxliste[jnv] = xliste[jnv+1]);\ /*nvsliste est la meme que juste au-dessus*/ \ decalagesuites = vector(N-1);\ decalagesuites[1] = 1;\ resula = resula - sum(tpr=0, -s[1],\ sum(tsec=tpr,-s[1],(-1)^(tsec-tpr)*binomial(tsec,tpr)*xliste[1]*\ subst(taba[tsec+1,-s[1]+1],xlibre,xliste[1])) * \ la(N-1,nvsliste - tpr*decalagesuites,nvxliste)\ );\ );\ );\ resula /* ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- Maintenant : fichier dcpelsimples_revolution ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- */ /*Implemente jusqu'en prof 9 ; en particulier calculder est implemente en profondeur N quelconque certains details ont ete changes, e.g. 'varlibre remplac\'e par varlibre = 'varlibre, le 2/2/05 pour coller a la version de ce jour de dcpelsimples (qui est plus portable sur differentes machines, notamment). Et aussi prof 6 est passee a 9 (i.e. on a xvarsoul_9 maintenant). A ete teste le 8/2/05, et donne les memes resultats que le fichier dcpelsimples quant a la creation de deja_N3_ag222_sg222_plusclair Quelques tests suivent (test de la dcp en el. sples, test de valeurs de briques) (peut etre implementes slt si N <= 2) Dans calculder, parfois on renvoie une erreur si ca DV, alors que le calcul est tout a fait licite si x_i > 1... Mais le prog n'a pas ete suffisamment verifie (du pt de vue de la nv regularisation) dans ce cas. Attention, calculder est defini ici, mais utilise en fait la fct pochhammer qui est definie dans le fichier sorokin_andco.... */ xvarsoul_1 = 'xvarsoul_1; xvarsoul_2 = 'xvarsoul_2; xvarsoul_3 = 'xvarsoul_3; xvarsoul_4 = 'xvarsoul_4; xvarsoul_5 = 'xvarsoul_5; xvarsoul_6 = 'xvarsoul_6; xvarsoul_7 = 'xvarsoul_7; xvarsoul_8 = 'xvarsoul_8; xvarsoul_9 = 'xvarsoul_9; varsoul = [xvarsoul_1,xvarsoul_2,xvarsoul_3,xvarsoul_4,xvarsoul_5,xvarsoul_6,\ xvarsoul_7,xvarsoul_8,xvarsoul_9]; /*Les varsoul sont pour les fractions rationnelles ; ce sont les k_i sur lesquels on va sommer */ arrondir(fctquelconque) = fctquelconque; degpoly(polynomedp,vardp) =\ local(resuldp);\ if(polynomedp == 0, \ resuldp = - 10^10,\ resuldp = poldegree(polynomedp,vardp)\ );\ resuldp factorielle(nbfacto) = (nbfacto!); /*Ainsi definie elle renvoie un entier, alors que factorial(.) renvoie un nombre reel */ /*Attention les ssoul sont des indices dans les pochhammer, rien a voir avec La_{s} */ /*Dcp en el. sples, suivant bricksor.dvi du 9/11/03 (mais partie entiere y est oubliee) et manuscrit p. 31-32, de R(k_1,..,k_N) = P(k_1,..,k_N) / (prod_{i=1}^N (k_i + r_i)_{s_i +1} ^{a_i}) sous la forme \sum_{\sigmasoul = 0} ^\ssoul \sum{\alsoul = 0} ^\asoul c(qsp) / (prod_{i=1}^N (k_i + r_i + sigma_i) ^{\al_i}) La fonction coeffc renvoie le coeff c(qsp), qui depend de N, asoul, rsoul, ssoul, P, alsoul, sigmasoul Pour tenir compte de la partie entiere, on convient qu'on autorise sigma_i a valoir -1, et que dans ce cas a_i est remplace par max(0,deg_{k_i}(P)-(s_i+1)a_i) a_i - \al_i par max(..)+\al_i on a -max(...) \leq \al_i \leq 0 (alors que dans les autres cas on exige \al_i \geq 1) k_i par 1/k_i dans evaluation de R k_i + r_i + sigma_i par 1/k_i dans le produit le point -r_i-sigma_i ou on calcule c par zero */ /*Hyp : les soul sont des N-uplets, avec 0 \leq alphasoul \leq asoul et 0 \leq \sigmasoul \leq ssoul */ /*Test : P(xvarsoul_1) = (xvarsoul_1)^4; coeffc(1,P(xvarsoul_1),[0],[0],[1],[1],[0]); */ varlibre; /* Test : coeffc(2, Pnumr, [0,0],[1,1],[2,1],[2,1],[1,1]) doit donner 2 (page 34 manuscrit) qd Pnumr = (xvarsoul_2 - 1)*xvarsoul_1 + (-xvarsoul_2^2 + 2*xvarsoul_2 - 1) */ coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul) = \ local(polintermedcoeffc,remplaceai, nvpolyaux);\ prynt("On commence coeffc avec N = ", N, " et Pnumr = ",Pnumr, " et enfin rsoul,ssoul,asoul,alphasoul, sigmasoul = ",\ rsoul,ssoul,asoul,alphasoul, sigmasoul);\ polintermedcoeffc = Pnumr;\ varlibre = 'varlibre;\ if(N > wmax, prynt("ATTENTION !!!!! Depassement de wmax ; pb dcpelsimples_b"); 0/0;);\ for(icoeffc = 1, N,\ if((polintermedcoeffc != 0)&(degpoly(polintermedcoeffc,varsoul[icoeffc])<0),\ prynt("IL Y A UN PB dans dcpelsimples_a");\ 0/0;\ );\ polintermedcoeffc = polintermedcoeffc *\ (pochhammer(varsoul[icoeffc]+rsoul[icoeffc],\ ssoul[icoeffc]+1))^(-asoul[icoeffc]);\ prynt("Juste avant le distinguo, icoeffc = ",icoeffc," et polintermed = ",polintermedcoeffc);\ if(sigmasoul[icoeffc] >= 0,\ polintermedcoeffc = polintermedcoeffc *\ (varsoul[icoeffc]+rsoul[icoeffc]+\ sigmasoul[icoeffc])^(asoul[icoeffc]);\ polintermedcoeffc = (1/factorielle(asoul[icoeffc]-alphasoul[icoeffc]))*\ polintermedcoeffc;\ prynt("Avant derivation, on a polintermed = ",polintermedcoeffc);\ for(jcoeffc = 1, asoul[icoeffc]-alphasoul[icoeffc],\ polintermedcoeffc = deriv(arrondir(polintermedcoeffc),varsoul[icoeffc]);\ );\ prynt("Apres derivation, on a polintermed = ",polintermedcoeffc);\ nvpolyaux = subst(polintermedcoeffc, varsoul[icoeffc], varlibre);\ prynt("Apres substitution on a nvpolyaux = ", nvpolyaux);\ varlibre = -rsoul[icoeffc] - sigmasoul[icoeffc];\ prynt("Avant eval en ", -rsoul[icoeffc] - sigmasoul[icoeffc]," on a polintermed = ",polintermedcoeffc);\ polintermedcoeffc = eval(nvpolyaux);\ prynt("Apres eval on a polintermed = ",polintermedcoeffc);\ varlibre = 'varlibre;\ );\ if(sigmasoul[icoeffc] == -1,\ if(polintermedcoeffc != 0,\ remplaceai = max(0,\ degpoly(Pnumr,varsoul[icoeffc])-asoul[icoeffc]*(1+ssoul[icoeffc]));\ prynt("On a remplaceai = ", remplaceai);\ polintermedcoeffc = subst(polintermedcoeffc,varsoul[icoeffc],1/varlibre);\ polintermedcoeffc = polintermedcoeffc * varlibre^(remplaceai);\ polintermedcoeffc = (1/factorielle(remplaceai+alphasoul[icoeffc]))*\ polintermedcoeffc;\ for(jcoeffc = 1, remplaceai + alphasoul[icoeffc],\ polintermedcoeffc = deriv(polintermedcoeffc,varlibre);\ );\ varlibre = 0;\ prynt("Avant eval en varlibre = 0, on a polintermed = ",polintermedcoeffc);\ polintermedcoeffc = eval(polintermedcoeffc);\ prynt("Apres eval on a polintermed = ",polintermedcoeffc);\ varlibre = 'varlibre;\ );\ );\ );\ polintermedcoeffc /*Polynome sommebernoulli(s,xbern) = P_s(k,x) de TR, p. 13 ; avec s >= 0. C'est vu comme un polynome en la variable varbern, qui correspond a k (et a laquelle on pourrait donner des valeurs >= 0), a coefficients qui dependent de xlibre (cf. l'implementation de taba, dans le fichier polylog)(mais on remplace xlibre par xbern dans le calcul de sommebernoulli). Par def P_s(k,x) = sum_{l=1} ^k l^s x^(-l) On a P_s(k,x) = x^(-k) A_{x,s}(k) - a_{0,x,s} pour x \neq 0,1 o\`u A_{x,s}(k) = sum_{t=0} ^s a_{t,x,s} k^t est un polynome en k de degre s, dont les coeffs sont des fractions rationnelles en x avec 1 pour seul pole. Si x=1 l'evaluation de P_s(k,1) se fait autrement, par formule avec polynomes de Bernoulli ; on pose sommebernoullifacile(s) = P_s(k,1) qui est un polynome en varbern=k si s=0 c'est varbern (i.e. la var qui joue le role de k dans sommebernoulli, y compris dans sa version facile) si s=1 : T*(T+1)/2 avec T= varbern si s=2 : T(T+1)(2T+1)/6, etc. En tout cas c'est un polynome de degre s+1, sans terme constant. Si x!=1, sommebernoulli(s,xbern) renvoie un tableau a 2 cases (si x=1 c'est le cas aussi, mais la 2eme est vide). La premiere est le coeff de x^(-k), la deuxieme le coeff de 1. Les deux sont des polynomes en varbern, a coefficients qui dependent de xbern. */ varbern = 'varbern; sommebernoullifacile(s)=\ local(nvresubernfac);\ if(s == 0, nvresubernfac = varbern,\ nvresubernfac = (1/(s+1))*(polynomebernoulli(s+1,varbern+1)-bernfrac(s+1)));\ nvresubernfac kroneckerdelta(akro,brko) = \ local(resulkro);\ if(akro == bkro, resulkro = 1, resulkro = 0);\ resulkro /*Polynome de Bernoulli http://mathworld.wolfram.com/BernoulliPolynomial.html (confirme par Tanguy) B_n(x) = (B+x)^n ou on interprete B^k comme B_k */ polynomebernoulli(nbern,xintermedbern) = \ sum(kbern=0,nbern,bernfrac(kbern)*binomial(nbern,kbern)*\ xintermedbern^(nbern-kbern)) sommebernoulli(s,xbern) = \ local(resubern);\ if(xbern == 1,\ resubern = [sommebernoullifacile(s),0];\ ,\ resubern = \ [sum(tbern=0,s, subst(taba[tbern+1,s+1],xlibre,xbern)*varbern^tbern), \ -subst(taba[1,s+1],xlibre,xbern)];\ );\ resubern /*cf. manuscrit page 71-75 */ /*Implemente slt si imaxdes = N, car c'est le seul cas utile pour calculder*/ /*Si N=2 : renvoie un tableau a 2 cases ; chacune est une frac rat. en varsoul[1],... varsoul[N-1], dont les coeffs dependent de xliste[N] */ /*imaxdes,ssouldes,famjdes correspondent a la brique qu'on est en train de considerer; rsoulder,ssoulder,asoulder correspondent au calculder qu'on est en train de faire; N et xliste sont communs aux deux. */ /* 2 2 [ 0, 0] [2, 2] [2, 2] [1, 1] [0, 0] [ 1, 1] */ contributiondunDV(N,imaxdes,rsoulder,ssoulder,asoulder,ssouldes,famjdes,xliste)=\ local(resucontriDV,resucontriDVtemp);\ if(N != imaxdes, prynt("ContributiondunDV implementee slt qd imaxdes=N ; pb dcpelsimples_c");0/0);\ if(N <= 1, prynt("ContributiondunDV non implementee pour N = ",N," ; pb dcpelsimples_d");0/0);\ prynt("On commence ContributiondunDV avec ",N,imaxdes,rsoulder,ssoulder,asoulder,ssouldes,famjdes,xliste);\ if(N == 2,\ resucontriDV = \ subst(sommebernoulli(-ssouldes[N],xliste[N]),varbern,varsoul[N-1])/\ prod(icDV=1,N-1,(varsoul[icDV]+famjdes[icDV])^ssouldes[icDV]);\ );\ prynt("Milieu de ContributiondunDV avec ",N,imaxdes,rsoulder,ssoulder,asoulder,ssouldes,famjdes,xliste);\ if((N >= 3) & (ssouldes[N] <= 0),\ resucontriDVtemp = \ subst(sommebernoulli(-ssouldes[N],xliste[N]),varbern,varsoul[N-1])/\ prod(icDV=1,N-1,(varsoul[icDV]+famjdes[icDV])^ssouldes[icDV]);\ resucontriDV = vector(N);\ resucontriDV[N-1] = resucontriDVtemp[1];\ resucontriDV[N] = resucontriDVtemp[2];\ );\ if((N >= 3) & (ssouldes[N] >= 1),\ indnegDV = N;\ resucontriDV = vector(N);\ while((indnegDV >= 2) && (ssouldes[indnegDV] >= 1),indnegDV=indnegDV-1);\ if(indnegDV <= 1, prynt("Il y a un pb de DV.... dans dcpelsimples_e"); 0/0);\ resucontriDV[indnegDV-1] = \ subst(sommebernoulli(-ssouldes[indnegDV],xliste[indnegDV])[1],\ varbern,varsoul[indnegDV-1]);\ resucontriDV[indnegDV] = -xliste[indnegDV] * \ subst(sommebernoulli(-ssouldes[indnegDV],xliste[indnegDV])[1],\ varbern,varsoul[indnegDV]-1);\ for(kindDV = indnegDV-1, indnegDV, \ resucontriDV[kindDV] = resucontriDV[kindDV] /\ prod(icDV=1,indnegDV-1,\ (varsoul[icDV]+famjdes[icDV])^ssouldes[icDV])/\ prod(icDV=indnegDV+1,N,\ (varsoul[icDV-1]+famjdes[icDV])^ssouldes[icDV]);\ );\ );\ prynt("On finit ContributiondunDV avec ",N,imaxdes,rsoulder,ssoulder,asoulder,ssouldes,famjdes,xliste);\ resucontriDV /*Cf. manuscrit p. 76 Le probleme est que dans contributiondunDV, l'un des k_j disparait, donc on se retrouve avec des nouveaux k_j, pour 2 <= j <= N-1, et au denom on a (k_j + rsoul[j])_{ssoul[j]+1} ^asoul[j] ou bien (k_j + rsoul[j+1])_{ssoul[j+1]+1} ^asoul[j+1] Cette fct renvoie rsoulfu,ssoulfu,asoulfu tq (k_j + rsoulfu[j])_{ssoulfu[j]+1} ^asoulfu[j] soit un multiple commun de ces deux denom possibles. Pour j = 1, pas besoin : poss prendre rsoulfu[1] = rsoul[1] etc. */ calculdesfuturescoordr(N,rsoul,ssoul,asoul)=\ local(rsoulfu,ssoulfu,asoulfu);\ rsoulfu = vector(N-1);\ ssoulfu = vector(N-1);\ asoulfu = vector(N-1);\ rsoulfu[1] = rsoul[1];\ ssoulfu[1] = ssoul[1];\ asoulfu[1] = asoul[1];\ for(jfu = 2, N-1, \ rsoulfu[jfu] = min(rsoul[jfu],rsoul[jfu+1]);\ ssoulfu[jfu] = max(rsoul[jfu]+ssoul[jfu],rsoul[jfu+1]+ssoul[jfu+1])\ - rsoulfu[jfu];\ if((rsoul[jfu]+ssoul[jfu] < rsoul[jfu+1]) || (rsoul[jfu+1]+ssoul[jfu+1] \ < rsoul[jfu]),\ asoulfu[jfu] = max(asoul[jfu],asoul[jfu+1])\ ,\ asoulfu[jfu] = asoul[jfu] + asoul[jfu+1]\ );\ );\ [rsoulfu,ssoulfu,asoulfu] /*Pnumr est un polynome en les N variables de varsoul*/ /*Implemente pour tout N (mais peu teste); voir l'autre fichier dcpelsimples pour l'implementation habituelle quand N <= 3 */ /*Avec cette la nv reg, il faut l'appliquer en xliste = [1,..1] si on utilise une regularisation (en relisant cette phrase, me semble bien mysterieux... de toute facon, a peu pres rien n'a ete teste serieusement quand les variables ne valent pas 1) */ /*On pourrait peut-etre acceler un peu (a peine...) le calculder dans le cas CV si on etait sur que sigmasoul[1] = -1 est impossible en cas de CV. */ /*Calcul de \sum_{k_1 >= ... >= k_N >= 1} x_1^{-k_1}...x_N^{-k_N} R(ksoul) avec R(k_1,..,k_N) = P(k_1,..,k_N) / (prod_{i=1}^N (k_i + r_i)_{s_i +1} ^{a_i}) */ /*On suppose A PRIORI qu'on peut se restreindre aux sigma tels que sigmasoul[1] >= 0, sinon on est trivialt dans un cas grossierement DV. Ceci est fait juste avant la premiere occurrence de : if(touspositifs, \ sommedesCV = sommedesCV + coer * \ briqbn(N,N,alphasoul,rsoul + sigmasoulpourb,xliste);\ et permet d'accelerer le calcul, en sachant a priori que des coeffs du dvlpt en el. sples. seront nuls par CV. Mais ceci suppose qu'on est au point 1,1,..,1, et aussi qu'on traite une serie CV. Aucun effort n'a ete fait pour supprimer ces hyp. */ /*Le while((icoura > 0) && ... a ete bidouille le 2/2/05 pour essayer de faire en sorte que ca tourne...*/ /*pour un test : calculder(3, 1, [0,0,1], [2,2,1], [2,2,2], [1,1,1]) */ calculder(N,Pnumr,rsoul,ssoul,asoul,xliste) =\ local(alphasoul,sigmasoul,sigmasoulpourb,sommedesCV,sommedesDV,\ partieDV,numerateurdepartieDV,coer,nvxliste,resulcalculder);\ local(futurescoordr,rsoulfutur,ssoulfutur,asoulfutur,onestenun);\ prynt("On utilise calculder avec N = ",N," et xliste = ",xliste);\ capadcp = vector(N);\ for(i = 1, N, \ capadcp[i] = 1 + max(asoul[i]*(1+ssoul[i]), degpoly(Pnumr,varsoul[i]));\ );\ resulcalculder = 0;\ sommedesCV = 0;\ partieDV = vector(N);\ sigmasoul = vector(N);\ sigmasoulpourb = vector(length(rsoul));\ alphasoul = vector(N);\ for(i=1, N, sigmasoul[i] = -1);\ for(i=1, N, alphasoul[i] = -(capadcp[i] - 1 - asoul[i]*(1+ssoul[i])));\ \ prynt("Etape 1a de calculder avec", sigmasoul, alphasoul,capadcp);\ for(inddcp = 1, prod(i=1,N,capadcp[i]),\ touspositifs = 1;\ for(i=1, N, \ if(sigmasoul[i] == -1, sigmasoulpourb[i] = -rsoul[i]; touspositifs = 0; ,\ sigmasoulpourb[i] = sigmasoul[i]\ );\ );\ prynt("Etape 1b de calculder avec", sigmasoul, alphasoul,capadcp,inddcp);\ if(sigmasoul[1] >= 0,\ coer = coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul,sigmasoul);\ prynt("Etape 1c de calculder avec N = ",N," et xliste = ",xliste,\ " et touspositifs = ",touspositifs,\ " et ",rsoul,ssoul,asoul,alphasoul,rsoul + sigmasoulpourb,xliste);\ if(coer != 0,\ if(touspositifs, \ sommedesCV = sommedesCV + coer * \ briqbn(N,N,alphasoul,rsoul + sigmasoulpourb,xliste);\ ,\ partieDV = partieDV + coer*\ contributiondunDV(N,N,rsoul,ssoul,asoul,alphasoul,rsoul + sigmasoulpourb,xliste);\ );\ );\ );\ prynt("On est au milieu des pre-calculs de calculder avec N,inddcp,capadcp = ",N,inddcp,capadcp);\ tempinddcp = inddcp;\ icoura = N;\ while((icoura > 0) && (Mod(tempinddcp,capadcp[icoura]) == Mod(0,capadcp[icoura])),\ prynt("Etape 1d de calculder avec N = ",N," et xliste = ",xliste);\ tempinddcp = tempinddcp / capadcp[icoura] ;\ icoura = icoura - 1;\ );\ for(i = icoura + 1, N, \ sigmasoul[i] = -1;\ alphasoul[i] = -(capadcp[i] - 1 - asoul[i]*(1+ssoul[i]));\ );\ prynt("Etape 1e de calculder avec N = ",N," et xliste = ",xliste);\ if(icoura > 0,\ prynt("Etape 1f de calculder avec N = ",N," et xliste = ",xliste);\ if(sigmasoul[icoura] == -1,\ if(alphasoul[icoura] < 0,\ alphasoul[icoura] = alphasoul[icoura] + 1;\ ,\ if(alphasoul[icoura] == 0,\ sigmasoul[icoura] = 0;\ alphasoul[icoura] = 1;\ );\ );\ ,\ if(alphasoul[icoura] < asoul[icoura],\ alphasoul[icoura] = alphasoul[icoura] + 1;\ ,\ if(alphasoul[icoura] == asoul[icoura],\ sigmasoul[icoura] = sigmasoul[icoura] + 1;\ alphasoul[icoura] = 1;\ prynt("Etape 1j de calculder avec", sigmasoul, ssoul, icoura);\ if(sigmasoul[icoura] > ssoul[icoura], \ print("Probleme de valuation dcpelsimples");0/0;\ );\ );\ );\ );\ );\ prynt("Etape 1z de calculder avec N = ",N," et xliste = ",xliste);\ );\ \ \ if((N == 1) & (xliste[1] == 1), sommedesDV = 0);\ if((N == 1) & (xliste[1] != 1), \ prynt("On n'a pas implemente le cas DV, meme quand x_i \neq 1");\ prynt("en plus ici N = 1... donc peut etre ca ne DV jamais ? pb dcpelsimples_j");0/0;\ prynt("En tout cas il faut definir numerateurdepartieDV maintenant");\ /*Faire ceci est probablement nec pour que ca tourne qd x_i != 1*/ \ sommedesDV = calculder(N-1,numerateurdepartieDV[2],\ rsoul,ssoul,asoul,xliste);\ nvxliste = xliste;\ nvxliste[N-1] = xliste[N-1] * xliste[N];\ sommedesDV = sommedesDV + calculder(N-1,numerateurdepartieDV[1],\ rsoul,ssoul,asoul,nvxliste);\ );\ \ \ \ if(N >= 2,\ futurescoordr = calculdesfuturescoordr(N,rsoul,ssoul,asoul);\ rsoulfutur = futurescoordr[1];\ ssoulfutur = futurescoordr[2];\ asoulfutur = futurescoordr[3];\ numerateurdepartieDV = partieDV * prod(inumDV = 1, N-1, \ (pochhammer(varsoul[inumDV] + rsoulfutur[inumDV], \ 1 + ssoulfutur[inumDV]))^(asoulfutur[inumDV]));\ onestenun = 1;\ for(ioeeu = 2, N, \ if(xliste[ioeeu] != 1, onestenun = 0);\ );\ if(onestenun,\ sommedesDV = calculder(N-1,\ sum(inu = 1, N, numerateurdepartieDV[inu]),\ rsoulfutur,ssoulfutur,asoulfutur,xliste);\ ,\ prynt("L'implementation du cas DV, meme quand x_i \neq 1, n'a pas ete testee suffisamment");\ for(iDV = 1, N, \ nvxliste = xliste;\ if(iDV <= N-1,\ nvxliste[iDV] = xliste[iDV] * xliste[iDV+1];\ for(jDV = iDV + 1, N-1, nvxliste[jDV] = xliste[jDV + 1]);\ );\ sommedesDV = calculder(N-1,numerateurdepartieDV[iDV],\ rsoul,ssoul,asoul,nvxliste);\ );\ );\ );\ resulcalculder = sommedesCV + sommedesDV;\ resulcalculder /* -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------- */ /* VRAIES VALEURS CALCULEES PAR SOMMES DE SERIES ------------------------------------------------- */ /*Implemente slt quand N= 2 */ vraievaleurder(N,Pnumr,rsoul,ssoul,asoul,xlistevrair) =\ local(denomr,resulvrair);\ resulvrair = 0;\ quotientr = Pnumr / prod(ivrair = 1, N, \ (pochhammer(varsoul[ivrair] + rsoul[ivrair], 1 + ssoul[ivrair])) \ ^(asoul[ivrair]));\ if(N == 2, \ resulvrair = suminf(kvrairun = 1, \ xlistevrair[1]^(-kvrairun) *\ sum(kvrairde = 1, kvrairun, \ xlistevrair[2]^(-kvrairde) *\ subst(subst(quotientr,xvarsoul_1,kvrairun),xvarsoul_2,kvrairde)\ ));\ );\ if(N == 3, \ resulvrair = suminf(kvrairun = 1, \ xlistevrair[1]^(-kvrairun) *\ sum(kvrairde = 1, kvrairun, \ xlistevrair[2]^(-kvrairde) *\ sum(kvrairtr = 1, kvrairde, \ xlistevrair[3]^(-kvrairtr) *\ subst(subst(subst(quotientr,xvarsoul_1,kvrairun),xvarsoul_2,kvrairde),\ xvarsoul_3,kvrairtr)\ )));\ );\ resulvrair /* PARTIE QUI SERT SLT AUX TESTS ------------------------------------------------------------------ */ testcalculder(N,Pnumr,rsoul,ssoul,asoul,xliste) =\ local(alphasoul,sigmasoul,resulcalculder);\ prynt("On commence calculder");\ prynt("avec Pnumr = ",Pnumr,", asoul = ",asoul);\ prynt("rsoul =",rsoul,"ssoul = ",ssoul);\ sigmasoul = vector(N);\ alphasoul = vector(N);\ if(N != 2, prynt("PROBLEME : IMPLEMENTE SLT QUAND N=2"));\ if(N == 2,\ sigmasoul[1] = -1;\ prynt("Etape 1 de calculder");\ for(alpha_1 = \ -max(0,degpoly(Pnumr,varsoul[1])-asoul[1]*(1+ssoul[1])) , 0 ,\ alphasoul[1] = alpha_1;\ sigmasoul[2] = -1;\ prynt("Etape 2 de calculder");\ for(alpha_2 =\ -max(0,degpoly(Pnumr,varsoul[2])-asoul[2]*(1+ssoul[2])) , 0 ,\ alphasoul[2] = alpha_2;\ prynt("Etape 3 de calculder, avec Pnumr = ",Pnumr," et alphasoul = ",alphasoul);\ prynt("et sigmasoul = ",sigmasoul);\ nbtest = nbtest + 1;\ ldtest[nbtest] = [asoul,ssoul,alphasoul,sigmasoul];\ tabresulttest[nbtest] = coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul);\ prynt("coeffc donne ",coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul));\ prynt("Avec nbtest = ",nbtest);\ prynt("Etape 3bis de calculder");\ );\ for(sigma_2 = 0, ssoul[2],\ sigmasoul[2] = sigma_2;\ for(alpha_2 = 1, asoul[2],\ alphasoul[2] = alpha_2;\ prynt("Etape 4 de calculder");\ nbtest = nbtest + 1;\ ldtest[nbtest] = [asoul,ssoul,alphasoul,sigmasoul];\ tabresulttest[nbtest] = coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul);\ prynt("coeffc donne ",coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul));\ );\ );\ );\ for(sigma_1 = 0, ssoul[1],\ sigmasoul[1] = sigma_1;\ for(alpha_1 = 1, asoul[1],\ alphasoul[1] = alpha_1;\ sigmasoul[2] = -1;\ for(alpha_2 =\ -max(0,degpoly(Pnumr,varsoul[2])-asoul[2]*(1+ssoul[2])) , 0 ,\ alphasoul[2] = alpha_2;\ prynt("etape 5a, avec alphasoul = ",alphasoul);\ nbtest = nbtest + 1;\ ldtest[nbtest] = [asoul,ssoul,alphasoul,sigmasoul];\ tabresulttest[nbtest] = coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul);\ prynt("coeffc donne ",coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul));\ );\ for(sigma_2 = 0, ssoul[2],\ prynt("etape 5b");\ sigmasoul[2] = sigma_2;\ for(alpha_2 = 1, asoul[2],\ prynt("Etape 6 de calculder");\ alphasoul[2] = alpha_2;\ nbtest = nbtest + 1;\ ldtest[nbtest] = [asoul,ssoul,alphasoul,sigmasoul];\ tabresulttest[nbtest] = coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul);\ prynt("coeffc donne ",coeffc(N,Pnumr,rsoul,ssoul,asoul,alphasoul, sigmasoul));\ );\ );\ );\ );\ );\ prynt("Etape 7 de calculder");\ /*Ce programme verifie que la dcp en elements simples qu'on calcule pour l'integrale de Sorokin est correcte. Et ca semble etre le cas, au moins pour exposantn = 1*/ verifelementssimples_sorokinzetatrois(exposantn) =\ local(paramlocal,polynomenumsoro,verifsoro,autreverif);\ polynomenumsoro = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ prynt("Le numerateur de Sorokin est ",polynomenumsoro);\ nbtest = 0;\ ldtest = vector(100);\ tabresulttest = vector(100);\ testcalculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,1],ldvar);\ verifsoro = 0;\ ssoul = [exposantn,exposantn];\ asoul = [2,1];\ for(iveri = 1, nbtest,\ sigmasoul = ldtest[iveri][4];\ alphasoul = ldtest[iveri][3];\ paramlocal = 1;\ if(sigmasoul[1] == -1, /*ici rsoul = 0 */ \ paramlocal = paramlocal * varsoul[1]^(-alphasoul[1]),\ paramlocal = paramlocal * (varsoul[1]+sigmasoul[1])^(-alphasoul[1])\ );\ if(sigmasoul[2] == -1, /*ici rsoul = 0 */ \ paramlocal = paramlocal * varsoul[2]^(-alphasoul[2]),\ paramlocal = paramlocal * (varsoul[2]+sigmasoul[2])^(-alphasoul[2])\ );\ verifsoro = verifsoro + tabresulttest[iveri]*paramlocal;\ );\ autreverif = polynomenumsoro / prod(ivrair = 1, 2, \ (pochhammer(varsoul[ivrair], 1 + ssoul[ivrair])) \ ^(asoul[ivrair]));\ [verifsoro,autreverif] /*Pr exposantn=1, donne deux fois la meme chose, donc la dcp en elements simples est correcte! */ /*Ce programme verifie la valeur de chacune des briques internes utilisees pour calculer sorokin */ verifdesbriquesinternes_sorokinzetatrois(exposantn) =\ local(paramlocal,polynomenumsoro,verifsoro,autreverif);\ polynomenumsoro = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ prynt("Le numerateur de Sorokin est ",polynomenumsoro);\ nbtest = 0;\ ldtest = vector(100);\ tabresulttest = vector(100);\ testcalculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,1],ldvar);\ ssoul = [exposantn,exposantn];\ asoul = [2,1];\ tabvb = vector(nbtest);\ for(iveri = 1, nbtest,\ sigmasoul = ldtest[iveri][4];\ alphasoul = ldtest[iveri][3];\ tabvb[iveri] = vdeux_pourverifier_briqbn(2,2,alphasoul,sigmasoul);\ );\ for(iveri = 1, nbtest,\ if(tabvb[iveri] > 10^(-10),\ prynt("IL Y A UN PB pour iveri =", iveri);\ );\ );\ /* POur test : for(i=1, 100, \ if((ldtest[i][4][1] == -1) & (ldtest[i][4][1] == -1), prynt(i));\ ); ldtest[nbtest] = [asoul,ssoul,alphasoul,sigmasoul];\ */ /* ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- Maintenant : fichier briques (qui, initialement, appelait les 3 precedents) ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- */ longueuractu = 0; /*Hyp : rsoulconnu[i] >= 0 et tous les x_i = 1 (calcul geometrique non implemente sinon)*/ /*les variables globales longueuractu, ldactuDg, ldactudimedg, ldactubaseedg, ldactumatzetag ldactuN, ldactursoul, ldactuasoul, ldactussoul doivent avoir une valeur; Si longueuractu <= 0 (ou si les longueuractu premiers elements de actuN sont <= 1), alors on n'utilisera jamais le calcul geometrique donc il n'est pas imperatif de supposer rsoulconnu[i] >= 0 et tous les x_i = 1 */ /*Pour utiliser le precalcul geom, il faut avoir appris ce qu'est la fonction p2c qui se trouve dans le fichier initgeom*/ /*OK en profondeur quelconque*/ /* Attention, calculder est defini ici, mais utilise en fait la fct pochhammer qui est definie dans le fichier sorokin_andco.... */ /* \r initialisation_poids19plusclair \r polylogs \r dcpelsimples_revolution */ /*_revolution*/ /*Ce fichier est valable quel que soit wmax. wmax intervient seulement dans l'initialisation ci-dessus */ /*remplirtabvraisla; C'est un precalcul utile pour valeurs num des polylogs ; il prend un peu de temps en poids 3 (quelques secondes, en \p 9) cf fichier initialisation_poids3_polylogs (c'est l'endroit ou il sert, mais slt pour verif). Quand on initialise avec initialisation_poids3_numerique, ce fichier est execute automatiquement. En fait, init_polylogs enchaine avec une vraievaleurduntab donne theoriquement la meme chose que si on initialise avec init_numerique. Les problemes sont pratiques : init_polylogs prend plus de memoire (utilise des tableaux a 38 cases au lieu de reels), init_numerique pose des problemes du type *** Warning: normalizing a polynomial with 0 leading term. cf. au debut de init_numerique pour des exemples numeriques (et pour une mise en garde sur la precision obtenue...) */ /*partout on autorise s a depasser au dela de N cases, mais les suivantes sont ignorees ; idem pour les var x*/ /* pour tester : briqbn(3,1,[1,-2,-4],[1,4,5],ldvar) fait exploser la pile initiale avec initialisation polylogs ! */ /*cf. page 109-110 ; hyp pas de modulation ; de plus pour N < Nconnu il semble peu judicieux d'utiliser un calcul geometrique*/ /* Attention, hyp les rsoulconnu[i] sont >= 0 */ facileacalculerbn(N,imax,s,famj,xliste,Nconnu, asoulconnu, rsoulconnu, ssoulconnu)=\ local(resufaccalcbn);\ resufaccalcbn = (N > 1) && (N == Nconnu) && (imax >= N) ;\ if(resufaccalcbn,\ for(icbn = 1, N, resufaccalcbn = resufaccalcbn && (s[icbn] <= asoulconnu[icbn]) && \ (rsoulconnu[icbn] <= famj[icbn]) && (famj[icbn] <= rsoulconnu[icbn] + ssoulconnu[icbn]);\ );\ );\ resufaccalcbn /*Ceci est l'ancienne version, qui admet N < Nconnu*/ /* facileacalculerbn(N,imax,s,famj,xliste,Nconnu, asoulconnu, rsoulconnu, ssoulconnu)=\ local(resufaccalcbn);\ resufaccalcbn = (N > 1) && (N <= Nconnu) && (imax >= N) ;\ if(resufaccalcbn,\ for(icbn = 1, N, resufaccalcbn = resufaccalcbn && (s[icbn] <= asoulconnu[icbn]) && \ (rsoulconnu[icbn] <= famj[icbn]) && (famj[icbn] <= rsoulconnu[icbn] + ssoulconnu[icbn]);\ );\ for(icbn = N, Nconnu - 1, resufaccalcbn = resufaccalcbn && \ (rsoulconnu[icbn] == 0) && (asoulconnu[icbn] >= 1);\ );\ if((famj[N] == 0) && (N < Nconnu), resufaccalcbn = resufaccalcbn && (s[N] + 1 <= asoulconnu[N]));\ );\ resufaccalcbn */ /*Dans briqbnconnu on suppose que facileacalculerbn = true */ /*On utilise la fonction p2c*/ /*Attention, le calcul geometrique n'est implemente que quand tous les x_i sont egaux a 1, donc cette procedure ne l'est que dans ce cas*/ calculbriqbnconnu(N,imax,s,famj,xliste,Nconnu, asoulconnu, rsoulconnu, ssoulconnu, Dsoulconnu, \ donnedimedconnu, donnebaseedconnu, matzetaconnu) = \ local(resubriqbnconnu, polybnconnu, vecteurligneconnu);\ polybnconnu = prod(icbn = 1, N, pochhammer(varsoul[icbn]+rsoulconnu[icbn], \ ssoulconnu[icbn]+1)^(asoulconnu[icbn]) / ((varsoul[icbn] + famj[icbn])^(s[icbn])));\ vecteurligneconnu = p2c(polybnconnu,Nconnu, Dsoulconnu,donnedimedconnu, donnebaseedconnu);\ resubriqbnconnu = mattranspose(matzetaconnu*(vecteurligneconnu~));\ resubriqbnconnu /*Ceci est l'ancienne version, qui admet N < Nconnu*/ /* calculbriqbnconnu(N,imax,s,famj,xliste,Nconnu, asoulconnu, rsoulconnu, ssoulconnu, Dsoulconnu, \ donnedimedconnu, donnebaseedconnu, matzetaconnu) = \ local(resubriqbnconnu, polybnconnu, vecteurligneconnu);\ polybnconnu = prod(icbn = 1, N, pochhammer(varsoul[icbn]+rsoulconnu[icbn], \ ssoulconnu[icbn]+1)^(asoulconnu[icbn]) / ((varsoul[icbn] + famj[icbn])^(s[icbn])));\ polybnconnu = polybnconnu *\ prod(icbn = N+1, Nconnu, pochhammer(varsoul[icbn]+rsoulconnu[icbn], \ ssoulconnu[icbn]+1)^(asoulconnu[icbn]) / (varsoul[icbn-1]));\ vecteurligneconnu = p2c(polybnconnu,Nconnu, Dsoulconnu,donnedimedconnu, donnebaseedconnu);\ resubriqbnconnu = mattranspose(matzetaconnu*(vecteurligneconnu~));\ resubriqbnconnu */ /* voir p. 24 manuscrit : c'est B_N(s, famj, moduln) ou famj[t] = j_t de TR pour t <= imax, j_t de TR = 0 pour t > imax, et m_t de TR = 0 sauf pour t = imax+1 <= N, et dans ce cas vaut famj[imax+1] Dans les arguments de briquebn on suppose N >= 1 et imax >=0, avec les s[i] dans Z et les famj[t] dans N On neglige les termes de s apres le Neme, et ceux de famj apres le (imax+1)eme Pas de modulation ssi imax >= N */ /*probablement OPS imax >= 1 car m_1 est toujours nul */ /*calcul par Prop 2 de TR, p. 4, si on n'utilise pas le pre-calcul geometrique */ /*tests possibles : briqbn(3,1,[1,2,4],[1,4,5], [1,1,1]) briqbn(3,1,[3,2,2],[2,6,5], [1,1,1]) */ briqbn(N,imax,s,famj,xliste) = \ local(resulbriqbn,nvimax,facileacalc, iicalbn);\ facileacalc = 0;\ iicalbn = 1;\ while((iicalbn <= longueuractu) && (facileacalc == 0),\ if(facileacalculerbn(N,imax,s,famj,xliste,\ ldactuN[iicalbn], ldactuasoul[iicalbn], ldactursoul[iicalbn], ldactussoul[iicalbn]),\ facileacalc = 1;\ resulbriqbn = calculbriqbnconnu(N,imax,s,famj,xliste,\ ldactuN[iicalbn], ldactuasoul[iicalbn], ldactursoul[iicalbn], ldactussoul[iicalbn], \ ldactuDg[iicalbn], ldactudimedg[iicalbn], ldactubaseedg[iicalbn], ldactumatzetag[iicalbn]);\ ,\ iicalbn = iicalbn +1;\ );\ );\ if(facileacalc == 0,\ if(N == 0, resulbriqbn = valeurconstante(1),\ nvimax = min(imax,N);\ prwnt("Etape 1 de briqbn avec N,imax,s,famj,xliste = ",N,imax,s,famj,xliste);\ resulbriqbn = prod(tintermu=1,nvimax,xliste[tintermu]^(famj[tintermu])) *\ la(N,s,xliste);\ prwnt("Etape 2 de briqbn");\ resulbriqbn = resulbriqbn - sum(pintermu=1, nvimax, \ prod(tintermd=pintermu,nvimax,xliste[tintermd]^(famj[tintermd])) * \ briqqn(N,pintermu,s,famj[pintermu],xliste) *\ briqbn(pintermu-1,imax,s,famj,xliste)\ );\ prwnt("Etape 3 de briqbn");\ for(ppp=2, min(N,imax+1), \ resulbriqbn = resulbriqbn + \ sign(famj[ppp]-famj[ppp-1]) *\ prod(tintermd=ppp,nvimax,xliste[tintermd]^(famj[tintermd])) * \ sum(kpbriqbn=1+min(famj[ppp-1],famj[ppp]), max(famj[ppp-1],famj[ppp]),\ xliste[ppp]^(-kpbriqbn) *\ briqrn(N,ppp,s,famj,kpbriqbn,xliste)\ );\ );\ );\ );\ resulbriqbn /*memes remarques que pour briqbn, sauf qu'ici OPS que la moduln est tjs nulle. Donc famj[t] = j_t de TR pour t <= ppprn-1, et les cases suivantes de famj sont ignor\'ees. Hyp : 2 <= ppprn <= N et K >= 1 Attention, briqrn est une fonction de x[1],...,x[N] meme si elle est consideree dans TR comme fct de _p x_N */ briqrn(N,ppprn,s,famj, K, xliste) = \ local(resulbriqrn, nvxlisteppp,nvsrn,nvfamj);\ nvxlisteppp = xliste;\ nvxlisteppp[ppprn-1] = xliste[ppprn-1] * xliste[ppprn];\ for(pppinterm = ppprn, N-1, nvxlisteppp[pppinterm] = xliste[pppinterm+1]);\ if(K == famj[ppprn-1],\ nvsrn = s ;\ nvsrn[ppprn-1] = s[ppprn-1] + s[ppprn];\ for(pppin = ppprn, N-1, nvsrn[pppin] = s[pppin+1]);\ nvfamj = famj ;\ if(ppprn <= N-1, nvfamj[ppprn] = famj[ppprn-1]);\ resulbriqrn = briqbn(N-1,ppprn-1,nvsrn,nvfamj,nvxlisteppp);\ ,if((s[ppprn-1] <= 0) & (s[ppprn] >= 1),\ resulbriqrn = valeurconstante(0);\ for(urn = 0, -s[ppprn-1],\ nvsrn = s ;\ nvsrn[ppprn-1] = s[ppprn] - urn;\ for(pppin = ppprn, N-1, nvsrn[pppin] = s[pppin+1]);\ nvfamj = famj ;\ nvfamj[ppprn-1] = K;\ if(ppprn <= N-1, nvfamj[ppprn] = K);\ resulbriqrn = resulbriqrn + binomial(-s[ppprn-1],urn) * \ (famj[ppprn-1]-K)^(-s[ppprn-1]-urn) *\ briqbn(N-1,ppprn-1,nvsrn,nvfamj,nvxlisteppp);\ );\ ,if((s[ppprn-1] >= 1) & (s[ppprn] <= 0),\ resulbriqrn = valeurconstante(0);\ for(urn = 0, -s[ppprn],\ nvsrn = s ;\ nvsrn[ppprn-1] = s[ppprn-1] - urn;\ for(pppin = ppprn, N-1, nvsrn[pppin] = s[pppin+1]);\ nvfamj = famj ;\ if(ppprn <= N-1, nvfamj[ppprn] = K);\ resulbriqrn = resulbriqrn + binomial(-s[ppprn],urn) * \ (K - famj[ppprn-1])^(-s[ppprn]-urn) *\ briqbn(N-1,ppprn-1,nvsrn,nvfamj,nvxlisteppp);\ );\ ,if((s[ppprn-1] <= 0) & (s[ppprn] <= 0),\ resulbriqrn = valeurconstante(0);\ nvsrn = s ;\ for(pppin = ppprn, N-1, nvsrn[pppin] = s[pppin+1]);\ nvfamj = famj ;\ nvfamj[ppprn-1] = 0;\ if(ppprn <= N-1, nvfamj[ppprn] = K);\ for(uplusvrn = 0, -s[ppprn-1]-s[ppprn],\ nvsrn[ppprn-1] = - uplusvrn;\ resulbriqrn = resulbriqrn + \ (sum(urn = max(0,uplusvrn+s[ppprn]), min(uplusvrn,-s[ppprn-1]),\ binomial(-s[ppprn-1],urn) * \ binomial(-s[ppprn],uplusvrn - urn) * \ (K^(-s[ppprn]-(uplusvrn - urn))) * \ ((famj[ppprn-1])^(-s[ppprn-1]-urn))\ ))* briqbn(N-1,ppprn-1,nvsrn,nvfamj,nvxlisteppp);\ );\ , /*i.e. sinon, dans tous les cas restants*/ \ resulbriqrn = valeurconstante(0);\ nvsrn = s ;\ for(pppin = ppprn, N-1, nvsrn[pppin] = s[pppin+1]);\ nvfamj = famj ;\ if(ppprn <= N-1, nvfamj[ppprn] = K);\ for(urn = 1, s[ppprn-1],\ nvsrn[ppprn-1] = urn;\ resulbriqrn = resulbriqrn + \ fctbinoma(urn,s[ppprn-1],s[ppprn]) * \ (famj[ppprn-1]-K)^(urn-s[ppprn-1]-s[ppprn]) *\ briqbn(N-1,ppprn-1,nvsrn,nvfamj,nvxlisteppp);\ );\ nvfamj[ppprn-1] = K;\ for(vrn = 1, s[ppprn],\ nvsrn[ppprn-1] = vrn;\ resulbriqrn = resulbriqrn + \ fctbinoma(vrn,s[ppprn],s[ppprn-1]) * \ (K-famj[ppprn-1])^(vrn-s[ppprn-1]-s[ppprn]) *\ briqbn(N-1,ppprn-1,nvsrn,nvfamj,nvxlisteppp);\ );\ ))));\ resulbriqrn fctbinoma(unoma,enoma,fnoma) = \ (-1)^fnoma * binomial(enoma+fnoma-1-unoma,fnoma-1); /*Pour briqqn : Hyp : N >= 1 et K >= 0 et 1 <= pppqn <= N+1 */ /*On tient cpte seulement de s[t] pour p <= t <= N, idem pour xliste[t] */ /*Ce calcul r\'ecursif prend surement bcp trop de temps.... */ /*Briqqn renvoie une fraction rationnelle qui sera en facteur, et pas qqch de la forme valeurconstante(valeurlamb) */ briqqn(N,pppqn,s,K,xliste) =\ local(resulbriqqn);\ if(pppqn == N+1, \ resulbriqqn = 1;\ ,\ resulbriqqn = sum(kbriqqn = 1, K, xliste[pppqn]^(-kbriqqn) * \ kbriqqn^(-s[pppqn]) * briqqn(N,pppqn+1,s,kbriqqn,xliste));\ );\ resulbriqqn /* ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- Maintenant : fichier sorokin_andco (initialement, appelait briques ci-dessus) ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- */ /* \r briques */ /*implemente jusqu'en profondeur 9 ; ne depend pas de wmax Voir les conditions d'utilisation dans l'en-tete de briques... */ /* Sorokin pour zeta(3) (resp. pour zeta(2) et zeta(4)) sont implementes. Les vraies valeurs sont calculees dans les memes cas. */ evaluation(tabaeval) = \ local(tempcompev);\ xvar_1 = 1;\ xvar_2 = 1;\ xvar_3 = 1;\ xvar_4 = 1;\ xvar_5 = 1;\ xvar_6 = 1;\ xvar_7 = 1;\ xvar_8 = 1;\ xvar_9 = 1;\ tempcompev = eval(tabaeval);\ xvar_1 = 'xvar_1;\ xvar_2 = 'xvar_2;\ xvar_3 = 'xvar_3;\ xvar_4 = 'xvar_4;\ xvar_5 = 'xvar_5;\ xvar_6 = 'xvar_6;\ xvar_7 = 'xvar_7;\ xvar_8 = 'xvar_8;\ xvar_9 = 'xvar_9;\ tempcompev pochhammer(nbpochh,indicepochh) =\ prod(ipochh = 0, indicepochh-1, nbpochh+ipochh) /*Int. Sorokin pour zeta(3) : x^{n+1} \int \produit / (x-uv)^{n+1} (xy-uvw)^{n+1} en utilisant le lemme 14 de Jacky p. 23 (on omet le facteur (-1)^(3n) qui est zerbi...) Ici : k_1 = xvarsoul_1 et k_2 = xvarsoul_2 x = xvar_1 et y = xvar_2 */ /* Integrale de Sorokin pour pi^2, dans le cas de prof 2, telle que page 3/11 _{27} et aussi dans texte de Tanguy du 5/11/2004 C'est (5.23) page 1835 de Sorokin */ sorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) antisorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_2 - xvarsoul_1 + 1, exposantn) *\ pochhammer(xvarsoul_1 - exposantn, exposantn);\ calculder(2,polynomenumsoropid,[0,exposantn+1],[exposantn,exposantn],[2,2],[1,1]) synsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(1, exposantn) *\ pochhammer(xvarsoul_1 - exposantn, exposantn);\ calculder(1,polynomenumsoropid,[0],[2*exposantn+1],[2],[1]) decouplesorokinpideux(exposantn) = \ antisorokinpideux(exposantn) +\ sorokinpideux(exposantn) -\ synsorokinpideux(exposantn) /*Necessite d_n ^2 d_{2n}^2 */ varsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^3 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, 0);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) /*Suffit d_n ^4, semble t il*/ devarsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^4 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, 0) *\ pochhammer(xvarsoul_2 - exposantn, 0);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) /*Necessite d_{2n}^4, semble t il*/ /*Se calcule tres vite...*/ dualsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - xvarsoul_1 + 1, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) antidualsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - xvarsoul_1 + 1, exposantn);\ calculder(2,polynomenumsoropid,[0,exposantn+1],[exposantn,exposantn],[2,2],[1,1]) bidousorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn)\ * xvarsoul_1;\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) espoirsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_2 + xvarsoul_1 - exposantn , 2*exposantn) *\ pochhammer(xvarsoul_2 - exposantn, 0);\ calculder(2,polynomenumsoropid,[0,exposantn+1],[exposantn,exposantn],[2,2],[1,1]) poidstroissorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, 0);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,1],[1,1]) /*semble donner denom raisonnable, mais tous les MZV y sont*/ poidscinqsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^3 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,3],[1,1]) /*semble donner denom raisonnable, mais tous les MZV y sont*/ poidssixsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^4 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[3,3],[1,1]) /*semble donner denom raisonnable, mais tous les MZV y sont*/ poidssixbissorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^3 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn)^2 *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[3,3],[1,1]) /*semble donner denom raisonnable, mais tous les MZV y sont*/ translasorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 7, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) /*Resultats : vvz(ta) = ta[1] + ta[2]*zeta(2) + ta[4]*zeta(4); conv(ta) = \ local(nvta);\ nvta = vector(4);\ nvta[1] = ta[1];\ nvta[2] = ta[2]/6;\ nvta[4] = ta[4]/90;\ nvta (%R) gp > sorokinpideux(0) time = 10 ms. %270 = [0, 0, 0, 3/4, 0, 0, 0, 0, 0] (%R) gp > sorokinpideux(1) time = 200 ms. %258 = [42, -28, 0, 15/4, 0, 0, 0, 0, 0] soit 0.0005582546664279969573891975742 (%R) gp > sorokinpideux(2) time = 690 ms. %259 = [6293/8, -532, 0, 327/4, 0, 0, 0, 0, 0] soit 0.0000007926290829531084335675172203 (%R) gp > sorokinpideux(3) time = 2,790 ms. %260 = [5239423/216, -147448/9, 0, 9951/4, 0, 0, 0, 0, 0] soit 0.000000001481265019715439362 (%R) gp > sorokinpideux(4) time = 17,490 ms. %261 = [3067476167/3456, -5395880/9, 0, 364503/4, 0, 0, 0, 0, 0] soit 3.200045894347102560 E-12 avec 3456 = 2^7 * 3^3 (%R) gp > sorokinpideux(5) time = 2mn, 28,370 ms. %263 = [78183902320027/2160000, -5501167252/225, 0, 14863515/4, 0, 0, 0, 0, 0] soit 7.56153686 E-15 avec 2160000 = 2^7 * 3^3 * 5^4 (%R) gp > sorokinpideux(6) time = 16mn, 21,620 ms. %132 = [410668013880341/259200, -9631798876/9, 0, 650603775/4, 0, 0, 0, 0, 0] soit 1.90751819 E-17 avec 259200 = 2^7 3^4 5^2 En conclusion, d_n ^4 semble toujours convenir, et d_n^2 d_{2n}^2 semble completement inutile... */ /*Apres avoir compris comment demontrer le denom de Sorokin*/ bisorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2, exposantn+1) *\ pochhammer(xvarsoul_1, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) antibisorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_2 - xvarsoul_1, exposantn+1) *\ pochhammer(xvarsoul_2, exposantn);\ calculder(2,polynomenumsoropid,[0,exposantn+1],[exposantn,exposantn],[2,2],[1,1]) decouplebisorokinpideux(exposantn) = \ antibisorokinpideux(exposantn) +\ bisorokinpideux(exposantn) /*Car sur la diago c'est nul identiquement*/ /*Suffit d_n^2 d_{n+1}^2*/ /* decouplebisorokinpideux(1) %136 = [-155/4, 69, 0, -70, 0, 0, 0, 0, 0] (%R) gp > decouplebisorokinpideux(2) %137 = [-124009/36, 8777/2, 0, -6975/2, 0, 0, 0, 0, 0] (%R) gp > decouplebisorokinpideux(3) %138 = [-291841385/1296, 2490089/9, 0, -212440, 0, 0, 0, 0, 0] decouplebisorokinpideux(4) %142 = [-3950997071533/259200, 167281334/9, 0, -28329875/2, 0, 0, 0, 0, 0] */ okcvsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2, exposantn) *\ pochhammer(xvarsoul_1, exposantn);\ calculder(2,polynomenumsoropid,[exposantn+1,0],[exposantn,exposantn],[2,2],[1,1]) antiokcvsorokinpideux(exposantn) =\ local(polynomenumsoropid);\ polynomenumsoropid = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_2 - xvarsoul_1, exposantn) *\ pochhammer(xvarsoul_2, exposantn);\ calculder(2,polynomenumsoropid,[0,exposantn+1],[exposantn,exposantn],[2,2],[1,1]) decoupleokcvsorokinpideux(exposantn) = \ antiokcvsorokinpideux(exposantn) +\ okcvsorokinpideux(exposantn) /*Car sur la diago c'est nul identiquement*/ /*Suffit d_{n+1}^4, et peut etre meme un peu mieux, e.g. d_{n+1}^3 d_n, voire d_{n+1}^2 d_n^2*/ /*Avec la regularisation foireuse d'avant le 16/11/2003, il faut appliquer calculder en ldvar ; avec la nv reg, il faut l'appliquer en [1,1] */ sorokinzetatrois(exposantn) =\ local(polynomenumsoro);\ polynomenumsoro = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,1],[1,1]) /*Int. Sorokin pour zeta(4) (donne aussi du zeta(2)): x^{n+1} \int \produit / (x-uv)^{n+1} (xy-uvwz)^{n+1} en utilisant la prop. 1 de vas4new.dvi (avec l=k_2 et k=k_1) Ici : k_1 = xvarsoul_1 et k_2 = xvarsoul_2 x = xvar_1 et y = xvar_2 */ sorokinzetaquatre(exposantn) =\ local(polynomenumsoro);\ polynomenumsoro = \ (factorielle(exposantn))^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,2],[1,1]) bouclesorokinzetaquatre(exponmaxi)=\ local(totot);\ for(n=0,exponmaxi,\ totot = sorokinzetaquatre(n);\ write(resultatsoro,"Pour l'exposant ",n," le tableau est ",\ totot," et la valeur numerique est ",vraievaleurduntab(totot));\ ); boucletanguyvariante(exponmaxi)=\ local(totot);\ for(n=0,exponmaxi,\ totot = tanguyvariante_quandDV(n);\ write(resultatsoro,"Pour l'exposant ",n," le tableau est ",\ totot," et la valeur numerique est ",vraievaleurduntab(totot));\ ); /*Verif avec les coeffs p_n^{(1)}(5) de Rivoal-Krattenthaler, Prop. 2, p. 20 (cf. arxiv)*/ pnuncinq(expon) = \ sum(i_1 = 0, expon, \ sum(i_2 = i_1, expon, \ binomial(expon,i_2)^2*binomial(expon+i_2,expon)*\ binomial(expon,i_1)^2*binomial(expon + i_2 - i_1,expon)\ )\ ); comparaison(expon) = \ local(tempcomp);\ tempcomp = sorokinzetaquatre(expon)[4];\ xvar_1 = 1;\ xvar_2 = 1;\ tempcomp = [eval(tempcomp),(-1)^(expon)*(7/4)*pnuncinq(expon)];\ xvar_1 = 'xvar_1;\ xvar_2 = 'xvar_2;\ tempcomp /*Resultats de la variante suivante : calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[3,2],ldvar) ne fait pas apparaitre zeta(3) calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,3],ldvar) fait tout apparaitre calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[1,3],ldvar) fait tout apparaitre calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[3,1],ldvar) fait tout apparaitre, et DV a partir de l'exposant 2 */ sorokinvariante(exposantn) =\ local(polynomenumsoro,tempsovar);\ polynomenumsoro = \ (factorielle(exposantn))^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ tempsovar =\ calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,2],ldvar);\ xvar_1 = 1;\ xvar_2 = 1;\ tempsovar = eval(tempsovar);\ xvar_1 = 'xvar_1;\ xvar_2 = 'xvar_2;\ tempsovar sorokinvariante_quandDV(exposantn) =\ local(polynomenumsoro,tempsovar);\ polynomenumsoro = \ (factorielle(exposantn))^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ tempsovar =\ calculder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[3,1],ldvar);\ tempsovar /* Suggestion Tanguy du 13/11 : S_n=sum_(K\ge M\ge L\ge 1} N(K,L,M)/D(K,L,M) avec N(K,L,M)=n!^2(K-M+1)_n(M-L+1)_n(L-n)_n et D(K,L,M)=(K)_{n+1}^2(M)_{n+1}(L)_{n+1}. Ici : K = xvarsoul_1, M = xvarsoul_2, L = xvarsoul_3 Ceci est la version serieuse. ATTENTION, dans un mail du 18/11/03 Tanguy dit qu'il faut diviser par n! : c'est pas dur, la regle est que la somme des indices des Pochammer du numerateur doit valoir de la somme de ceux du denominateur. Mais on n'a pas tenu cpte de cette Rq ci-dessous. La version erronnee est quand on echange L et M dans les xvarsoul_i. Dans ce cas, ca DV quand exposantn >= 2. Donc peut etre que meme avec exposantn=1 on manque des termes... On obtient [9, -10, 30, -16, 0, 0] avec exposantn = 1, dans ce cas erronne. */ tanguyvariante_quandDV(exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ (factorielle(exposantn))^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - xvarsoul_3 + 1, exposantn) *\ pochhammer(xvarsoul_3 - exposantn, exposantn);\ tempsovar =\ calculder(3,polynomenumtr,[0,0,0],[exposantn,exposantn,exposantn],\ [2,1,1],[1,1,1]);\ tempsovar /*Curiosite : mail de Tanguy du 20/11/03 S_N=sum_{1\le L\le K} n!(K-n)_n(L-n)_n(K+L)_n/(K)_{n+1}^2(L)_{n+1}^2 Ici : K = xvarsoul_1, L = xvarsoul_2 */ curiosite(exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - exposantn, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn) *\ pochhammer(xvarsoul_1 + xvarsoul_2, exposantn);\ tempsovar =\ calculder(2,polynomenumtr,[0,0],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar bouclecuriosite(exponmaxi)=\ local(totot);\ for(n=0,exponmaxi,\ totot = curiosite(n);\ write(resultatcur,"Pour l'exposant ",n," le tableau est ",\ totot," et la valeur numerique est ",vraievaleurduntab(totot));\ ); /*En attendant : mail de Tanguy du 19/11/03 n!^2(K-L+1)_n(L-n)_n/(K)_{n+1}^3(L)_{n+1} Ici : K = xvarsoul_1, L = xvarsoul_2 */ enattendant(exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ factorielle(exposantn)^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ tempsovar =\ calculder(2,polynomenumtr,[0,0],[exposantn,exposantn],\ [3,1],[1,1]);\ tempsovar boucleenattendant(exponmaxi)=\ local(totot);\ for(n=0,exponmaxi,\ totot = enattendant(n);\ write(resultatatt,"Pour l'exposant ",n," le tableau est ",\ totot," et la valeur numerique est ",vraievaleurduntab(totot));\ ); /*Simplifie : mail de Tanguy du 19/11/03 n!^{a+b-1} (L-n)_n/(K)_{n+1}^a(L)_{n+1}^b avec a>= 1 et b >= 1 et a + b <= 4 Ici : K = xvarsoul_1, L = xvarsoul_2 */ simplifie(a,b,exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ factorielle(exposantn)^(a+b-1) * \ pochhammer(xvarsoul_2 - exposantn, exposantn);\ tempsovar =\ calculder(2,polynomenumtr,[0,0],[exposantn,exposantn],\ [a,b],[1,1]);\ tempsovar bouclesimplifie(exponmaxi)=\ local(totot);\ for(a=1,3,\ for(b=1,4-a,\ write(resultatsimplifie,"");\ write(resultatsimplifie,"");\ write(resultatsimplifie,"On teste maintenant a = ",a," et b = ",b);\ write(resultatsimplifie,"");\ write(resultatsimplifie,"");\ for(n=0,exponmaxi,\ totot = simplifie(a,b,n);\ write(resultatsimplifie,"Pour l'exposant ",n," le tableau est ",\ totot," et la valeur numerique est ",vraievaleurduntab(totot));\ ))); /* Non-intersection : mail de Tanguy du 18/11/03 R(K,L)= n! (K-L+1)_n(L-n)_n/(K)_{n+1}^2 (L+n+1)_{n+1} (on somme sur K\ge L\ge 1) (le n! est absent du mail de Tanguy, mais on le rajoute ici) (on prend aussi des Pochhammer n+1 au lieu de n au denom, contrairement au mail de Tanguy) Ici : K = xvarsoul_1, L = xvarsoul_2 */ nonintersec(exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn)*\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn+1],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar bouclenonintersec(exponmaxi)=\ local(totot);\ for(n=0,exponmaxi,\ totot = nonintersec(n);\ write(resultatsimplifie,"Pour l'exposant ",n);\ write(resultatsimplifie,"le tableau est ",totot);\ write(resultatsimplifie,"et la valeur numerique est ",vraievaleurduntab(totot));\ write(resultatsimplifie,"");\ );\ autrenumnonintersec(exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn)*\ pochhammer(xvarsoul_2 - exposantn, exposantn)*\ pochhammer(xvarsoul_1 + exposantn + 1, exposantn);\ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar nonintersec_decalage(exposantn,decalage) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, 0)*\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn-decalage],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar bouclenonintersec_decalage(exponmaxi)=\ local(totot);\ for(n=0,exponmaxi,\ for(dec=-3*n,3*n,\ totot = nonintersec_decalage(n,dec);\ write(resultatsimplifie,"Pour l'exposant ",n," et le decalage ",dec);\ write(resultatsimplifie,"le tableau est ",totot);\ write(resultatsimplifie,"et la valeur numerique est ",vraievaleurduntab(totot));\ write(resultatsimplifie,"");\ );\ );\ /* Nouveaunumerateur: mail de Tanguy du 24/11/03 n! (k-l+1)_n(l-n)_n(k+n+1)_n/(k)_{n+1}^2(l)_{n+1}^2 (on somme sur K\ge L\ge 1) Ici : K = xvarsoul_1, L = xvarsoul_2 */ nouveaunum(exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn)*\ pochhammer(xvarsoul_2 - exposantn, exposantn)*\ pochhammer(xvarsoul_1 + exposantn + 1, exposantn);\ tempsovar =\ calculder(2,polynomenumtr,[0,0],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar bouclenouveaunum(exponmaxi)=\ local(totot);\ for(n=0,exponmaxi,\ totot = nouveaunum(n);\ write(resultatsres,"Pour l'exposant ",n);\ write(resultatsres,"le tableau est ",totot);\ write(resultatsres,"et la valeur numerique est ",vraievaleurduntab(totot));\ write(resultatsres,"");\ );\ /* denomquitue num /(K)_{n+1}^2 (L+n-d)_{n+1}^2 (on somme sur K\ge L\ge 1), avec d=0 pour commencer au moins. Ici : K = xvarsoul_1, L = xvarsoul_2 cf. fichier denomquitue */ denomquitue(a,b,exposantn,d) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ pochhammer(xvarsoul_1, a) *\ pochhammer(xvarsoul_2, b); \ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn-d],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar boucledenomquitue(nfixe)=\ local(totot,n);\ n = nfixe;\ for(d=-4*n,0,\ for(a=0,n,\ for(b=0,2*n,\ totot = denomquitue(a,b,n,d);\ write(resultatsres,"Pour l'exposant ",n," et au numerateur (a,b,c) = (",\ a,",",b,",",c,")");\ write(resultatsres,"le tableau est ",totot);\ write(resultatsres,"et la valeur numerique est ",vraievaleurduntab(totot));\ write(resultatsres,"");\ if(totot[3] != 0, write(resultatsproblemes,a,b,d,totot));\ );\ );\ );\ /* autredenomquitue num /(K)_{n+1}^2 (L+n-d)_{n+1}^2 (on somme sur K\ge L\ge 1), avec d=-1 Ici : K = xvarsoul_1, L = xvarsoul_2 */ autredenomquitue(a,b,c,exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ pochhammer(xvarsoul_1, a) *\ pochhammer(xvarsoul_2, b) * \ pochhammer(xvarsoul_1 - xvarsoul_2 +1, c);\ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn+1],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar boucleautredenomquitue(cmax,nfixe)=\ local(totot,n);\ n = nfixe;\ for(a=0,n,\ for(b=0,2*n,\ for(c=0,min(min(cmax,2*n-a),4*n+1-a-b),\ totot = autredenomquitue(a,b,c,n);\ write(resultatsres,"Pour l'exposant ",n," et au numerateur (a,b,c) = (",\ a,",",b,",",c,")");\ write(resultatsres,"le tableau est ",totot);\ write(resultatsres,"et la valeur numerique est ",vraievaleurduntab(totot));\ write(resultatsres,"");\ );\ );\ );\ /* R_{d,P} = \sum_{K \geq L \geq 1} \frac{P(K,L)}{(K)_{n+1}^2 (L+n-d)_{n+1}^2} \begin{Conj} \label{conjtr} Supposons $d \leq 0$ ou $d > 2n$. Si $P(K,L)$ est de degre $\leq n$ en $K$ et de degre $\leq 2n$ en $L$ alors $R_{d,P}$ est une combinaison lin\'eaire (sur $\Q$) de $1$, $\zeta(2)$ et $\zeta(4)$. \end{Conj} Test de cette conj : */ testdenomquitue(a,b,exposantn,d) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ pochhammer(xvarsoul_1, a) *\ pochhammer(xvarsoul_2, b); \ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn-d],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar testboucledenomquitue(nfixe)=\ local(totot,n);\ n = nfixe;\ for(d=-4*n,0,\ for(a=0,n,\ for(b=0,2*n,\ totot = testdenomquitue(a,b,n,d);\ write(resultatsres,"Pour l'exposant ",n," et au numerateur (a,b,c) = (",\ a,",",b,",",c,")");\ write(resultatsres,"le tableau est ",totot);\ write(resultatsres,"et la valeur numerique est ",vraievaleurduntab(totot));\ write(resultatsres,"");\ if(totot[3] != 0, write(resultatsproblemes,a,b,d,totot));\ );\ );\ );\ essaidenomquitue(a,b,exposantn,d) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ xvarsoul_1 ^a *\ xvarsoul_2 ^b* \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, 1); \ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn-d],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar essaiboucledenomquitue(nfixe)=\ local(totot,n);\ n = nfixe;\ for(d=1,2*n,\ for(a=0,2,\ for(b=0,2,\ totot = essaidenomquitue(a,b,n,d);\ write(resultatsres,"Pour l'exposant ",n," et au numerateur (a,b,c) = (",\ a,",",b,",",c,")");\ write(resultatsres,"le tableau est ",totot);\ write(resultatsres,"et la valeur numerique est ",vraievaleurduntab(totot));\ write(resultatsres,"");\ if(totot[3] != 0, write(resultatsproblemes,a,b,d,totot));\ );\ );\ );\ /* numetdenomquituent num /(K)_{n+1}^2 (L+n-d)_{n+1}^2 (on somme sur K\ge L\ge 1), avec d=0 pour commencer au moins. Ici : K = xvarsoul_1, L = xvarsoul_2 cf. fichier denomquitue */ numetdenomquituent(exposantn) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ pochhammer(xvarsoul_2^2 + 2*xvarsoul_2 +1/2 , 1) * \ pochhammer(xvarsoul_1 - xvarsoul_2 +1, 0);\ tempsovar =\ calculder(2,polynomenumtr,[0,exposantn],[exposantn,exposantn],\ [2,2],[1,1]);\ tempsovar /* ------------------------------------------------------------------------ */ /* VRAIES VALEURS CALCULEES PAR SOMMES DE SERIES ------------------------------------------------- */ vraievaleursorokinzetatrois(exposantn,ldvarvraie) =\ local(polynomenumsoro);\ polynomenumsoro = \ factorielle(exposantn) * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ vraievaleurder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,1],ldvarvraie) vraievaleursorokinzetaquatre(exposantn,ldvarvraie) =\ local(polynomenumsoro);\ polynomenumsoro = \ (factorielle(exposantn))^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - exposantn, exposantn);\ vraievaleurder(2,polynomenumsoro,[0,0],[exposantn,exposantn],[2,2],ldvarvraie) vraievaleurtanguyvariante_quandDV(exposantn,ldvarvraie) =\ local(polynomenumtr,tempsovar);\ polynomenumtr = \ (factorielle(exposantn))^2 * \ pochhammer(xvarsoul_1 - xvarsoul_2 + 1, exposantn) *\ pochhammer(xvarsoul_2 - xvarsoul_3 + 1, exposantn) *\ pochhammer(xvarsoul_3 - exposantn, exposantn);\ vraievaleurder(3,polynomenumtr,[0,0,0],[exposantn,exposantn,exposantn],\ [2,1,1],ldvarvraie);\ /* ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- Maintenant : ajouts permettant de rendre l'algo plus convivial ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------- */ algo(N,Pnumr,rsoul,ssoul,asoul) = \ local(xlisteauxn);\ xlisteauxn = vector(N);\ for(iaux = 1, N, xlisteauxn[iaux] = 1);\ calculder(N,Pnumr,rsoul,ssoul,asoul,xlisteauxn) kv_1 = xvarsoul_1; kv_2 = xvarsoul_2; kv_3 = xvarsoul_3; kv_4 = xvarsoul_4; kv_5 = xvarsoul_5; kv_6 = xvarsoul_6; kv_7 = xvarsoul_7; kv_8 = xvarsoul_8; kv_9 = xvarsoul_9;