case=4 # chosen case 0,..,8 ; CAUTION: case n here is case n+1 in the paper # bits of precision to compute (interval arithmetic) m1, mr, Z1, Zr and eps bits=[8,8,8,8,9,8,8,10,16][case] # for the constants given in the paper #bits=53 # better but the constants are more complicated ##################### # radii and density # ##################### # decreasing sizes of discs allowing a compact packing (same order used everywhere) var('r') all_r=[ AA(ZZ[r](r^4-10*r^2-8*r+9).roots(AA)[0][0]), AA(ZZ[r](r^8 - 8*r^7 - 44*r^6 - 232*r^5 - 482*r^4 - 24*r^3 + 388*r^2 - 120*r + 9).roots(AA)[2][0]), AA(ZZ[r](8*r^3+3*r^2-2*r-1).roots(AA)[0][0]), AA(sqrt(2)-1), AA((2*sqrt(3)+1-2*sqrt(1+sqrt(3)))/3), AA(sin(pi/12)/(1-sin(pi/12))), AA((sqrt(17)-3)/4), AA(2/sqrt(3)-1), AA(5-2*sqrt(6)) ] # corresponding density/pi (algebraic value) all_d_opt=[ ZZ[r](27*r^4 + 112*r^3 + 62*r^2 + 72*r - 29).roots(AA)[1][0], ZZ[r](r^8 - 4590*r^6 - 82440*r^5 + 486999*r^4 - 1938708*r^3 + 2158839*r^2 - 1312200*r + 243081).roots(AA)[2][0], ZZ[r](1024*r^3 - 692*r^2 + 448*r - 97).roots(AA)[0][0], ZZ[r](2*r^2 - 4*r + 1).roots(AA)[0][0], ZZ[r](944784*r^4 - 3919104*r^3 - 2191320*r^2 - 1632960*r + 757681).roots(AA)[0][0], ZZ[r](144*r^4 + 9216*r^3 + 133224*r^2 - 127104*r + 25633).roots(AA)[2][0], ZZ[r](4096*r^4 + 2924*r^2 - 289).roots(AA)[1][0], ZZ[r](108*r^2 + 288*r - 97).roots(AA)[1][0], ZZ[r](144*r^4 - 4162200*r^2 + 390625).roots(AA)[2][0]] r=RIF(all_r[case]) d_opt=RIF(all_d_opt[case]*pi) ######################### # some useful functions # ######################### # area of a triangle with sides of length a, b and c area = lambda a,b,c: sqrt((a+b+c)*(a+b-c)*(a-b+c)*(b+c-a))/4 # angle opposite to edge a angle = lambda a,b,c: arccos((b^2+c^2-a^2)/(2*b*c)) # area of the triangle with sides of length a, b and c covered by the discs of radius ra, rb and rc cov = lambda a,b,c,ra,rb,rc: angle(a,b,c)*ra^2/2+angle(b,c,a)*rb^2/2+angle(c,a,b)*rc^2/2 # excess of the same triangle E=lambda a,b,c,ra,rb,rc: RIF(d_opt*area(a,b,c)-cov(a,b,c,ra,rb,rc)) # lexicographic enumeration of nonnegative integer vectors x such that, weighted by u, the sum is at most U def odometer(u,U): x=[0 for i in u] z=0 while true: yield [i for i in x] i=0 while iU: z-=x[i]*u[i] x[i]=0 i+=1 if i==len(u): raise StopIteration x[i]+=1 z+=u[i] ##################################################### # radius of the support circle and distance to edge # ##################################################### def radius(a,b,c,ra,rb,rc): S=(a+b+c)*(a+b-c)*(a-b+c)*(b+c-a) A=4*((ra^2-ra*rb-ra*rc+rb*rc)*a^2 + (rb^2-ra*rb+ra*rc-rb*rc)*b^2 + (rc^2+ra*rb-ra*rc-rb*rc)*c^2)-S B=2*(a^2*(a^2*ra-b^2*(ra + rb)+(2*ra^3-ra^2*rb-ra*rb^2-ra^2*rc+rb^2*rc-ra*rc^2+rb*rc^2)) + b^2*(b^2*rb-c^2*(rb + rc)+(2*rb^3-rb^2*rc-rb*rc^2-rb^2*ra+rc^2*ra-rb*ra^2+rc*ra^2)) + c^2*(c^2*rc-a^2*(rc + ra)+(2*rc^3-rc^2*ra-rc*ra^2-rc^2*rb+ra^2*rb-rc*rb^2+ra*rb^2))) C=a^2*b^2*c^2 + a^2*ra^4 + b^2*rb^4 + c^2*rc^4 + (a^4 - a^2*b^2 - a^2*c^2)*ra^2 - (a^2*b^2 - b^4 + b^2*c^2 + (a^2 + b^2 - c^2)*ra^2)*rb^2 + (c^4 - (a^2 + b^2)*c^2 - (a^2 - b^2 + c^2)*ra^2 + (a^2 - b^2 - c^2)*rb^2)*rc^2 if (A.contains_zero() or C.contains_zero()): x=4*A*C/B^2 return (-C/B*(1+x) if x.upper()<0.78 else RIF(-Infinity,Infinity)) else: D=sqrt(4*S*(a+rb-rc)*(a-rb+rc)*(b+ra-rc)*(b-ra+rc)*(c+ra-rb)*(c-ra+rb)) return ((-B+D)/(2*A) if C<0 else (-B-D)/(2*A)) # signed distance to edge a of the center of the support circle of radius rs def dist(a,b,c,ra,rb,rc,R): S=(a+b+c)*(a+b-c)*(a+c-b)*(b+c-a) return ((b^2+c^2-a^2 + 2*R*(rb+rc-2*ra)+rb^2+rc^2-2*ra^2)*a^2 + (2*R+rb+rc)*(rb-rc)*(b^2-c^2))/(2*a*sqrt(S)) ######################################## # vertex potentials in tight triangles # ######################################## # base vertex potential : 6 values V111 Vrrr Vr1r V1rr V1r1 V11r to be defined # M is the matrix which encodes conditions on base vertex potentials # V111 Vrrr Vr1r V1rr V1r1 V11r # first condition: the sum of vertex potential equal the excess in each of the 4 triangles 111, rrr, 1rr, 11r M=matrix([ [3,0,0,0,0,0], # 3*V111=E111 [0,3,0,0,0,0], # 3*Vrrr=Errr [0,0,1,2,0,0], # Vr1r+2*V1rr=E1rr [0,0,0,0,1,2]]) # V1r1+2*V11r=E11r # second condition: the sum of potentials around a circle is zero (case dependent) small_corona=[ [0,0,0,2,3,0], # 3*V1r1+2*V1rr [0,1,0,2,2,0], # Vrrr+2*V1r1+2*V1rr [0,0,0,4,1,0], # 4*V1rr+V1r1 [0,0,0,0,4,0], # 4*V1r1 [0,2,0,2,1,0], # 2*Vrrr+2*V1rr+V1r1 [0,1,0,4,0,0], # Vrrr+4*V1rr [0,0,0,2,2,0], # 2*V1rr+2*V1r1 [0,0,0,0,3,0], # 3*V1r1 [0,1,0,2,1,0] # Vrrr+2*V1rr+V1r1 ] # large coronas are actually redundant - this is just for information large_corona=[ [0,0,1,0,0,6], # 6*V11r+Vr1r [2,0,1,0,0,4], # 4*V11r+2*V111+Vr1r [0,0,4,0,0,4], # 4*Vr1r+4*V11r [0,0,0,0,0,8], # 8*V11r [0,0,3,0,0,6], # 6*V11r+3*Vr1r [0,0,12,0,0,0], # 12*Vr1r [0,0,2,0,0,8], # 8*V11r+2*Vr1r [0,0,0,0,0,12], # 12*V11r [0,0,6,0,0,12] # 12*V11r+6*Vr1r ] M=matrix(M.rows()+[small_corona[case]]) # 5 conditions on 6 parameters -> add the apex of r1r or 1r1 as a parameter (choice s.t. M is invertible) if matrix(M.rows()+[[0,0,1,0,0,0]]).rank()==6: # case<>5: Vr1r=0 M=matrix(M.rows()+[[0,0,1,0,0,0]]) else: # case 5 : V1r1=0 M=matrix(M.rows()+[[0,0,0,0,1,0]]) ##### c5 (case=4) ############# if case==4: M=matrix([ [3,0,0,0,0,0], [0,2,0,0,0,0], # V'rrr=Errr/2 [0,0,1,2,0,0], [0,0,0,0,1,2], [0,2,0,2,1,0], # r-corona 11rrr : 2*V'rrr+2*V1rr+V1r1 [0,0,1,0,0,0]]) ########################## ME=vector([ E(2,2,2,1,1,1), # E111 E(2*r,2*r,2*r,r,r,r), # Errr E(1+r,1+r,2*r,r,r,1), # E1rr E(r+1,r+1,2,1,1,r), # E11r 0, # 0 around a small circle 0 # 0 for Vr1r (case<>5) or V1r1 (case=5) ]) # parametrized base potentials # in order V111 Vrrr Vr1r V1rr V1r1 V11r V2=M.inverse()*ME # convenient to use a dictionary pos={ (1,1,1):0, (r,r,r):1, (r,1,r):2, (1,r,r):3, (r,r,1):3, (1,r,1):4, (1,1,r):5, (r,1,1):5 } # case=4 (c5) : Vrrr is always at least 0 (and Errr/2 in the center of an r-corona 11rrr) V=lambda i,j,k: RIF(V2[pos.get((i,j,k))]) if case<>4 or (i,j,k)<>(r,r,r) else 0 ############################# # lower bounds on m1 and mr # ############################# # A(i,j,k): tight angle in j A=lambda i,j,k: angle(i+k,j+i,j+k) # L(ri,rj,rk): lower bound on angle in rj in any triangle which can appears in a FM-triangulation of a saturated packing (ri, rj and rk: sizes of circles) L=lambda ri,rj,rk: RIF(arcsin(min(ri/(rj+ri+2*r),rk/(rj+rk+2*r)))) # does exist a corona with a/b/c contacts between circles of sizes 1 and 1 / 1 and r / r and r coronable=lambda a,b,c: mod(b,2)==0 and (a*c==0 or b>0) and (a+b+c>=3) # smallest mx such that positivity around a vertex x is ensured def smallest_m(x): m=RIF(0) # checks all the ways to have angles around x with sum at most 2pi for k in odometer([L(i,x,k) for (i,k) in [(1,1),(1,r),(r,r)]],RIF(2*pi)): if coronable(*k): # lower bound on m given by the configuration k deno=RIF(abs(2*pi-(k[0]*A(1,x,1)+k[1]*A(1,x,r)+k[2]*A(r,x,r)))) if not deno.contains_zero(): mk=-(k[0]*V(1,x,1)+k[1]*V(1,x,r)+k[2]*V(r,x,r))/deno m=max(m,mk) else: # must be non-negative (except for the x-corona in the target packing, theoretically ensured) if case==4 and x<>1 and k==[1,2,2]: # singular vertex c5 : Vrrr=Errr/2 print(k,(k[0]*V(1,x,1)+k[1]*V(1,x,r)+k[2]*E(2*r,2*r,2*r,r,r,r)/2)) else: print(k,(k[0]*V(1,x,1)+k[1]*V(1,x,r)+k[2]*V(r,x,r))) return RealIntervalField(bits)(m).upper() m1=smallest_m(1) mr=smallest_m(r) ######################### # Capping the potential # ######################### # increase mq to satisfy condition (5) of Prop. 2 m1=max(m1,max([-RIF(V(x,1,y)/A(x,1,y)) for (x,y) in [(1,1),(1,r),(r,r)]]).upper()) mr=max(mr,max([-RIF(V(x,r,y)/A(x,r,y)) for (x,y) in [(1,1),(1,r),(r,r)]]).upper()) Z1=RealIntervalField(bits)(-2*pi*min([V(x,1,y)/A(x,1,y) for (x,y) in [(1,1),(1,r),(r,r)]])).upper() Zr=RealIntervalField(bits)(-2*pi*min([V(x,r,y)/A(x,r,y) for (x,y) in [(1,1),(1,r),(r,r)]])).upper() Z=lambda x:Z1 if x==1 else Zr ####################### # Vertex potential Uv # ####################### def Uv(a,b,c,ra,rb,rc): f=lambda x:m1 if x==1 else mr u =min(Z(ra),RIF(V(ra,rb,rc)+f(ra)*abs(angle(a,b,c)-angle(rb+rc,ra+rc,ra+rb)))) u+=min(Z(rb),RIF(V(rb,rc,ra)+f(rb)*abs(angle(b,c,a)-angle(ra+rc,ra+rb,rb+rc)))) u+=min(Z(rc),RIF(V(rc,ra,rb)+f(rc)*abs(angle(c,a,b)-angle(ra+rb,rb+rc,ra+rc)))) # c5: V(r,r,r)=0 as above (we do not know their exact values) but their sum is equal to Es if case==4 and (ra,rb,rc)==(r,r,r): u+=E(2*r,2*r,2*r,r,r,r) return u ###################### # upper bound on eps # ###################### var('a b c ra rb rc') var('M1 Mr') # partial derivative in e of the excess over eps-tight triangles dE=lambda ra,rb,rc,eps,e: RIF(derivative(d_opt*area(a,b,c)-cov(a,b,c,ra,rb,rc),e).subs({a:RIF(rb+rc,rb+rc+eps),b:RIF(ra+rc,ra+rc+eps),c:RIF(ra+rb,ra+rb+eps)})) # upper bound on the partial derivative in e of the vertex potential over eps-tight triangles def dUv(ra,rb,rc,eps,e): # weigth associated to each vertex of the triangle, depending on the radius of the disc at this vertex f=lambda x:M1 if x==1 else Mr # to substitute edge length to edge name after derivation dico={a:RIF(rb+rc,rb+rc+eps),b:RIF(ra+rc,ra+rc+eps),c:RIF(ra+rb,ra+rb+eps)} # sum over the three vertices z =f(ra)*abs(derivative(angle(a,b,c),e).subs(dico)) z+=f(rb)*abs(derivative(angle(b,c,a),e).subs(dico)) z+=f(rc)*abs(derivative(angle(c,a,b),e).subs(dico)) return z # for given m1 and mr, test whether eps suits def test_eps(m1,mr,eps): for (ra,rb,rc) in [(1,1,1),(1,1,r),(1,r,r),(r,r,r)]: for e in [a,b,c]: du=RIF(dUv(ra,rb,rc,eps,e).subs({M1:m1,Mr:mr})) de=RIF(dE(ra,rb,rc,eps,e).subs({M1:m1,Mr:mr})) if du.is_NaN() or de.is_NaN() or du.upper()>de.lower(): return false return true # refine interval eps to have an as large as possible eps def dicho(m1,mr,eps): if eps.absolute_diameter()<2^-(40): return RealIntervalField(bits)(eps).lower() if test_eps(m1,mr,eps.center()): return dicho(m1,mr,eps.bisection()[1]) else: return dicho(m1,mr,eps.bisection()[0]) eps=dicho(m1,mr,RIF(0,1)) ##################### # edge potential Ue # ##################### # parameters found "by hand" (with the help of curves) # le and qe for e in (1,1), (1,r), (r,r) # intuitively, good parameters are such that: # - U=Uv+Ue stays below E (mandatory!) # - for stretched, U slightly below E # - U is more or less "regular" MLQ=[ [(2.7,0.1),(2.3,0.1),(1.9,0.1)], [(2.6,0.2),(2.1,0.2),(1.6,0.2)], [(2.6,0.1),(2.2,0.2),(1.65,0.1)], [(2.5,0.15),(1.8,0.2),(1.2,0.2)], [(2.4,0.05),(1.8,0.05),(1.1,0.07)], [(2.5,0.2),(1.75,0.2),(1,0.2)], [(2.4,0.08),(1.6,0.05),(0.8,0.1)], [(2.24,0.02),(1.33,0.015),(0.44,0.02)], [(2.17,0.02),(1.21787,0.015),(0.285729,0.02)] ] LQ={(1,1):MLQ[case][0],(1,r):MLQ[case][1],(r,1):MLQ[case][1],(r,r):MLQ[case][2]} # Ue edge potential def Ue(a,b,c,ra,rb,rc,R): u=RIF(0) (l,q)=LQ.get((rb,rc)) if a>l: u+=q*dist(a,b,c,ra,rb,rc,R) (l,q)=LQ.get((ra,rc)) if b>l: u+=q*dist(b,c,a,rb,rc,ra,R) (l,q)=LQ.get((ra,rb)) if c>l: u+=q*dist(c,a,b,rc,ra,rb,R) return u ############################# # test on all the triangles # ############################# # the four triangles, with maximal range of edge length in a staturated FM-triangulation T=[(RIF(y+z,y+z+2*r),RIF(x+z,x+z+2*r),RIF(x+y,x+y+2*r),x,y,z) for (x,y,z) in [(1,1,1),(1,1,r),(1,r,r),(r,r,r)]] # an FM-triangle should at least contain half a small disc -> lower bound on 16 times its squared area Smin=RIF(16*((pi/2)*r^2)^2).lower() # return numbers of (good, quasi-tight, non feasible) triangles def test(a,b,c,ra,rb,rc): if a<=rb+rc+eps and b<=ra+rc+eps and c<=ra+rb+eps: # quasi-tight triangle return vector([0,1,0]) S=(a+b+c)*(a+b-c)*(a+c-b)*(b+c-a) # 16 times the squared area if S=0: # otherwise the precision is too low -> need to refine R=radius(a,b,c,ra,rb,rc) if R.lower()>=r.upper(): # not feasible because the packing cannot be saturated return vector([0,0,1]) # one can assume that R is in [0,r] (otherwise not feasible and surely ruled out) R=RIF(max(R.lower(),0),min(R.upper(),r.upper())) e=E(a,b,c,ra,rb,rc) u=Uv(a,b,c,ra,rb,rc)+Ue(a,b,c,ra,rb,rc,R) if e>u: # good triangle return vector([1,0,0]) if e refine v=vector([0,0,0]) for (a2,b2,c2) in cartesian_product([x.bisection() for x in [a,b,c]]): v+=test(a2,b2,c2,ra,rb,rc) return v print('m1='+str(m1)) print('mr='+str(mr)) print('Z1='+str(Z1)) print('Zr='+str(Zr)) print('eps='+str(eps)) import time t0=time.time() l=[test(*i) for i in T] t1=time.time() print(l) print('%d triangles in %d min %d s'%(sum([sum(i) for i in l]),floor((t1-t0)/60),mod(floor(t1-t0),60)))