import numpy as np, matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
ns=[1,2,3,4,6,8,9,12]; ex={1:(0,0),2:(1,0),3:(0,1),4:(2,0),6:(1,1),8:(3,0),9:(0,2),12:(2,1)}; sig=0.55
def P(t2,t3): return sum(n**-sig*np.exp(-1j*(ex[n][0]*t2+ex[n][1]*t3)) for n in ns)
g=np.linspace(0,2*np.pi,500); T2,T3=np.meshgrid(g,g); F=np.abs(P(T2,T3)); bad=F<0.25
L2,L3=np.log(2),np.log(3)
fig=plt.figure(figsize=(16,5.2))
def panel(k,x0,x1,m,tit):
    ax=fig.add_subplot(1,3,k); x=np.linspace(x0,x1,m); a=(x*L2)%(2*np.pi); b=(x*L3)%(2*np.pi)
    ax.imshow(F,origin='lower',extent=[0,2*np.pi,0,2*np.pi],cmap='viridis',alpha=0.85)
    ax.contourf(T2,T3,bad.astype(float),levels=[0.5,1.5],colors=['red'],alpha=0.75)
    if x1-x0<200:
        br=np.where((np.abs(np.diff(a))>3)|(np.abs(np.diff(b))>3))[0]+1
        for sa,sb in zip(np.split(a,br),np.split(b,br)): ax.plot(sa,sb,color='white',lw=1.0)
    else: ax.scatter(a[::8],b[::8],s=0.15,color='white',alpha=0.35)
    i=np.clip((a/(2*np.pi)*499).astype(int),0,499); j=np.clip((b/(2*np.pi)*499).astype(int),0,499)
    ax.scatter(a[bad[j,i]],b[bad[j,i]],s=4,color='yellow',zorder=5)
    ax.set_xlabel(r'$\theta_2=x\log 2\ (\mathrm{mod}\ 2\pi)$'); ax.set_ylabel(r'$\theta_3=x\log 3\ (\mathrm{mod}\ 2\pi)$'); ax.set_title(tit,fontsize=11)
    for f in (ax.set_xticks,ax.set_yticks): f([0,np.pi,2*np.pi])
    ax.set_xticklabels(['0',r'$\pi$',r'$2\pi$']); ax.set_yticklabels(['0',r'$\pi$',r'$2\pi$'])
panel(1,1000,1060,6000,'(a) two primes, window $x\\in[1000,1060]$')
panel(2,1000,9000,400000,'(b) two primes, $x\\in[1000,9000]$')
ax=fig.add_subplot(1,3,3)
pr=np.array([p for p in range(2,200) if all(p%q for q in range(2,int(p**0.5)+1))],float); lp=np.log(pr)
rng=np.random.default_rng(1); X0=1e4
for L,c in [(60,'C0'),(600,'C1'),(6000,'C3')]:
    ds=list(range(1,31)); med=[]
    for d in ds:
        dx=0.04/lp[d-1]; x=np.arange(X0,X0+L,dx); tg=rng.uniform(0,2*np.pi,(40,d)); best=np.full(40,np.inf)
        for k in range(0,len(x),20000):
            ph=(x[k:k+20000,None]*lp[None,:d])%(2*np.pi)
            for q in range(40):
                dd=np.abs(ph-tg[q]); best[q]=min(best[q],np.minimum(dd,2*np.pi-dd).max(1).min())
        med.append(np.median(best))
    ax.plot(ds,med,'o-',color=c,ms=3,label='window length %d'%L)
ax.axhline(np.pi,color='gray',ls=':'); ax.text(1,np.pi*1.02,r'maximal distance $\pi$',color='gray')
ax.set_xlabel('number of primes (dimension of the torus)'); ax.set_ylabel('median distance from curve to a random point')
ax.set_title('(c) the curve of a finite window in higher dimension',fontsize=11); ax.legend(loc='lower right'); ax.set_ylim(0,3.3); ax.grid(alpha=0.3)
plt.tight_layout(); plt.savefig('torus.png',dpi=150)