# -*- coding: utf-8 -*- # ============================================================================ # VERIF_CS.py -- verificateur AUTONOME des bornes superieures sur les # constantes de Chvatal-Sankoff gamma_k. # # A LANCER TEL QUEL, sans argument. Les parametres sont fixes ci-dessous. # Aucune dependance a un autre fichier, aucun fichier intermediaire. # Sortie au fil de l'eau. # # Ce que le script fait, dans l'ordre : # (1) reconstruit de zero la machine de canonicite M_h de Lueker # (JACM 56(3), art. 17, 2009, Figure 7) ; # (2) AUTO-TEST : verifie qu'elle reproduit la Table II de Lueker (k=2) ; # (3) pour chaque k demande, calcule la borne et VERIFIE le certificat # A_z v <= lambda v (son eq. 41) composante par composante en # arithmetique rationnelle EXACTE (module fractions) ; # (4) en deduit gamma_k <= 2 log(k z)/log lambda (son Theoreme 3.3). # # Le sens de la rigueur : y = 2 log(kz)/log(lambda) est CROISSANTE en lambda # (les deux logarithmes sont negatifs), donc majorer lambda majore y. Le # certificat de Collatz-Wielandt fournit un MAJORANT de lambda pour tout # vecteur x > 0, quelle qu'en soit la provenance. La borne est donc sure. # ============================================================================ # ------------------------- PARAMETRES (deja fixes) -------------------------- CAS = [(3,3),(4,3),(5,3),(6,3),(7,3),(8,3),(9,3),(10,3),(11,3),(12,3),(13,3),(14,3),(15,3),(3,4)] # (k, h) CAS += [(4,4)] # ~2 min 30 ; duree totale du script : quelques minutes # CAS += [(5,4)] # ~un quart d'heure. Ligne non lancee par defaut ; # # decommenter pour la verifier avec les memes routines # # (11 093 etats, 7 053 339 aretes, certificat exact). DENOM_X = 10**8 # denominateur du vecteur rationnel du certificat ITER_PUIS = 400 # iterations de la methode de la puissance (flottant) # ---------------------------------------------------------------------------- import itertools, time from math import comb, log from fractions import Fraction as F import numpy as np # ---- appariements minimaux (u a, v a) : a absent de u et de v, alphabets # ---- de u et de v disjoints ; tronques a h caracteres de chaque cote def minimal_matches(k, h): out = set() for a in range(k): rest = [c for c in range(k) if c != a] for lu in range(h): for u in itertools.product(rest, repeat=lu): allowed = [c for c in rest if c not in set(u)] for lv in range(h): for v in itertools.product(allowed, repeat=lv): out.add((tuple(u)+(a,), tuple(v)+(a,))) return sorted(out) def wkey(p): return (len(p[0]), len(p[1]), len(set(p[0][:-1])), len(set(p[1][:-1]))) # ---- poids w_{mu,z}(p) : serie generatrice des appariements minimaux dont la # ---- troncature donne p. Si les DEUX cotes sont au plafond, les deux # ---- prolongements ne peuvent pas partager une lettre libre : inclusion- # ---- exclusion, qui ne depend que des cardinaux. def weight(key, h, k, z): lu, lv, na, nb = key f = k - 1 - na - nb if lu < h and lv < h: return z**(lu+lv) if lu == h and lv < h: return z**(lu+lv)/(1-(na+f)*z) if lv == h and lu < h: return z**(lu+lv)/(1-(nb+f)*z) def e(base, pp): return sum(comb(pp,j)*(-1)**(pp-j)/(1-(base+j)*z) for j in range(pp+1)) EA = [e(na,p) for p in range(f+1)]; EB = [e(nb,q) for q in range(f+1)] tot = sum(comb(f,p)*comb(f-p,q)*EA[p]*EB[q] for p in range(f+1) for q in range(f-p+1)) return z**(2*h)*tot # ---- fonction de transition delta (Lueker, Figure 7) def delta(x, y, dx, dy, p, h): a = x+p[0]; b = y+p[1]; na = len(a); nb = len(b) T = [[0]*(nb+1) for _ in range(na+1)] T[h][h] = h for i in range(h,0,-1): T[i-1][h] = T[i][h]-dx[i-1] for j in range(h,0,-1): T[h][j-1] = T[h][j]-dy[j-1] for i in range(na+1): Ti = T[i]; Tm = T[i-1] if i>0 else None; ai = a[i-1] if i>0 else None for j in range(nb+1): if i>h or j>h: if i==0: Ti[j] = Ti[j-1] elif j==0: Ti[j] = Tm[j] elif ai==b[j-1]: Ti[j] = Tm[j-1]+1 else: u = Ti[j-1]; v = Tm[j]; Ti[j] = u if u>v else v i, j = na, nb # test de canonicite : le retro-parcours while i>h or j>h: # canonique doit ramener en (h,h) if T[i][j]==T[i-1][j]: i -= 1 elif T[i][j]==T[i][j-1]: j -= 1 else: i -= 1; j -= 1 if i!=h or j!=h: return None ndx = tuple(T[na-h+t][nb]-T[na-h+t-1][nb] for t in range(1,h+1)) ndy = tuple(T[na][nb-h+t]-T[na][nb-h+t-1] for t in range(1,h+1)) return (a[-h:], b[-h:], ndx, ndy) # ---- formes canoniques : les permutations de l'alphabet sont une symetrie # ---- (l'ECHANGE HAUT/BAS n'en est PAS une : Lueker, p. 17:18) def canon(q): x,y,dx,dy = q; ren = {}; n = 0 for ch in x+y: if ch not in ren: ren[ch] = n; n += 1 return (tuple(ren[c] for c in x), tuple(ren[c] for c in y), dx, dy) def rel(m, p): ren = {}; n = m; out = [] for part in p: cur = [] for ch in part: if ch < m: cur.append(ch) else: if ch not in ren: ren[ch] = n; n += 1 cur.append(ren[ch]) out.append(tuple(cur)) return (out[0], out[1]) def build(k, h): M0 = minimal_matches(k,h) q0 = canon((tuple([0]*h), tuple([0]*h), (1,)*h, (1,)*h)) idx = {q0:0}; order = [q0]; E = {}; qi = 0; CL = {} while qi < len(order): x,y,dx,dy = order[qi]; i0 = qi; qi += 1 m = (max(x+y)+1) if (x+y) else 0 if m not in CL: g = {} for p in M0: r = rel(m,p) if r in g: g[r][1] += 1 else: g[r] = [p,1] CL[m] = [(v[0],v[1]) for v in g.values()] for (p,mult) in CL[m]: r = delta(x,y,dx,dy,p,h) if r is None: continue r = canon(r); j = idx.get(r) if j is None: j = idx[r] = len(order); order.append(r) kk = (i0,j,wkey(p)); E[kk] = E.get(kk,0)+mult return len(M0), len(order), [(a,b,c,v) for (a,b,c),v in E.items()] def matvec(n, r, c, w, x): return np.bincount(r, weights=w*x[c], minlength=n) def analyse(k, h, zfix=None, iters=ITER_PUIS): nsym, n, edges = build(k,h) keys = sorted(set(e[2] for e in edges)); ki = {kk:i for i,kk in enumerate(keys)} r = np.array([e[0] for e in edges]); c = np.array([e[1] for e in edges]) lab= np.array([ki[e[2]] for e in edges]); mul = np.array([e[3] for e in edges], float) def lam(z, it): w = mul*np.array([float(weight(kk,h,k,z)) for kk in keys])[lab] x = np.ones(n) for _ in range(it): y = matvec(n,r,c,w,x); mx = y.max() if mx <= 0: return 0.0, None x = np.maximum(y/mx, 1e-300) y = matvec(n,r,c,w,x) return float((y/x).max()), x if zfix is None: # recherche du z optimal (flottant) best = None for z in np.linspace(0.015, 0.98/k, 22): l,_ = lam(z, 30) if 0 < l < 1: yy = 2*np.log(k*z)/np.log(l) if best is None or yy < best[1]: best = (z,yy) a,b = max(1e-4,best[0]-0.02), min(0.98/k-1e-6, best[0]+0.02) gr = (5**0.5-1)/2 for _ in range(22): cc = b-gr*(b-a); d = a+gr*(b-a) lc,_ = lam(cc,45); ld,_ = lam(d,45) fc = 2*np.log(k*cc)/np.log(lc) if 0= # la valeur publiee (sens de la rigueur) et s'en ecarter tres peu. good = (n==spub) and (float(ls)>=float(lpub)) and (float(ls)-float(lpub)<1e-6) \ and abs(y-float(bpub))<1e-5 ok_all &= good print(f" h={h} : etats={n:5d} (publie {spub:5d}) | lambda={float(ls):.9f} " f"(publie {lpub}) | borne={y:.6f} (publie {bpub}) -> {'OK' if good else 'ECHEC'}", flush=True) print(f"\n AUTO-TEST GLOBAL : {'REUSSI' if ok_all else 'ECHEC'}\n", flush=True) print("="*74, flush=True) print("BORNES SUPERIEURES CERTIFIEES", flush=True) print("="*74, flush=True) DANCIK = {3:0.765803,4:0.708236,5:0.664428,6:0.629316,7:0.600184,8:0.575407, 9:0.553937,10:0.534855,11:0.517842,12:0.502591,13:0.488800, 14:0.476198,15:0.464619} for (k,h) in CAS: t0 = time.time() nsym,n,ne,z0,ls,ok = analyse(k,h) y = 2*log(k*float(z0))/log(float(ls)) ref = DANCIK.get(k) print(f" k={k:2d} h={h} | {nsym:7d} symboles, {n:6d} etats, {ne:8d} aretes", flush=True) print(f" z0={z0} lambda_sup={float(ls):.12f} certificat exact : {ok}", flush=True) print(f" gamma_{k} <= {y:.9f}" + (f" (Dancik 1994 : {ref})" if ref else "") + f" [{time.time()-t0:.1f}s]\n", flush=True)