# restart; read `c:/Users/Asus/Google Drive/Aek/Floor/Floor.txt`; with(combinat): ################################################## # Base on the paper: # The Curious Bounds of Floor Function Sums (2017) # From: Thotsaporn "Aek" Thanatipanonda and # Elaine Wong # Written August 10, 2017 ################################################## # Section 1: Generate F and S directly # Functions: # Formu, GenF, S, # Section 2: Values of S fast! # Functions: # DirectS1, MyS1, # DirectS2, MyS2, RecS2, # MyS3, RecS, # Section 3: Min/Max of S for each n and m # (Observation for Conjectures) # Functions: # MaxS, MaxSK, SubMaxSK, # MinS, MinSK, SubMinSK, # Section 4: Helping with the proofs # of case n=2,3 in Theorem 2.4 and 3.4 # Functions: # Del, CheckDel2, CheckDel3, HelpDel3 # Section 5: Conjecture of the missing bounds # Functions: # ConjBound, CalBound, RecBounds ######################################## # Section 1: Generate F and S directly ######################################## # Calculate the term inside the sum # given by Tverberg's formula directly # Input: numeric n, and symbolic k,m and a # Output: the term inside the sum # Try: Formu(3,k,m,a); Formu := proc(n,k,m,a) option remember; local i,t,C,A; A := [seq(a[i],i=1..n)]; add(add( (-1)^(nops(A)-i)*floor((add(t,t in C)+k)/m) ,C in choose(A,i)),i=0..n); end: # Evavuate Formu at the given numbers # Input: numeric m, set of numbers A and numeric k # Ouput: Evaluate the term at the given point # Note that n = nops(A) # Try: [seq([seq(GenF(11,[a,b],5),a=0..10)],b=0..10)]; GenF := proc(m,A,k) option remember; local i,b; simplify(subs([seq(b[i]=A[i],i=1..nops(A))] ,Formu(nops(A),k,m,b))); end: # Evaluate the main formula S_m # which is the sum of F. # Input: numeric m, set of numbers A and numeric K # Ouput: the sum S_m as mentioned in the paper # Try: S(9,[1,3,4],7); S := proc(m,A,K) option remember; local k; add(GenF(m,A,k),k=0..K); end: ############################################ # Section 2: Calculate the values of S fast! ############################################ # The sum S of one variable # calculated directly from definition 1.1 # Input: numeric m, a and K # Ouput: the sum S_m({a},K) # Try: [seq([seq(DirectS1(25,a,K),K=0..24)],a=0..24)]: DirectS1 := proc(m,a,K) option remember; local k; add( floor((a+k)/m)-floor(k/m) ,k=0..K); end: # Formula for S by Proposition 1.2 # Input: numeric m, a and K # Ouput: the sum S_m({a1},K) # Try: A:=[seq([seq(MyS1(25,a,K),K=0..24)],a=0..24)]: # B:=[seq([seq(DirectS1(25,a,K),K=0..24)],a=0..24)]: # A-B; MyS1 := proc(m,a,K) option remember; floor(a/m)*(K+1) +max(0,a+(K mod m)-(floor(a/m)+1)*m+1) + floor(K/m)*(a mod m); end: # The sum S of two variables # calculated directly from definition 2.1 # Input: numeric m,a,b and K # Ouput: the sum S_m({a,b},K) # Try: DirectS2(10,3,7,6); DirectS2 := proc(m,a,b,K) option remember; local k; add( floor((a+b+k)/m)-floor((a+k)/m)-floor((b+k)/m)+floor(k/m) ,k=0..K); end: # Formula for S by Proposition 2.2 # Input: numeric m,a,b and K # Ouput: the sum S_m({a,b},K) # Try: A:=[seq([seq(MyS2(25,a,11,K),K=0..24)],a=0..24)]: # B:=[seq([seq(DirectS2(25,a,11,K),K=0..24)],a=0..24)]: # A-B; MyS2 := proc(m,a,b,K) option remember; if K >= m then ERROR("NotInRange"); fi: (floor((a+b)/m)-floor(a/m)-floor(b/m))*(K+1) +max(0,((a+b) mod m)+K+1-m) -max(0,(a mod m)+K+1-m)-max(0,(b mod m)+K+1-m); end: # Formula for S recursively # i.e. equation in the proof of Proposition 2.2 # Input: numeric m,a,b and K # Ouput: the sum S_m({a,b},K) # Try: # A:=[seq([seq(RecS2(15,a,11,K),K=0..24)],a=0..24)]: # B:=[seq([seq(DirectS2(15,a,11,K),K=0..24)],a=0..24)]: # A-B; RecS2 := proc(m,a,b,K) option remember; local a1,b1; if K >= m then ERROR("NotInRange"); fi: MyS1(m,a+b,K)-MyS1(m,a,K)-MyS1(m,b,K); end: # The sum S of three variables # Formula for S by Proposition 3.2 # (Much faster than the direct calculation) # Input: numeric m,a,b,c and K # Ouput: the sum S_m({a,b,c},K) # Try: A:=[seq(MyS3(8,2,3,5,K),K=0..7)]: # B:=[seq(S(8,[2,3,5],K),K=0..7)]: # A-B; MyS3 := proc(m,a,b,c,K) option remember; if K >= m then ERROR("NotInRange"); fi: (floor((a+b+c)/m)-floor((a+b)/m) -floor((b+c)/m)-floor((a+c)/m) +floor(a/m)+floor(b/m)+floor(c/m))*(K+1) +max(0,((a+b+c) mod m)+K+1-m) -max(0,((a+b) mod m)+K+1-m) -max(0,((b+c) mod m)+K+1-m) -max(0,((a+c) mod m)+K+1-m) +max(0,(a mod m)+K+1-m) +max(0,(b mod m)+K+1-m) +max(0,(c mod m)+K+1-m); end: # Formula for S of any number # of variables recursively # Input: numeric m, set A={a,b,c,...} and K # Ouput: the sum S_m(A,K) # Try: # A:=[seq([seq(RecS(15,[a,11,4,9,7],K),K=0..14)],a=0..24)]: # B:=[seq([seq(S(15,[a,11,4,9,7],K),K=0..14)],a=0..24)]: # A-B; RecS := proc(m,A,K) option remember; local n,B; n := nops(A); if n =2 then return(MyS2(m,A[1],A[2],K)); fi: B := [op(1..n-2,A)]; RecS(m, sort([op(B),A[n]+A[n-1]]),K) -RecS(m,sort([op(B),A[n-1]]),K) -RecS(m,sort([op(B),A[n]]),K); end: ################################## # Section 3: Find Max/Min of S ################################## # Maximum value of S for some fixed m and n # Input: numeric integers m and n # Output: Maximum values of S_m(A,K) where |A|=n # over 0<= a1,...,an,K <= m-1 # and the set of positions that it occurs # Try: MaxS(12,2); # [seq(MaxS(m,2),m=1..10)]; MaxS := proc(m,n) option remember; local K,CMax,PosMax; CMax := MaxSK(m,0,n)[1]; PosMax := {[op(MaxSK(m,0,n)[2]),0]}; for K from 1 to m-1 do if MaxSK(m,K,n)[1] > CMax then CMax := MaxSK(m,K,n)[1]; PosMax := {[op(MaxSK(m,K,n)[2]),K]}; elif MaxSK(m,K,n)[1] = CMax then PosMax := PosMax union {[op(MaxSK(m,K,n)[2]),K]}; fi: od: return([CMax,PosMax]); end: # Maximum value of S for some fixed m,K and n # Input: numeric integers m,K and n where 0<=K<=m-1. # Output: Maximum values of S_m(A,K) where |A|=n # over 0<= a1,...,an <= m-1 # and the set of positions that it occurs # Try: MaxSK(12,3,2); # [seq(MaxSK(12,K,2),K=0..11)]; MaxSK := proc(m,K,n) option remember; SubMaxSK(m,K,[],n); end: # Sub function of MaxSK for recursive purpose SubMaxSK := proc(m,K,A,s) option remember; local i,a,bou,CMax,PosMax,SS; if nops(A) = 0 then bou := m-1; else bou := A[nops(A)]; fi: if nops(A) < s then SS := [seq(SubMaxSK(m,K,[op(A),a],s),a=0..bou)]; else SS := [seq([RecS(m,A,K),{A}],a=0..bou)]; fi: CMax := SS[1][1]; PosMax := SS[1][2]; for i from 2 to nops(SS) do if SS[i][1] > CMax then CMax := SS[i][1]; PosMax := SS[i][2]; elif SS[i][1] = CMax then PosMax := PosMax union SS[i][2]; fi: od: return([CMax,PosMax]); end: ############################################ # Minimum value of S for some fixed m and n # Input: numeric integers m and n # Output: Minimum values of S_m(A,K) where |A|=n # over 0<= a1,...,an,K <= m-1 # and the set of positions that it occurs # Try: MinS(12,2); # [seq(MinS(m,2),m=1..10)]; MinS := proc(m,n) option remember; local K,CMin,PosMin; CMin := MinSK(m,0,n)[1]; PosMin := {[op(MinSK(m,0,n)[2]),0]}; for K from 1 to m-1 do if MinSK(m,K,n)[1] < CMin then CMin := MinSK(m,K,n)[1]; PosMin := {[op(MinSK(m,K,n)[2]),K]}; elif MinSK(m,K,n)[1] = CMin then PosMin := PosMin union {[op(MinSK(m,K,n)[2]),K]}; fi: od: return([CMin,PosMin]); end: # Minimum value of S for some fixed m,K and n # Input: numeric integers m,K and n where 0<=K<=m-1. # Output: Minimum values of S_m(A,K) where |A|=n # over 0<= a1,...,an <= m-1 # and the set of positions that it occurs # Try: MinSK(12,3,2); # [seq(MinSK(12,K,2),K=0..11)]; MinSK := proc(m,K,n) option remember; SubMinSK(m,K,[],n); end: # Sub function of MaxSK for recursive purpose SubMinSK := proc(m,K,A,s) option remember; local i,a,bou,CMin,PosMin,SS; if nops(A) = 0 then bou := m-1; else bou := A[nops(A)]; fi: if nops(A) < s then SS := [seq(SubMinSK(m,K,[op(A),a],s),a=0..bou)]; else SS := [seq([RecS(m,A,K),{A}],a=0..bou)]; fi: CMin := SS[1][1]; PosMin := SS[1][2]; for i from 2 to nops(SS) do if SS[i][1] < CMin then CMin := SS[i][1]; PosMin := SS[i][2]; elif SS[i][1] = CMin then PosMin := PosMin union SS[i][2]; fi: od: return([CMin,PosMin]); end: ####################################### # Section 4: Helping with the proofs # of case n=2,3 in Theorem 2.4 and 3.4 ####################################### # Delta function defined inside the # proof of theorem 2.4 # Input: integers m,a,b,K # Output: S2(m,a,b,K+1)-S2(m,a+1,b,K) # Try: Del(15,12,4,9); Del := proc(m,a,b,K) option remember; local a1,b1; a1 := a mod m; b1 := b mod m; if a1+b1 < m and b1+K+2-m <=0 then return(0); elif a1+b1 < m and b1+K+2-m > 0 then return(-1); elif a1+b1 >= m and b1+K+2-m <=0 then return(1); elif a1+b1 >= m and b1+K+2-m > 0 then return(0); else ERROR("NotInRange",[m,a,b,K]); fi: end: # Check that the delta function "Del" # is correct for each given m # Input: positive integer m # Output: print list [a,b,K] if Del and # the actual difference are not the same # Try: CheckDel2(17); CheckDel2 := proc(m) option remember; local a,b,K,s,t,ca; for K from 0 to m/2-1 do for a from 0 to 2*m-2 do for b from 0 to 2*m-2 do s := MyS2(m,a,b,K+1)-MyS2(m,a+1,b,K); t := Del(m,a,b,K); if s <> t then print("diff is",s,t, "at", [a,b,K]); fi: od:od:od: return("Done"); end: ##################### # Three Variables ##################### # Check that the formula of square function # (on top of page 10): # S({a,b,c},K)-S({a+1,b,c},K-1)= # Del({a,b+c},K)-Del({a,b},K)-Del({a,c},K) # is correct for each given m # Input: positive integer m # Output: print list [a,b,c,K] if # the actual difference and square # are not the same # Try: CheckDel3(13); CheckDel3 := proc(m) option remember; local a,b,c,K,s,t; for K from 0 to m/2-1 do for a from 0 to m-2 do for b from 0 to a do for c from 0 to b do s := S(m,[a,b,c],K+1)-S(m,[a+1,b,c],K); t := Del(m,a,b+c,K)-Del(m,a,b,K)-Del(m,a,c,K); if s <> t then print("values are",s,t,"at",[a,b,c,K]); fi: od:od:od:od: return("Done"); end: # Verify the value of S_m in # case B of the proof of Theorem 3.4 # Input: positive integer m and ca in {1,2} # Output: the difference between S and the # case formula B1 or B2. # Try: HelpDel3(14,2); HelpDel3 := proc(m,ca) option remember; local a,b,c,K,s,t,S; S :={}; for K from floor(m/3) to floor(m/2)-1 do for a from 1 to m-1 do for b from 0 to a do for c from 0 to b do # Case 1 if ca=1 and a+b+c >=m and b+c+K+1 <= m and a+b < m then s := RecS(m,[a,b,c],K); t := (K+1)-max(0,a+b+K+1-m) -max(0,a+c+K+1-m)+max(0,a+K+1-m); S := S union {s-t}; # Case 2 elif ca=2 and a+b+c >= 2*m and b+c+K+1 <=2*m and c+K+1-m > 0 then s := RecS(m,[a,b,c],K); t := -(K+1)-max(0,a+b+K+1-2*m) -max(0,a+c+K+1-2*m)+(a+b+c)+3*(K+1-m); S := S union {s-t}; fi: od:od:od:od: return(S); end: ############################################## # Section 5: Conjecture of the missing bounds ############################################## # Conjecture the value of the uppper # bound if n is odd and the lower # bound if n is even. # Here we consider the bound to be # in the form a/b, a < b where a,b runs from # 1 to n-1 # Input: positive integer n and symbolic m # Output: maximum values if n is odd, # minimum values if n is even. # Try: ConjBound(3,7); # A:=[0,seq(ConjBound(2*n,m),n=2..36)]; # Guess1D(A,4,2,1,n,N); # A:=[seq(ConjBound(2*n-1,m),n=2..36)]; # Guess1D(A,4,2,2,n,N); # A:=[seq(ConjBound(n,m),n=3..67)]: # Guess1D(A,9,2,3,n,N); ConjBound := proc(n,m) option remember; local a,b,PosMax,PosMin,CMax,CMin; if n<=2 then ERROR("N is too small"); fi: PosMax := {}; PosMin := {}; CMax := -10^6; CMin := 10^6; for a from 1 to n-1 do for b from a+1 to n do if CalBound(n,a/b,m)/m > CMax then CMax := CalBound(n,a/b,m)/m; PosMax := {a/b*m}; elif CalBound(n,a/b,m)/m = CMax then PosMax := PosMax union {a/b*m}; fi: if CalBound(n,a/b,m)/m < CMin then CMin := CalBound(n,a/b,m)/m; PosMin := {a/b*m}; elif CalBound(n,a/b,m)/m = CMin then PosMin := PosMin union {a/b*m}; fi: od:od: if n mod 2 =1 then #[CMax*m,PosMax]; return(CMax); else #[CMin*m,PosMin]; return(CMin); fi: end: # Sub-function of ConjBound # Input: positive integer n, point c where # a1=a2=...=an=m*c and K=m*c-1, # and symbolic m # Output: calculated values of S_m(A,K) # Try: CalBound(3,1/3,m); CalBound := proc(n,c,m) option remember; local t; add( (-1)^(n-t)*floor(t*c)*binomial(n,t),t=0..n)*c*m +add( (-1)^(n-t)*binomial(n,t) *max(0,t*c-floor(t*c)+c-1),t=0..n)*m; end: # Check the recurrence of f(n) # in conjecture 5.2 page 15 # Input: positive integer n >=4 # Output: Value of f(n) # Try: [seq(RecBound(n),n=4..50)]; RecBound := proc(n) option remember; if n <= 1 then ERROR("nTooSmall"); elif n=2 then return(0); elif n=3 then return(1/3); elif n=4 then return(-1); elif n=5 then return(2); elif n=6 then return(-3); elif n=7 then return(8); elif n=8 then return(-18); elif n=9 then return(36); elif n=10 then return(-65); fi: (10*(n^2+n-8)*RecBound(n-1) -4*(2*n^2-10*n+3)*RecBound(n-2)-24*(2*n-11)*RecBound(n-3) -32*(2*n^2-10*n-1)*RecBound(n-4)-192*(n-1)*(n-5)*RecBound(n-5) +64*(2*n^2-22*n+51)*RecBound(n-6)+384*(2*n-13)*RecBound(n-7) -256*(n-3)*(n-8)*RecBound(n-8)+512*(n-9)*(n-8)*RecBound(n-9)) /(-5*(n+3)*(n-2)); end: