congruent conditional prime distribution¶

Hiroyuki OGAWA (Osaka Univ.)¶

The original version of this program was written in Mathematica and made for a talk at the RIMS conference.The original version of this program was written in Mathematica and made for a talk at the RIMS conference on 2013/11/5.¶

In [1]:

import primes.primes as primes
import matplotlib.pyplot as plt
from math import log10

In [2]:

def chop(x, e=3):
    return(int(100*x*10**e+1/2)/10**e)

In [3]:

def ccpd(X, d, m, residue=()):
    if type(m)==int:
        m = (m, m)
    if len(m)<2:
        m = m*2
    if len(residue) < 2:
        residue = (primes.irreducible_residues(m[0]), 
                   primes.irreducible_residues(m[1]))
    prsm = [(p%m[0], p%m[1]) for p in primes.table((X, X+d))]
    ccM = [[0]*m[1] for _ in range(m[0])]
    s = prsm[0]
    for t in prsm[1:]:
        ccM[s[0]][t[1]] += 1
        s = t
    n = {k:sum(ccM[k]) for k in residue[0]}
    return({'modulo':m, 'residue':residue, 'X':(X,d), 
            'data':{k:{j:(ccM[k][j], chop(ccM[k][j]/n[k])) 
                       for j in residue[1]} 
                    for k in residue[0]}})

In [4]:

def listTranspose(a):
    return([list(x) for x in zip(*a)])

In [5]:

def ccpdAll(m, Xs, residue=()):
    if type(m)==int:
        m = (m, m)
    if len(m)<2:
        m = m*2
    if len(residue) < 2:
        residue = (primes.irreducible_residues(m[0]), 
                   primes.irreducible_residues(m[1]))
    wk = {X:ccpd(*X[-2:], m, residue=residue)['data'] for X in Xs}
    return({'modulo':m, 'residue':residue, 'Xs':Xs, 
            'data':[[((i,j),listTranspose([wk[X][i][j] for X in Xs])) 
                     for j in residue[1]] 
                    for i in residue[0]]})

In [6]:

def ccpdPlot(m, Xs, residue=(), size=10):
    if type(m)==int:
        m = (m, m)
    if len(m)<2:
        m = m*2
    if len(residue) < 2:
        residue = (primes.irreducible_residues(m[0]), 
                   primes.irreducible_residues(m[1]))
    wk = ccpdAll(m, Xs, residue=residue)['data']
    idx = [log10(max(1,X[-2])+X[-1]/2) for X in Xs]
    #
    for ds in wk:
        c = 100/len(ds)
        plt.plot(idx, [c]*len(idx), c='red')
        for ij,data in ds:
            plt.scatter(idx, data[1], s=size, 
                        label='{} ({})'.format(ij[1],m[1]))
        plt.yticks([c-5,c-1,c,c+1,c+5], 
                   ['-5%','-1%','mean','+1%','+5%'])
        plt.legend(title='{} ({}) →'.format(ij[0],m[0]),
                   loc="upper left", bbox_to_anchor=(1.02, 1.0,), 
                   borderaxespad=0, handletextpad = 0)
        plt.show()

In [7]:

Xs = sum([[(4**i+j*10**5, 10**5) for j in range(10)] for i in [0,8,9]]
         +[[(4**i+j*10**6, 10**6) for j in range(10)] for i in [0,10,11]]
         +[[(4**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(12, 14))]
         +[[(2**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(28, 51))], [])

In [8]:

Xs8 = sum([[(i, 2**i+2*j*10**8, 10**8) for j in range(5)] for i in range(28, 51)], [])

In [9]:

#Xs = sum([[(4**i+j*10**5, 10**5) for j in range(10)] for i in [0,8,9]]
#         +[[(4**i+j*10**6, 10**6) for j in range(10)] for i in [0,10,11]]
#         #+[[(4**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(12, 14))]
#         #+[[(2**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(28, 51))]
#         , [])

In [10]:

ccpdPlot(4, Xs)

In [11]:

ccpdPlot(6, Xs)

In [12]:

ccpdPlot(5, Xs)

In [13]:

ccpdPlot(7, Xs)

In [14]:

ccpdPlot(8, Xs)

In [15]:

ccpdPlot(9, Xs)

In [16]:

ccpdPlot(11, Xs)

In [17]:

ccpdPlot((4,6), Xs)
ccpdPlot((6,4), Xs)

In [18]:

ccpdPlot((4,5), Xs)
ccpdPlot((5,4), Xs)

In [19]:

ccpdPlot((6,5), Xs)
ccpdPlot((5,6), Xs)

In [20]:

ccpdPlot((5,7), Xs)
ccpdPlot((7,5), Xs)

In [21]:

ccpdPlot((5,7), Xs8)
ccpdPlot((7,5), Xs8)

In [22]:

ccpdPlot((4,8), Xs)
ccpdPlot((8,4), Xs)

In [23]:

ccpdPlot((6,8), Xs)
ccpdPlot((8,6), Xs)

In [24]:

ccpdPlot((6,7), Xs)
ccpdPlot((7,6), Xs)

In [ ]:

In [ ]: