# restart; read `G:/My Drive/Aek/Card3New/GenMoment.txt`; with(combinat): # First version: April 20, 2022 # This program accompanies # the paper "No-feedback Card # Guessing Game: Moments and # Distributions under the Optimal Strategy" # by Tipaluck Krityakierne, Poorich, Chaloemkiat # and myself # Sections go according # to the paper. ################################# # Section 2: Introduction # to Card Guessing: # Data of n cards # k shufflings empirically # # Functions: # Decks(n,k), ShufDeck(D), # CountMul(S), PosA(b,S), # # EmM(n,k), BestGuess(n,k), # ########################### # Input: positive integer n and k # Output: the set of deck of n cards # after riffle shuffling k times # Try: # Decks(3,1); # Decks(4,4); Decks := proc(n,k) option remember; local i,j,a,p,B,S,T; S := [[seq(i,i=1..n)]]; T := [1]; for i from 1 to k do B := [0$nops(T)]; for j from 1 to nops(T) do for a in ShufDeck(S[j]) do p := PosA(a[1],S); if p > 0 then B := [op(1..p-1,B),B[p]+a[2]*T[j], op(p+1..nops(B),B)]; else S := [op(S),a[1]]; B := [op(B),a[2]*T[j]]; fi: od: od: T := B; od: {seq( [S[i],T[i]],i=1..nops(T))}; end: # Input: Deck D as a list # Output: the set of all # possibilities after one # riffle shuffling of D. # Try: # ShufDeck([1,2,3]); ShufDeck := proc(D) option remember; local i,c,t,n,A,S,T; n := nops(D); T := []; for t from 0 to n do for S in choose(n,t) do A := [0$n]; for i from 1 to nops(S) do A[S[i]] := D[i]; od: c := t+1; for i from 1 to n do if A[i] = 0 then A[i] := D[c]; c := c+1; fi: od: T := [op(T), A]; od: od: CountMul(sort(T)); end: # Input: List of possibilities, S # Output: Multiplicity of # each deck in the list S # Try: # S:=[[1,2],[1,2],[1,2],[2,1]]: # CountMul(S); CountMul := proc(S) option remember; local i,T,c; T := {}; c := 1; for i from 2 to nops(S) do if S[i]=S[i-1] then c := c+1; else T := T union {[S[i-1],c]}; c := 1; fi: od: T union {[S[nops(S)],c]}; end: # Input: positive integer b # and the list S # Output: the position of b # in the list S, if no b then # return 0 # Try: # PosA(3,[1,4,2,3,5]); PosA := proc(b,S) option remember; local i; for i from 1 to nops(S) do if b = S[i] then return(i); fi: od: return(0); end: ###################################### # Empirical matrix of size n-by-n # Input: positive integers n and k # Output: the list of size n-by-n # on frequency of number j # in position i after k-time # riffle shuffles # Try: # EmM(2,2); # EmM(4,1); EmM := proc(n,k) option remember; local i,s,A; A := [[0$n]$n]; for s in Decks(n,k) do for i from 1 to nops(s[1]) do A[i][s[1][i]] := A[i][s[1][i]]+s[2]; od: od: Matrix(A); end: # Input: positive integer n and k # Output: list of numbers # as the most occuring in each # position for k-shuffling of n cards. # Try: # BestGuess(10,1); # [seq(BestGuess(11,1)[i],i=1..11)]; BestGuess := proc(n,k) option remember; local i,j,m,M,S; S := Array([{}$n]); M := EmM(n,k); for i from 1 to n do m := max(M[i]); for j from 1 to n do if M[i][j] = m then S[i] := S[i] union {j}; fi: od: od: S; end: ############################### # Section 2.1: Generating # function empirial results # # Functions: # GenSlow(n,q), # MyBG1(n), MatchPre(A,B), # ############################### # Empirical generating function # for 1-shuffling # Input: positive integer n and # symbolic q # Output: generating function F # counting the number of correct # guess in each element # Try: # GenSlow(3,q); # seq(print(GenSlow(n,q)),n=1..10); # seq(subs(q=1,diff(GenSlow(n,q),q))/2^n,n=1..10); # seq(Mean1(n),n=1..10); GenSlow := proc(n,q) option remember; local P,A; A := add(P[2]*q^MatchPre(P[1],MyBG1(n)) ,P in Decks(n,1)); factor(sort(A)); end: # Formula of best guess S # for 1-shuffle # Input: positive integer n # Output: Best guess S # Try: # MyBG1(10); MyBG1 := proc(n) option remember; local i,A,B; A := [seq(floor(i/2)+1,i=1..ceil(n/2))]; B := [seq(n+1-A[floor(n/2)+1-i],i=1..floor(n/2))]; [op(A),op(B)]; end: # Number of correct guess of A to B # Input: list A and B of same length # Output: number of match element of # A and B # Try: MatchPre([1,2,4,3],[1,2,3,4]); MatchPre := proc(A,B) option remember; local i,cor; cor := 0; for i from 1 to nops(A) do if A[i]=B[i] then cor := cor+1; fi: od: cor; end: ################################ # Section 3.1: Fast and Faster # generating functions # # Functions: # GenFast(n,q), # RecFast(a1,a2,s,q), # Kro(a,b), # # GenFastest(n,q), # FA(h,q), RecFastest(a1,a2,q), # ################################ # Calculate generating function # by breaking to two pieces # This is very fast!!! # Input: positive integer n and # symbolic q # Output: generating function F # Try: # GenFast(12,q); GenSlow(12,q); GenFast := proc(n,q) option remember; local h,a,b,F; h := ceil(n/2); F := add(add( RecFast(a,h-a,n-h-b+a+1,q) *RecFast(b,n-h-b,h-a+b+1,q) ,b=0..n-h),a=0..h); sort(expand(F)); end: # Rec of G (half interval): # Generating function for # first half (or equivalently # second half) without counting # but by recurrence!! # Input: a1,a2 the length of # first and second sequence in the first half # and starting number of second sequence s, # with symbolic q # Output: Corresponding generating function # Try: # RecFast(3,2,5,q); RecFast := proc(a1,a2,s,q) option remember; local h,S; if a1 >= s then ERROR(); elif a1 = 0 and a2 = 0 then return(1); fi: h := floor((a1+a2)/2)+1; S := 0; if a1 > 0 then S := S+q^Kro(h,a1)*RecFast(a1-1,a2,s,q); fi: if a2 > 0 then S := S+q^Kro(h,s+a2-1)*RecFast(a1,a2-1,s,q); fi: expand(S); end: # Input: positive integers a, b # Output: 1 if a=b and 0 otherwise # (Kronecker delta function) # Try: Kro(13,10); Kro := proc(a,b) option remember; if a=b then 1; else 0; fi: end: ############################################ # New faster idea!! # Calculate generating function # by breaking to two pieces # Input: positive integer n and # symbolic q # Output: generating function F # Try: # GenFastest(120,q); GenFast(120,q); GenFastest := proc(n,q) option remember; local h,B,F; h := ceil(n/2); F := FA(h,q)*FA(n-h,q); if n >=4 then B := -2*q^2-2*q^3+4*q^4; elif n=3 then B := -q-2*q^2+3*q^3; elif n=2 then B := -2*q+2*q^2; elif n=1 then B := -1+q; else B := 0; fi: expand(F+B); end: # Input: numeric h and symbolic q # Output: Generating function of the # first half of length h. # Try: # FA(42,q); FA := proc(h,q) option remember; local a,A; A := add(RecFastest(a,h-a,q),a=0..h); expand(A); end: # Rec of G (half interval): # Generating function for # first half (or equivalently # second half) without counting # but by recurrence!! # Count only correction of # the first sequence!!! # Input: a1,a2 the length of # first and second sequence in the first half # and starting number of second sequence s, # with symbolic q # Output: Corresponding generating function # Try: # RecFastest(3,2,q); RecFastest := proc(a1,a2,q) option remember; local h,G; if a1<0 or a2 <0 then return(0); elif a1=0 and a2=0 then return(1); fi: h := floor((a1+a2)/2)+1; G := q^Kro(h,a1)*RecFastest(a1-1,a2,q) +RecFastest(a1,a2-1,q); expand(G); end: ########################### # Section 3.2: # Numerical moments for k=1 # # Functions: # NumMoF(n,r), # MoFA(h,r), # ########################### # Numerical moments for # generating function F # Input: numeric n and r # Output: the rth moment, E[X^r], # from the deck of n cards. # Try: # NumMoF(10,0); # [seq(NumMoF(n,1)/2^n,n=1..10)]; NumMoF := proc(n,r) option remember; local i,q,F; F := GenFastest(n,q); for i from 1 to r do F := q*diff(F,q); od: subs(q=1,F); end: # Numeric moment for # generating function FA # (on page 13) # Input: numeric h and r # Output: the rth moment, E[X^r], # from the half deck of length h. # Try: # [seq(MoFA(h,1),h=1..10)]; MoFA := proc(h,r) option remember; local a1,i,FA; FA := add(RecFastest(a1,h-a1,q),a1=0..h); for i from 1 to r do FA := expand(q*diff(FA,q)); od: subs(q=1,FA); end: ########################### # Section 3.3: # Still Numeric results # # Functions: # CountYA(n,h,r), NumMoYA(n,h,r), # NumMoX(n,r), # NumMoXMean(n,r), # MyC(n,r), # ########################### # Corollary 4, page 13 # Numeric value of # counting YA # Input: symbolic n and numeric h,r # Output: Counting of YA^r in equaiton 7 # Try: # CountYA(n,11,2); CountYA := proc(n,h,r) option remember; 2^(n-h)*MoFA(h,r); end: # Numeric value of moments from YA # Input: symbolic n and numeric h,r # Output: Moments from YA^r # Try: # NumMoYA(n,11,2); NumMoYA := proc(n,h,r) option remember; local A; A := CountYA(n,h,r)/2^n; simplify(A); end: # Numeric value of moments of X # Input: numeric n,r # Output: Moments from X^r # (the whole thing) # Try: # A:=[seq(NumMoX(n,1),n=1..15)]; # B:=[seq(NumMoF(n,1)/2^n,n=1..15)]; NumMoX := proc(n,r) option remember; local i,h; h := ceil(n/2); add(binomial(r,i)*NumMoYA(n,h,i)*NumMoYA(n,n-h,r-i) ,i=0..r)+MyC(n,r)/2^n; end: # Numeric value of moments # about the mean of X # Input: numeric n,r # Output: Moments about # the mean of X^r # Try: # [seq(NumMoXMean(n,1),n=1..15)]; # [seq(NumMoXMean(n,2),n=1..15)]; NumMoXMean := proc(n,r) option remember; local i,mu; mu := NumMoX(n,1); add(binomial(r,i)*NumMoX(n,i) *(-mu)^(r-i),i=0..r); end: MyC := proc(n,r) option remember; if n >= 4 then 4*4^r-2*(3^r+2^r); elif n=3 then 3*3^r-2*2^r-1; elif n=2 then 2*2^r-2; elif n=1 then 1; fi: end: ########################### # Section 4: Calculation # for higher Moments ########################### ########################### # Section 4.3: # Fit sub-function from # overlapping stages # # Functions: # FitOverlap(h1,a1,r), # ########################### # Test for fit in SubSumCYr. # Compare this result to # section 4.3 E[X^2] in the paper # Input: binary (0 or 1) of # h1 and a1, and numeric r # Output: Fit SubSumCYr after # r iterations # Try: # Section 4.3, Middle level, Case 1: # FitOverlap(0,0,1,L,S); # Section 4.3, Middle level, Case 2: # FitOverlap(0,1,1,L,S); # A:=FitOverlap(1,0,2,L,S); # [seq(eval(subs({S=i,L=15},A)),i=1..13)]; # [seq(SubSumCYr(2*i,30-1,2),i=1..13)]; FitOverlap := proc(h1,a1,r,LL,ss) option remember; local L,s,eq,var,P,Q,x,a,b; P := x-> add(a[s]*x^s,s=0..ceil(r/2)); Q := x-> add(b[s]*x^s,s=0..floor(r/2)); eq := {seq(seq( P(L-s)*binomial(2*(L-s),L-s)/4^L +Q(L-s)/4^s=SubSumCYr(2*s-a1,2*L-h1,r) ,s=4..L-1),L=4..r+8)}; var := {seq(a[s],s=0..ceil(r/2)) ,seq(b[s],s=0..floor(r/2))}; var := solve(eq,var); # In L and s form subs(var union {x=LL-ss}, P(x)*binomial(2*x,x)/4^LL+Q(x)/4^ss); end: ########################### # Section 4.4: # Numeric values of # Overlapping stages # (sum of Lower Level) # # Functions: # SumCYr(n,h,r), SubSumCYr(a,h,r), # Eq19(n,h,r), # ########################### # Sum of C[Y1Y2..Yr] from # Lower Level and Middle Level # Input: numeric n,h,r # Output: Sum of C[Y1Y2..Yr] with n cards # Try: # [seq(SumCYr(2*h,h,1),h=1..15)]; SumCYr := proc(n,h,r) option remember; local i; 2^n*add( binomial(i-1,floor(i/2)) *SubSumCYr(i,h,r-1),i=1..h); end: # Sub function of Sum of C[Y1Y2..Yr] # Input: numeric i,h,r # Output: Sum of C[Y1Y2..Yr] # starting from position i # Try: # [seq(SubSumCYr(2,h,1),h=1..15)]; SubSumCYr := proc(i,h,r) option remember; local j; if r=0 then return(2^(-i)); fi: add(binomial(j-i-1,floor(j/2)-floor(i/2)-1) *SubSumCYr(j,h,r-1),j=i+1..h); end: # Section 4.4, Proposition 9 # Combine M_m(h) for C[YA^r] # Input: numeric n,h,r # Output: C[YA^r] of n cards # (using calculation in proposition 9) # Try: # [seq(CountYA(2*h,h,4),h=1..15)]; # [seq(Eq19(2*h,h,4),h=1..15)]; Eq19 := proc(n,h,r) option remember; local k; if r=0 then return(2^n); fi: add( k!*stirling2(r,k)*SumCYr(n,h,k),k=1..r); end: ########################### # Section 4.4, 4.5: # Symbolic Guessing # for SumCYr in Corollary 10 # and results in section 4.5 # # Functions: # Cor10(h1,r,n,L), # FitMoFA(h1,r,n,L), # ForMoX(n1,r,L), # SimpleForMoX(n1,r,L), # ForMoXMean(n1,r,L), # ########################### # Input: binary h1 = L mod 2, numeric r # and symbolic n,L # Output: Formula of E[Y_A^r] in n and L # for case h is even or odd # according to corollary 10 # Try: # Cor10(0,2,n,L); # seq(print(Cor10(1,r,n,L)),r=1..5); Cor10 := proc(h1,r,n,L) option remember; local L1,i,a,b,P,Q,eq,var; if member(h1,{0,1}) = false then ERROR("Wrong h"); fi: eq := {}; P := L -> add(a[i]*L^i,i=0..ceil(r/2)); Q := L -> add(b[i]*L^i,i=0..floor(r/2)); for L1 from 3 to r+10 do eq := eq union {P(L1)*binomial(2*L1,L1)/4^L1 +Q(L1)=Eq19(n,2*L1-h1,r)/2^n}; od: var := {seq(a[i],i=0..ceil(r/2)) ,seq(b[i],i=0..floor(r/2))}; var := solve(eq,var); 2^n*subs(var,P(L)*binomial(2*L,L)/4^L+Q(L)); end: # This is equation 21 # Calculate formula for C[Y_A^r] # through fitting polynomial # from data of generating function # Input: binary h1 = L mod 2, numeric r # and symbolic n,L # Output: Formula of E[Y_A^r] in n and L # for case h is even or odd # Try: # seq(FitMoFA(0,r,n,L),r=1..5); # t:=time():A:=FitMoFA(0,200,n,L): time()-t; FitMoFA := proc(h1,r,n,L) option remember; local i,j,a,b,deg1,deg2,A,B,S,eq,var; if member(h1,{0,1}) = false then ERROR("Wrong h"); fi: deg1 := ceil(r/2); deg2 := floor(r/2); A := L -> add(a[i]*L^i,i=0..deg1)*binomial(2*L,L); B := L -> add(b[i]*L^i,i=0..deg2)*4^L; S := [seq(MoFA(2*j-h1,r),j=1..r+5)]; eq := {seq(A(j)+B(j)=S[j],j=1..nops(S))}; var := {seq(a[i],i=0..deg1),seq(b[i],i=0..deg2)}; var := solve(eq,var); 2^(n-h)*subs(var,A(L)+B(L)); end: # Formula of E[X^r] at the end # of section 4.4 # Input: numeric n1 in {-1,0,1,2} i.e. n mod 4 # numeric r and symbolic L # Output: Formula in L of E[X^r] # where n mod 4 = n1 # Note: Good for n >=4. # Try: # ForMoX(0,1,L); # A:=[seq( eval( subs(L=i,ForMoX(0,1,L))),i=1..10)]; # B:=[seq(NumMoX(4*L,1),L=1..10)]; ForMoX := proc(n1,r,L) option remember; local i,n,C,A,B; if member(n1,{-1,0,1,2})= false then ERROR("Bad n1"); fi: n := 4*L+n1; C := 4*4^r-2*(3^r+2^r); if n1 = -1 then A:= add(binomial(r,i)*Cor10(0,i,n,L)/2^n *Cor10(1,r-i,n,L)/2^n,i=0..r)+C/2^n; elif n1 = 0 then A:= add(binomial(r,i)*Cor10(0,i,n,L)/2^n *Cor10(0,r-i,n,L)/2^n,i=0..r)+C/2^n; elif n1 = 1 then A:= add(binomial(r,i)*Cor10(1,i,n,L+1)/2^n *Cor10(0,r-i,n,L)/2^n,i=0..r)+C/2^n; elif n1 = 2 then A:= add(binomial(r,i)*Cor10(1,i,n,L+1)/2^n *Cor10(1,r-i,n,L+1)/2^n,i=0..r)+C/2^n; fi: B := factor(coeff(A,binomial(2*L,L),2))*binomial(2*L,L)^2 + factor(coeff(A,binomial(2*L,L),1))*binomial(2*L,L)^1; B+simplify(A-B); end: # No C # Try: # SimpleForMoX(0,4,L); SimpleForMoX := proc(n1,r,L) option remember; local i,n,A,B; if member(n1,{-1,0,1,2})= false then ERROR("Bad n1"); fi: n := 4*L+n1; if n1 = -1 then A:= add(binomial(r,i)*Cor10(0,i,n,L)/2^n *Cor10(1,r-i,n,L)/2^n,i=0..r); elif n1 = 0 then A:= add(binomial(r,i)*Cor10(0,i,n,L)/2^n *Cor10(0,r-i,n,L)/2^n,i=0..r); elif n1 = 1 then A:= add(binomial(r,i)*Cor10(1,i,n,L+1)/2^n *Cor10(0,r-i,n,L)/2^n,i=0..r); elif n1 = 2 then A:= add(binomial(r,i)*Cor10(1,i,n,L+1)/2^n *Cor10(1,r-i,n,L+1)/2^n,i=0..r); fi: B := sort(factor(coeff(A,binomial(2*L,L),2)))*binomial(2*L,L)^2 + sort(factor(coeff(A,binomial(2*L,L),1)))*binomial(2*L,L)^1; B+sort(simplify(A-B)); end: # Symbolic moment about mean # Formula of E[(X-mu)^r] # in section 4.5 # Input: numeric n1 in {-1,0,1,2} i.e. n mod 4 # numeric r and symbolic L # Output: Formula in L of E[(X-mu)^r] # where n mod 4 = n1 # Note: Good for n >= 4. # Try: # ForMoXMean(0,1,L); # A:=[seq( eval( subs(L=i,ForMoXMean(1,5,L))),i=1..15)]; # B:=[seq(NumMoXMean(4*L+1,5),L=1..15)]: A-B; ForMoXMean := proc(n1,r,L) option remember; local i,mu,A,B; if member(n1,{-1,0,1,2})= false then ERROR("Bad n1"); fi: mu := SimpleForMoX(n1,1,L); A := add(binomial(r,i) *SimpleForMoX(n1,i,L) *(-1*mu)^(r-i),i=0..r); # For cleaner output! B := add( sort(factor(simplify( (4^(2*r+1-i))^(L+n1)* coeff(A,binomial(2*L,L),2*r+1-i)))) *binomial(2*L,L)^(2*r+1-i) /(4^(L+n1))^(2*r+1-i) ,i=1..2*r); print("Ignore the smaller terms"); B+sort(simplify(A-B)); end: