Ask Your Question
1

Speed a function with hex grid and primes numbers

asked 2018-11-30 17:31:22 +0200

wisher gravatar image

updated 2018-12-04 14:09:22 +0200

Hello, I am trying to solve the question 175 of the Turing Challenge. We create a network of hexagonal tiles containing the following integers 1, 2, 3, 4, ... I represent this network in a square matrix by placing the number 1 in the center of the matrix. We must calculate the difference between each tile numbered n and each of its six neighbors, we will call DP (n) the number of these differences # which are a prime number. My problem is this: I have to calculate a lot of times the function below. We save in DP3 the n which are surrounded by exactly 3 prime numbers. Here is all my script. My idea is this: every time I have completed a turn I can calculate the DP3 of the previous turn. The program stops when len (DP3) = 2016. As the spirals become bigger and bigger you have to make more and more calls to the DP3 () search function. My question is: is it possible to speed up this function?

PS: I'm sorry for my bad english. I also note that I am a beginner in programming.

# TURING 175 - Tuiles hexagonales (essai VII)  DERNIER 26/11/18   OK ma!s trop lent
# Une tuile hexagonale portant le numéro 1 est entourée d'un anneau de six tuiles numérotées de 2 à 7, la première tuile de l'anneau étant # posée à midi et en tournant dans le sens inverse des aiguilles d'une montre. 
# De nouveaux anneaux sont ajoutés de la même façon, les anneaux suivants étant numérotés 8 à 19, de 20 à 37, de 38 à 61, et ainsi de 
# suite. Le diagramme ci-contre montre les trois premiers anneaux. 

# En calculant la différence entre la tuile numérotée n et chacune de ses six voisines, nous appellerons DP(n) le nombre de ces différences # qui sont un nombre premier. 
# Par exemple, autour de la tuile numéro 8, les différences sont 12, 29, 11, 6, 1 et 13. Donc DP(8)=3. 
# De la même manière, autour de la tuile numéro 17, les différences sont 1, 17, 16, 1, 11 et 10, d'où DP(17)=2. 
# On peut montrer que la valeur maximale de DP(n) est 3. 
# Quand toutes les tuiles pour lesquelles DP(n)=3 sont énumérées dans l'ordre croissant pour former une suite, la 10ème tuile
# porte le numéro 271. 
# Trouver le numéro de la 2016ème tuile de cette suite.

temps = walltime()

var('n,x,y')

def rechercheDP3(n,x,y):
    n = M[x, y]
    d = [ M[x-1, y], M[x, y+1], M[x+1, y+1], M[x+1, y], M[x, y-1], M[x-1, y-1] ]    # les 6 cases qui entourent n
    dp = len( [t for t in d if is_prime(abs(n-t))] )
    if dp == 3:
        DP3.append(n)
    return

lim  = 17                  # Attention, il faut que lim soit impair
M = matrix(ZZ, lim, lim)     # création matrice nulle M    ATTENTION:  X est vertical vers le bas &  Y est horizontal vers la droite 
x, y = lim//2, lim//2        # centre de M on compte à partir de 0 
print "Centre:", x, y
M[x, y]     = 1              # le nombre 1 est placé au centre puis le 1er tour (2,3,4,5,6,7)
M[x-1, y]   = 2
M[x-1, y-1] = 3
M[x, y-1]   = 4 
M[x+1, y]   = 5
M[x+1, y+1] = 6
M[x, y+1]   = 7
n = 7
x_debut_prec, y_debut_prec = x-1, y

DP3 = [1]
maxtour = lim//2
tour = 1                    # on a déjà placé le 1er tour 2...7
grandeur_deplact = 1        # plus la spirale est extérieure plus le déplacement est grand  
print "maxtour = ", maxtour 
print "*",walltime(temps) 
#print M        # valable si n est petit
#print

# on remplit la table en spirale
while tour < maxtour and len(DP3)<2016:        # 2016 est la valeur demandée dans l'énoncé
    tour += 1
    grandeur_deplact +=1
    # on se place au-dessus de la case de départ (2,8,20,...) -> 6 trajets de 1 au 1er tour, de 2 au 2d tour, de 3 au 3e tour, ....
    x = x_debut_prec - 1
    n += 1                   #n = dernier_terme_tour + 1

    M[x, y] = n
    x_debut_suivant, y_debut_suivant = x, y     # coord case départ tour suivant

    for i in xrange(grandeur_deplact):     # trajet 1 horizontal gauche
        y = y-1
        n += 1
        M[x, y] = n

    for i in xrange(grandeur_deplact):     # trajet 2 vertical bas
        x = x+1
        n += 1
        M[x, y] = n

    for i in xrange(grandeur_deplact):     # trajet 3 oblique bas droit     
        x = x+1  
        y = y+1        
        n += 1
        M[x, y] = n

    for i in xrange(grandeur_deplact):     # trajet 4 horizontal droite       
        y = y+1
        n += 1
        M[x, y] = n        

    for i in xrange(grandeur_deplact):     # trajet 5 vertical haut        
        x = x-1        
        n += 1
        M[x, y] = n  

    for i in xrange(grandeur_deplact-1):     # trajet 6 oblique haut gauche -1 car on retombe sur le 1er -> le tour est fini
        x = x-1 
        y = y-1       
        n += 1
        M[x, y] = n          

    dernier_tour_n = n                        # on mémorise le tour qu'on vient de faire
    xfin_tour = x
    yfin_tour = y   


    # recherche les dp3 -> on recommence le tour précédent et on cherche ses DP3
    grandeur_deplact -= 1       
    # on se place au-début du tour précédent (2,8,20,38,62,92,...)
    x,y = x_debut_prec, y_debut_prec
    np = M[x,y]    # np comme cela on ne touche pas à n
    rechercheDP3(np,x,y)

    for i in xrange(grandeur_deplact):     # trajet 1 horizontal gauche
        y = y-1
        np += 1
        rechercheDP3(np,x,y)

    for i in xrange(grandeur_deplact):     # trajet 2 vertical bas
        x = x+1
        np += 1
        rechercheDP3(np,x,y)

    for i in xrange(grandeur_deplact):     # trajet 3 oblique bas droit     
        x = x+1  
        y = y+1        
        np += 1
        rechercheDP3(np,x,y)

    for i in xrange(grandeur_deplact):     # trajet 4 horizontal droite       
        y = y+1
        np += 1
        rechercheDP3(np,x,y)       

    for i in xrange(grandeur_deplact):     # trajet 5 vertical haut        
        x = x-1        
        np += 1
        rechercheDP3(np,x,y)

    for i in xrange(grandeur_deplact-1):     # trajet 6 oblique haut gauche
        x = x-1 
        y = y-1       
        np += 1
        rechercheDP3(np,x,y)

    # retour aux bonnes valeurs 
    grandeur_deplact += 1   # on revient à la valeur exacte pour le tour suivant
    x = x_debut_suivant  
    y = y_debut_suivant
    x_debut_prec, y_debut_prec = x_debut_suivant, y_debut_suivant

if len(DP3) >= 2016:
    print DP3(2015)
else:
    print "len(DP3) = ", len(DP3)

print "Durée = ", walltime(temps)
edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted
1

answered 2018-11-30 23:35:30 +0200

slelievre gravatar image

updated 2018-11-30 23:35:56 +0200

Using (n-t).abs().is_pseudoprime() instead of is_prime(abs(n-t)) is likely to result in a major speedup, at the cost of using only a strong pseudo-primality test instead of a provable primality test.

edit flag offensive delete link more

Comments

Thanks for your answer but your idea is slower. I have made a little try.

  • with my code: dp = len( [t for t in d if is_prime(abs(n-t))] ) I found an answer after 7.8 s

  • with your code: dp = len( [t for t in d if (n-t).abs().is_pseudoprime()] ) I found an answer after 10.4 s

wisher gravatar imagewisher ( 2018-12-04 00:10:08 +0200 )edit

It is really strange that it gets slower... Could you test something like abs(n-t).is_prime(proof=False) and see if it also takes 10 seconds?

Hilder Vitor Lima Pereira gravatar imageHilder Vitor Lima Pereira ( 2018-12-04 16:06:03 +0200 )edit

I tested your second idea in more depth. As soon as we increase the value of "lim" it is your second idea that is the most powerful.

lim         DP3      time in s
1001      78           4.6
2001      119        17.5
3001      152          38
4001      180          69
5001      208        108
6001      240        157
8001      288        282
16001    433       1151
30001    644       4077

The curve that adjusts the pairs of points (DP3, time) seems parabolic or exponential. So it's not pleasant since you have to reach a DP3 = 2016. Rather than the brute force of the computer, there must be a simpler mathematical process that goes beyond me. It is simply frustrating to feel incapable of not knowing the right simplification techniques. Thank you

wisher gravatar imagewisher ( 2018-12-05 10:47:09 +0200 )edit

I'm going to run the program for more value from "lim" quite to running my computer all night. I will provide the results later.

wisher gravatar imagewisher ( 2018-12-05 10:52:19 +0200 )edit

I did two more tests, here is the result:

lim       DP3           time in s
60001     1057          18 440
90001     1454          37 690  !
100001  -> the Sage software stops after a few seconds without an error.

I think the matrix has become too big. We must find another method of resolution, either mathematical or computer. I'll think about it. Anyway thank Mr. Lelièvre you for your help that allowed me to learn new concepts about prime numbers.

wisher gravatar imagewisher ( 2018-12-06 09:43:04 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-11-30 17:31:22 +0200

Seen: 449 times

Last updated: Dec 04 '18