# qcode / j.sage

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76``` ```# Askey-Wilson polynomials # (a;q)_n:=prod_{k=0}^n-1 (1-aq^k) def br_n(q,n,a): return prod([1-a*q**k for k in range(n)]) # (a_1,a_2,...,a_k;q)_n:=prod_{j=1}^k(a_j;q)_n def brlist_n(q,n,*a): def br(x): return prod([1-x*q**k for k in range(n)]) return prod([br(aj) for aj in a]) # a and b are sequences def q_hyper(q,z,a,b,ulim=oo): s = len(a) if s-1 == len(b): if ulim==oo: # do a symbolic summation var('kkk') return sum(z**kkk*brlist_n(q,kkk,*a)/brlist_n(q,kkk,*(b+[q])), kkk, 0, ulim) else: return \ reduce(lambda x,y: x+y, \ [z**k*brlist_n(q,k,*a)/brlist_n(q,k,*(b+[q])) for \ k in range(ulim+1)]) # q-binomial coefficient def qbinomial(q,v,j): if j < 0: return 0 return q^(v*j-j*(j-1)/2)*(-1)^j*br_n(q,j,q^(-v))/br_n(q,j,q) def Eber(v,k,i,u): if 2*k<=v: return sum( [(-1)^t*binomial(u,t)*binomial(k-u,i-t)*binomial(v-k-u,i-t) for t in range(i+1)]) return Eber(v,v-k,i,u) # P and Q for J(v,k), the Johnson scheme def JP(v,k): def ff(i,j): return Eber(v,k,i,j) return matrix(k+1,k+1,ff) def JQ(v,k): return binomial(v,k)*JP(v,k)^-1 def delsarte_bound_J(n, k, d, \ isinteger=False, return_data=False, solver="PPL"): from sage.numerical.mip import MixedIntegerLinearProgram if not isinstance(d, list): d = range(1, d) p = MixedIntegerLinearProgram(maximization=True, solver=solver) A = p.new_variable(integer=isinteger) # A>=0 is assumed p.set_objective(sum([A[r] for r in range(k+1)])) p.add_constraint(A[0]==1) for i in d: p.add_constraint(A[i]==0) Q = JQ(n,k) Q = diagonal_matrix([y.denominator() for y in Q])*Q # we can scale, # as the RHS==0 for j in range(1,k+1): rhs = sum([Q[j,r]*A[r] for r in range(k+1)]) p.add_constraint(0*A[0] <= rhs) try: bd=p.solve() except sage.numerical.mip.MIPSolverException, exc: print "Solver exception: ", exc, exc.args if return_data: return A,p,False return False if return_data: return A,p,bd else: return bd ```