Prime circles

prime_circles.py

#!/usr/bin/env python3
from pprint import pprint
from itertools import compress
from math import sqrt, inf, gcd
from random import randrange
import sys

# prime_circles.py: Find primes that form circles
# https://math.stackexchange.com/questions/2886024/a-conjecture-involving-prime-numbers-and-circles
# https://www.reddit.com/r/math/comments/98ie5z/a_conjecture_involving_prime_numbers_and_circles/

''' Usage:
./prime_circles.py <base> <limit>
Default values:
./prime_circles.py 10 100000
Will not work for very large values because of the floating point math
used when calculating the circle center and radius.

Known counterexamples (N = 1_000_000):

Base 3:
    (5, 11)

Base 4:
    (5, 13)
    (5, 37)
    (7, 11)
    (7, 23)
    (7, 59)
    (11, 19)
    (11, 31)
    (13, 29)
    (19, 23)
    (19, 47)
    (23, 43)

Base 8:
    (11, 19)
    (11, 43)
'''

# Values set in main function
PRIME_LIMIT = None
BASE = None
prime_y = None
all_primes_greater_than = None
all_the_primes = None

# https://stackoverflow.com/a/46635266
def rwh_primes1v1(n):
    """ Returns  a list of primes < n for n > 2 """
    sieve = bytearray([True]) * (n//2)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = bytearray((n-i*i-1)//(2*i)+1)
    return [2,*compress(range(3,n,2), sieve[1:])]

# All primes greater than 7
def all_primes_up_to(n):
    return [x for x in rwh_primes1v1(n) if x > all_primes_greater_than]

# https://stackoverflow.com/a/50974391
def define_circle(p1, p2, p3):
    """
    Returns the center and radius of the circle passing the given 3 points.
    In case the 3 points form a line, returns (None, infinity).
    """
    temp = p2[0] * p2[0] + p2[1] * p2[1]
    bc = (p1[0] * p1[0] + p1[1] * p1[1] - temp) / 2
    cd = (temp - p3[0] * p3[0] - p3[1] * p3[1]) / 2
    det = (p1[0] - p2[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p2[1])

    if abs(det) < 1.0e-6:
        return (None, inf)

    # Center of circle
    cx = (bc*(p2[1] - p3[1]) - cd*(p1[1] - p2[1])) / det
    cy = ((p1[0] - p2[0]) * cd - (p2[0] - p3[0]) * bc) / det

    radius = sqrt((cx - p1[0])**2 + (cy - p1[1])**2)
    return ((cx, cy), radius)

def verify_circle(a, b, c, d):
    c1, r1 = define_circle(a, b, c)
    c2, r2 = define_circle(d, b, c)

    if abs(r1 - r2) < 1e-6:
        return True
    else:
        return False

''' (notes for base 10)
Trivial case: x/10 == y/10
Since (11, 13, 17, 19) are all prime, any pair in the form
(n*10 + a, n*10 + b) forms a circle with (10 + a, 10 + b).

Similarly, (17, 29) are (1, 2) units apart, (1 * 10 + 2 == 12)
so any pair which is separated by (-1 * 10 + 2 == 8) units
forms a circle with the initial pair, for example (59, 67).
So if a % 10 == 7 and abs(a - b) == 12, we can pair (a, b) with
(59, 67). This can be extended to larger numbers, perhaps as a form
of cache. (unimplemented)

Algorithm for the general case:
    Choose a third prime number and form a circle
    Check intersects with the lines y=1, y=3, y=7, y=9 for primes
    If no primes found, repeat.
'''

# https://gist.github.com/bnlucas/5857478
def miller_rabin(n, k=10):
    if n == 2 or n == 3:
        return True
    if not n & 1:
        return False

    def check(a, s, d, n):
        x = pow(a, d, n)
        if x == 1:
            return True
        for i in range(s - 1):
            if x == n - 1:
                return True
            x = pow(x, 2, n)
        return x == n - 1

    s = 0
    d = n - 1

    while d % 2 == 0:
        d >>= 1
        s += 1

    for i in range(k):
        a = randrange(2, n - 1)
        if not check(a, s, d, n):
            return False
    return True

def is_prime(n):
    if n < PRIME_LIMIT:
        return n in all_the_primes
    else:
        return miller_rabin(n)

def circle4primes(a, b):
    ax = a // BASE
    ay = a % BASE
    bx = b // BASE
    by = b % BASE

    # Trivial case
    if BASE == 10:
        if ax == bx:
            if ax == 1:
                return (100 + ay, 100 + by)
            else:
                return (10 + ay, 10 + by)

    # General case
    for c in all_the_primes:
        # We cannot assume that a < b < c < d
        if c == a or c == b:
            continue
        #print("-------------------------------")
        #print("Checking",a,",",b,",",c,":")
        cx = c // BASE
        cy = c % BASE
        dd = check_intersects((ax, ay), (bx, by), (cx, cy))
        #print("Checking",a,",",b,",",c,":", dd)
        if dd is None:
            continue
        #elif len(dd) > 1:
        #    return (c, dd[0])
        else:
            return (c, dd)

def check_intersects(a, b, c):
    ax = a[0]
    ay = a[1]
    bx = b[0]
    by = b[1]
    cx = c[0]
    cy = c[1]

    center, r = define_circle(a, b, c)
    #print("Form a circle with center",center)
    #print("and radius",r)
    if center is None:
        return None

    for y in prime_y:
        x1, x2 = circle_intersects_with_y(center, r, y)
        #print("Got intersect with y =",y,":",x1,",",x2)
        if x1 is None:
            pass
        elif x1 * BASE + y < all_primes_greater_than:
            pass
        elif (x1, y) == (ax, ay) or (x1, y) == (bx, by) or (x1, y) == (cx, cy):
            pass
        elif not is_prime(x1 * BASE + y):
            pass
        elif verify_circle((ax, ay), (bx, by), (cx, cy), (x1, y)):
            return x1 * BASE + y

        if x2 is None:
            pass
        elif x2 * BASE + y < all_primes_greater_than:
            pass
        elif (x2, y) == (ax, ay) or (x2, y) == (bx, by) or (x2, y) == (cx, cy):
            pass
        elif not is_prime(x2 * BASE + y):
            pass
        elif verify_circle((ax, ay), (bx, by), (cx, cy), (x2, y)):
            return x2 * BASE + y

def circle_intersects_with_y(center, r, y):
    cx = center[0]
    cy = center[1]
    # Find x such that dx^2 + dy^2 == r^2
    # dx = abs(cx-x)
    dy = abs(cy-y)
    dx2 = r*r - dy*dy
    if dx2 < 0:
        # No intersection
        return (None, None)

    dx = sqrt(dx2)
    x1 = cx-dx
    x2 = cx+dx
    # Round to nearest integer
    ix1 = int(x1+0.5)
    ix2 = int(x2+0.5)
    #print("x1:",x1,",ix1:",ix1)
    #print("x2:",x2,",ix2:",ix2)

    # Check that intersect points are integers, otherwise set them to None
    if abs(x1-ix1) > 1.0e-6:
        ix1 = None
    if abs(x2-ix2) > 1.0e-6:
        ix2 = None
    return (ix1, ix2)

#p1 = [x for x in all_the_primes if x % 10 == 1]
#p3 = [x for x in all_the_primes if x % 10 == 3]
#p7 = [x for x in all_the_primes if x % 10 == 7]
#p9 = [x for x in all_the_primes if x % 10 == 9]

if __name__ == "__main__":
    # Is this how constants work?
    BASE = 10 if len(sys.argv) < 2 else int(sys.argv[1])
    PRIME_LIMIT = 100000 if len(sys.argv) < 3 else int(sys.argv[2])

    # For base 10, check y=[1, 3, 7, 9], for other bases
    # check y coprime to base
    prime_y = [x for x in range(1, BASE) if gcd(x, BASE) == 1]
    print("Base:",BASE,"  Checking y:",prime_y)

    # Only check for primes greater than n in base n
    all_primes_greater_than = BASE
    all_the_primes = all_primes_up_to(PRIME_LIMIT)

    for a in all_the_primes:
        for b in all_the_primes:
            if a >= b:
                continue
            x = circle4primes(a, b)
            if x is None:
                print("Found counterexample: a:",a,"b:",b, flush=True)
                #exit(1)
            #else:
                #pprint((a, b, *x))
    '''
    #candidates = [[1061, 1117], [1229, 1291], [1559, 1621], [2039, 2131], [3299, 3361], [5557, 5573], [6329, 6451], [6337, 6473]]
    #candidates = [[2179, 31587561361]]
    #candidates = [[5, 7], [5, 11], [5, 13], [7, 11], [11, 13]]
    #candidates = [[11, 71]]
    #candidates = [[5, 13], [5, 37], [7, 11], [7, 23], [7, 59], [11, 19], [11, 31], [13, 29], [19, 23], [19, 47], [23, 43]]
    candidates = [[11, 19], [11, 43]]
    for ab in candidates:
        a = ab[0]
        b = ab[1]
        x = circle4primes(a, b)
        if x is None:
            print("Found counterexample: a:",a,"b:",b, flush=True)
        else:
            pprint((a, b, *x))
    '''