# restart; read `G:/My Drive/Aek/Coupon/Coupon.txt`; with(combinat,choose): with(combinat,powerset): ############################## # Functions: # # Section 2: Main results # RecMax(S,n), GenMax(S,r,n), THM1(S,r,n), # RecMin(S,n), GenMin(S,r,n), Cor3(S,r,n), # # Section 3: Continuous version # ContMax(S,n), ContMin(S,n), # # Section 4: Probability # RecPL(s,T), GenPL(S,r,n), Prop4(s,T), # RecPW(s,T), GenPW(S,r,n), Cor5(s,T), # ############################## ############################## # Section 2: Main results ############################## # Input: The list S of number of # still missing coupons of each player, # numeric n the total number of coupons # needed to be collected. # Output: The average number of rounds # (each player opens one box per round) # to be opened until all of the players # complete their collections. # Try: # RecMax([10,10],10); RecMax := proc(S,n) option remember; local s,p,P,k,T,left,right,term; if {op(S)} = {0} or min(op(S)) < 0 then return(0); fi: k := nops(S); left := 1-mul( 1-s/n,s in S); right := 1; for P in (powerset(k) minus {{}}) do T := S; for p in P do T[p] := T[p]-1; od: term := 1; for p from 1 to k do if member(p,P) then term := term*S[p]/n; else term := term*(1-S[p]/n); fi: od: right := right+term*RecMax(sort(T),n); od: right/left; end: # Input: The list S of number of # still missing coupons of each player, # numeric r which is number of leading # terms we want to calculate and # symbolic n the total number of coupons # needed to be collected. # Output: The average number of rounds # (each player opens one box per round) # ,as a function of n, to be opened # until all of the players complete # their collections. # Try: # GenMax([2,2,3],2,n); # seq(GenMax([s,2,3],2,n),s=2..10); # evalf(GenMax([60,60],23,n)); GenMax := proc(S,r,n) option remember; local i,j,s,k,S1,Soy,del,R,MyF,P,Q,T,p; if {op(S)} = {0} or min(op(S)) < 0 then return(0); fi: k := nops(S); S1 := add(s, s in S); del := expand(mul(1-S[i]/n,i=1..k)-(1-S1/n)); R := expand(n/S1*del); MyF := expand(n/S1*add(R^i,i=0..r)); #print(S,"MyR",R,"MyF",MyF); Soy := 1; for P in (powerset(k) minus {{}}) do Q := {seq(i,i=1..k)} minus P; T := S; for p in P do T[p] := T[p]-1: od: Soy := Soy+mul(S[i]/n,i in P) *mul(1-S[j]/n, j in Q)*GenMax(sort(T),r,n); od: #print("Soy",S,Soy); Soy := expand(MyF*Soy); add(coeff(Soy,n,2-i)*n^(2-i),i=1..r); end: # Input: The list S of number of # still missing coupons of each player, # numeric r which is number of leading # terms we want to calculate (works only # for r=1 or 2) and symbolic n the total # number of coupons needed to be collected. # Output: The pre-calculated formulas of # average number of rounds ,as a function # of n, to be opened until all of the players # complete their collections. # Try: # seq(THM1([s,2,3],1,n),s=2..10); # A:=[seq(GenMax([s,2,3],2,n),s=2..10)]; # B:=[seq(THM1([s,2,3],2,n),s=2..10)]; # A-B; THM1 := proc(S,r,n) option remember; local i,s,t,x,H1,H2,all,po; all := add(s, s in S); po := add(mul(S[i], i in C),C in choose(nops(S),2)); H1 := x -> add(1/i,i=1..x); #H2 := x -> add(1/i^2,i=1..x); if r = 1 then n*H1(all); elif r = 2 then n*H1(all)-(H1(all)-1)*po/all/(all-1); else ERROR("NoCaseYet"); fi: end: # Section 2.4: Answer to Question 2 # Digits:=30: # A:=evalf([seq(THM1([n,n],2,n),n=30000..30002)]); # B:=evalf([seq(n*log(n)+(log(2)+gamma)*n-n*log(n)/2/(2*n-1) # -n*(log(2)+gamma-1)/2/(2*n-1)+1/4,n=30000..30002)]); # A-B; ############################################# # Input: The list S of number of # still missing coupons of each player, # numeric n the total number of coupons # needed to be collected. # Output: The average number of rounds # (each player opens one box per round) # to be opened until some player # completes the collection. # Try: # RecMin([10,10],10); # A:=evalf(subs(n=170,[seq(GenMin([s,s],2,n),s=30..32)])); # B:=[seq(evalf(RecMin([s,s],170)),s=30..32)]; # B-A; RecMin := proc(S,n) option remember; local s,p,P,k,T,left,right,term; if min(op(S)) = 0 then return(0); fi: k := nops(S); left := 1-mul( 1-s/n,s in S); right := 1; for P in (powerset(k) minus {{}}) do T := S; for p in P do T[p] := T[p]-1; od: term := 1; for p from 1 to k do if member(p,P) then term := term*S[p]/n; else term := term*(1-S[p]/n); fi: od: right := right+term*RecMin(sort(T),n); od: right/left; end: # Input: The list S of number of # still missing coupons of each player, # numeric r which is number of leading # terms we want to calculate and # symbolic n the total number of coupons # needed to be collected. # Output: The average number of rounds # (each player opens one box per round) # ,as a function of n, to be opened # until the FIRST player completes the # collection. # Try: # seq(GenMin([s,2,3],2,n),s=2..10); GenMin := proc(S,r,n) option remember; local i,j,s,k,Soy,Left,Lead,F1,MyF,P,Q,T,p; if min(op(S)) = 0 then return(0); fi: k := nops(S); Lead := add(s, s in S)/n; Left := expand(1-mul(1-S[i]/n,i=1..k)); F1 := expand(1-Left/Lead); MyF := expand(add(F1^i,i=0..r)/Lead); Soy := 1; for P in (powerset(k) minus {{}}) do Q := {seq(i,i=1..k)} minus P; T := S; for p in P do T[p] := T[p]-1: od: Soy := Soy+mul(S[i]/n,i in P) *mul(1-S[j]/n, j in Q)*GenMin(T,r,n); od: Soy := expand(MyF*Soy); add(coeff(Soy,n,2-i)*n^(2-i),i=1..r); end: # Input: The list S of number of # still missing coupons of each player, # numeric r which is number of leading # terms we want to calculate (works only # for r=1 or 2) and symbolic n the total # number of coupons needed to be collected. # Output: The pre-calculated formulas of # average number of rounds ,as a function # of n, to be opened until the first player # completes the collection. # Try: # Cor3([10,10],2,10); Cor3 := proc(S,r,n) option remember; local i,j,C,k,M,T,P; k := nops(S); P := 0; for i from 1 to k do M :=0; for C in choose(k,i) do T := [seq(S[j],j in C)]; M := M+THM1(T,r,n); od: P := P + (-1)^(i-1)*M; od: return(P); end: # Section 2.4: Answer to Question 2 # Digits:=30: # A:=evalf([seq(Cor3([n,n],2,n),n=30000..30002)]); # B:=evalf([seq(n*log(n)-(log(2)-gamma)*n+n*log(n)/2/(2*n-1) # +n*(log(2)+gamma-1)/2/(2*n-1)+3/4,n=30000..30002)]); # A-B; ############################## # Section 3: Continuous version ############################## # Input: The list S of number of # still missing coupons of each player, # symbolic n the total number of coupons # needed to be collected # Output: Values generated from the # recurrence (2) # Note: the output should be the same as # the leading term from THM1 # Try: # seq(ContMax([s,2,3],n),s=2..10); # seq(THM1([s,2,3],1,n),s=2..10); ContMax := proc(S,n) option remember; local i,s,k,all; if {op(S)} = {0} or min(op(S)) < 0 then return(0); fi: k := nops(S); all := add(s, s in S); add(S[i]/all*ContMax([op(1..i-1,S),S[i]-1,op(i+1..k,S)],n) ,i=1..k)+n/all; end: # Input: The list S of number of # still missing coupons of each player, # symbolic n the total number of coupons # needed to be collected # Output: Values generated from the # recurrence (2) (with initial # condition for Min(E)) # Note: the output should be the same as # the leading term from Cor3 # Try: # seq(ContMin([s,2,3],n),s=2..10); # seq(Cor3([s,2,3],1,n),s=2..10); ContMin := proc(S,n) option remember; local i,s,k,all; if min(op(S)) = 0 then return(0); fi: k := nops(S); all := add(s, s in S); add(S[i]/all*ContMin([op(1..i-1,S) ,S[i]-1,op(i+1..k,S)],n),i=1..k)+n/all; end: ############################## # Section 4: Probability ############################## # Input: numeric s, the number of card # missing for the considered player, # the list T, the list of # number of cards missing for the # rest of the player # Output: Prob lost i.e. last place # calculated from the recurrence # Try: # seq(RecPL(s,[5],n),s=1..10); # A:=RecPL(4,[2,3],n): # taylor(A,n=infinity,4); # GenPL([4,2,3],4,n); RecPL := proc(s,T,n) option remember; local i,t,k,all,R,right,left,p,P,term; if s>=0 and {op(T)} = {0} then return(1); elif s=0 or min(op(T)) = -1 then return(0); fi: all := s+add(t, t in T); k := nops(T)+1; left := 1-(1-s/n)*mul( 1-t/n,t in T); right := 0; for P in (powerset(k) minus {{}}) do R := [s,op(T)]; term := 1; for p from 1 to k do if member(p,P) then term := term*R[p]/n; else term := term*(1-R[p]/n); fi: od: for p in P do R[p] := R[p]-1; od: right := right+term*RecPL(R[1], sort([op(2..nops(R),R)]),n); od: factor(right/left); end: # Input: The list S of number of # still missing coupons of each player, # numeric r which is number of leading # terms we want to calculate and # symbolic n the total number of coupons # needed to be collected. # Output: The probability (as a # function of n) that the # first player will finish last. # Try: # seq(GenPL([s,2,3],2,n),s=2..10); # A:=[seq(coeff(GenPL([s,2],2,n),n,-1),s=2..20)]; GenPL := proc(S,r,n) option remember; local i,j,s,k,Soy,Left,Lead,F1,MyF,P,Q,T,p; if {op(1..nops(S),S)} = {0} then return(1); elif S[1]=0 or min(op(2..nops(S),S)) = -1 then return(0); elif {op(2..nops(S),S)} = {0} then return(1); fi: k := nops(S); Lead := add(s, s in S)/n; Left := expand(1-mul(1-S[i]/n,i=1..k)); F1 := expand(1-Left/Lead); MyF := expand(add(F1^i,i=0..r)/Lead); Soy := 0; for P in (powerset(k) minus {{}}) do Q := {seq(i,i=1..k)} minus P; T := S; for p in P do T[p] := T[p]-1: od: Soy := Soy+mul(S[i]/n,i in P) *mul(1-S[j]/n, j in Q)*GenPL(T,r,n); od: Soy := expand(MyF*Soy); add(coeff(Soy,n,1-i)*n^(1-i),i=1..r); end: # Input: numeric s, the number of card # missing for the considered player, # the list T, the list of number of # cards missing for the rest of the # player # Output: Formula of the leading # term of prob lost i.e. last place # Try: # evalf(RecPL(6,[4,4],70000)); # evalf(Prop4(6,[4,4])); Prop4 := proc(s,T) option remember; local t; s/(s+add(t, t in T)); end: #################################################### # Input: numeric s, the number of card # missing for the considered player, # the list T, the list of # number of cards missing for the # rest of the player # Output: Prob wins i.e. first place # calculated from the recurrence # Try: # seq(RecPW(s,[5],n),s=1..10); # A:=RecPW(4,[2,3],n): # taylor(A,n=infinity,4); # GenPW([4,2,3],4,n); RecPW := proc(s,T,n) option remember; local i,t,k,all,R,right,left,p,P,term; if s=0 then return(1); elif min(op(T)) = 0 then return(0); fi: all := s+add(t, t in T); k := nops(T)+1; left := 1-(1-s/n)*mul( 1-t/n,t in T); right := 0; for P in (powerset(k) minus {{}}) do R := [s,op(T)]; term := 1; for p from 1 to k do if member(p,P) then term := term*R[p]/n; else term := term*(1-R[p]/n); fi: od: for p in P do R[p] := R[p]-1; od: right := right+term*RecPW(R[1], sort([op(2..nops(R),R)]),n); od: factor(right/left); end: # Input: The list S of number of # still missing coupons of each player, # numeric r which is number of leading # terms we want to calculate and # symbolic n the total number of coupons # needed to be collected. # Output: The probability (as a # function of n) that the # first player will finish first. # Try: # seq(GenPW([s,2,3],2,n),s=2..10); GenPW := proc(S,r,n) option remember; local i,j,s,k,Soy,Left,Lead,F1,MyF,P,Q,T,p; if S[1] = 0 then return(1); elif min(op(2..nops(S),S)) = 0 then return(0); fi: k := nops(S); Lead := add(s, s in S)/n; Left := expand(1-mul(1-S[i]/n,i=1..k)); F1 := expand(1-Left/Lead); MyF := expand(add(F1^i,i=0..r)/Lead); Soy := 0; for P in (powerset(k) minus {{}}) do Q := {seq(i,i=1..k)} minus P; T := S; for p in P do T[p] := T[p]-1: od: Soy := Soy+mul(S[i]/n,i in P) *mul(1-S[j]/n, j in Q)*GenPW(T,r,n); od: Soy := expand(MyF*Soy); add(coeff(Soy,n,1-i)*n^(1-i),i=1..r); end: # Input: numeric s, the number of card # missing for the considered player, # the list T, the list of number of # cards missing for the rest of the # player # Output: Formula of the leading # term of prob lost i.e. last place # Try: # evalf(RecPW(6,[4,4],70000)); # evalf(Cor5(6,[4,4])); Cor5 := proc(s,T) option remember; local i,j,C,k,P; k:= nops(T); P := 0; for i from 0 to k do P := P + (-1)^i* add(s/(s+add(T[j], j in C)) ,C in choose(k,i)); od: return(P); end: ####################################### # See the recurrence relation # Output: the value of the terms # Try: # Term([5,7],-1); Term := proc(S,pow) option remember; local A,n; A := GenMax(S,2-pow,n); coeff(A,n,pow); end: # Try: # NewS([5,7,3],{1,3}); NewS := proc(S,R) option remember; local r,T; T := S; for r in R do T[r] := T[r]-1; od: T; end: # Solution to PDE version of A_{-t} # Try: # SolU(-1,[x,y]); SolU := proc(t,var) option remember; local m,i,x,s,X,S,Sol,d,A; X := add(x, x in var); if t = -1 then return(ln(X)); fi: Sol := 0; for m from 1 to t+1 do for S in choose(var,m+1) do A := mul(s, s in S)*diff(SolU(t-m,var),S); d := degree(A,ln(X)); A := A*ln(X)/(d+1); Sol := Sol+A; od: od: factor(Sol); end: