39개 나오는거 같네요.

간략한 풀이 과정입니다.

a1 a2 a3

a4 a5 a6

a7 a8 a9

의 각각의 합을 b 라고 하면

a1+a2+a3=b

a4+a5+a6=b

a7+a8+a9=b

f1) sum(a[1:9])=3*b 가 성립하고

a1+a5+a9=b

a2+a5+a8=b

a3+a5+a7=b

a4+a5+a6=b

f2) 3a5+sum(a[1:9])=4b

가 성립하게 됩니다.

f1)-f2) 를 하면

f3) a5=b/3 이게 되죠

이걸 가지고 마방진에 다시 입력을 해보면

b/3+d1 b/3+d2 b/3+d3

b/3+d4 b/3+0 b/3+d6

b/3+d7 b/3+d8 b/3+d9

이렇게 적어도 큰 문제는 없을거구요

d1+b/3 + 0+b/3 +d9+b/3 = b

d2+b/3 + 0+b/3 +d8+b/3 = b

d3+b/3 + 0+b/3 +d7+b/3 = b

d4+b/3 + 0+b/3 +d6+b/3 = b

이 성립하게되고

d9=-d1 d8=-d2 d7=-d3 d6=-d4 가 됩니다.

다시 적어보면

b/3+d1 b/3+d2 b/3+d3

b/3+d4 b/3+0 b/3-d4

b/3-d3 b/3-d2 b/3-d1

이 변화에는 큰 문제 없어보입니다. 이제 양측 끝쪽 가로 세로의 식을 적어보면

f3) b/3+d1 + b/3+d2 + b/3+d3 = b

f4) b/3+d3 + b/3-d4 + b/3-d1 = b 이고

d3=d4+d1 이 되고

-d2=d4+2*d1 가 되는것을 확인해볼수 있습니다.

즉 편차값 2개(d4,d1)가 있다면

이 두개를 더한값(d3=d4+d1)과

이 두개중 하나를 더 더한값 (-d2=d4+2*d1) 으로 편차값은 구성될수 있다고 볼수 있습니다.

그 내용을 python 3 로 표현하였습니다.

#마방진 테스트 함수 입니다. 
def is_mabang(test):
    _1st_line=sum(test[0:3])
    _2nd_line=sum(test[3:6])
    _3rd_line=sum(test[6:9])
    _1st_col=sum(test[0::3])
    _2nd_col=sum(test[1::3])
    _3rd_col=sum(test[2::3])
    RAO=sum(test[0::4])
    LAO=sum(test[2:7:2])
    result=[_1st_line,_2nd_line,_3rd_line,_1st_col,_2nd_col,_3rd_col,RAO,LAO]
    return len(set(result))==1

#솟수 생성함수입니다. 
def core03(l):
    if l<20:
        return [2, 3, 5, 7, 11, 13, 17, 19]
    SqrtL=int(l**0.5)
    TargetSet =set(range(SqrtL,l))
    TinyPrimeList=core03(int(l**0.50)+1)
    for i in TinyPrimeList:
        TargetSet = TargetSet -set(range(i*2,l,i))
    result= sorted(list(set(TinyPrimeList + list(sorted(TargetSet)))))
    return result
#500까지의 prime number 를 구합니다. 
total_prime_list=core03(500)
#print(len(total_prime_list))
#갯수는 95개군요. 
#이번에 알게된 재미난 툴 itertools 입니다. 
import itertools
'''
#뭔가.. 헛짓을 했더랬죠 이걸 실행시키면 hell gate 를 여는 느낌입니다. 
test01=itertools.permutations(total_prime_list,9)
k=0    
for j in test01:
    k+=1
    print(j)
    if k==200:
        break

test01=itertools.combinations(total_prime_list,9)
k=0    
for j in test01:
    k+=1
print(k)
'''
result01=[]
#중간 숫자를 기준으로 평균값이 중간이 나오는 숫자쌍을 찾아봅니다. 
for i in total_prime_list[4:-4]:
    temp=[]
    temp.append(i)
    if len(total_prime_list[total_prime_list.index(i)+1:])>=len(total_prime_list[:total_prime_list.index(i)]):
        for j in total_prime_list[:total_prime_list.index(i)]:
            if i-j+i in total_prime_list:
                temp.append(j)
                temp.append(i-j+i)
    else:
        for j in total_prime_list[total_prime_list.index(i)+1:]:
            if i-j+i in total_prime_list:
                temp.append(j)
                temp.append(i-j+i)
    #그러고는 그게 9개가 넘는지 확인하는 과정입니다. 
    if len(temp)>9:
        result01.append(sorted(temp))
#편차 a b c d(오름차순) 와 중간값 m 으로 9개의 숫자로 이루어진 마방진을 구성한다치면
#c==a+b 이며 d==a*2+b 또는 d==a+b*2일 경우 마방진이 성립하더군요.. 이외의 방법이 더 있을수 있지만 
#제가 가지고 있는 마방진의 지식은 이것이 전부였습니다.
#그래서 찾습니다. 꾸역꾸역   
#itertools.combinations 란 함수는 리스트 내부에서 주어진 숫자의 갯수만큼의 조합을 찾아내 주는...(말이 어렵습니다) 
#경우의 수 를 모두 찾아주지만 순서는 신경쓰지 않는... 아시죠? 중학교땐가 고등학교때인가.. ..

result02=[]
import itertools
for i in result01:
    temp01=[]
    temp02=sorted([i[(len(i)-1)//2]-x for x in i[:(len(i)-1)//2]])
    for j in itertools.combinations(temp02[:-2],2):
        if sum(j) in temp02:
            if sum(j)+j[0] in temp02:
                result02.append([i[(len(i)-1)//2],j[0],j[1],sum(j),sum(j)+j[0]])
            #한주영님의 지적으로 뒤늦게 추가. 또한번 더 수정!
            if sum(j)+j[1] in temp02:
                result02.append([i[(len(i)-1)//2],j[0],j[1],sum(j),sum(j)+j[1]])

#검산과정입니다 한주영님의 지적으로 버그 수정을 완료하였습니다. 
result03=[]
result04=[]
for i in result02:
    result03.append(sorted([i[0],i[0]-i[1],i[0]+i[1],i[0]-i[2],i[0]+i[2],i[0]-i[3],i[0]+i[3],i[0]-i[4],i[0]+i[4]]))
for i in result03:
    temp1=[]
    temp2=[]
    #d==a*2+b(a<b) 인 경우는 오름차수의 리스트가 다음과 같은 변형을 줄때 마방진이 형성됩니다. 
    for j in [1,6,5,8,4,0,3,2,7]:
        temp1.append(i[j])
    #d==b*2+a(a<b) 인 경우는 오름차수의 리스트가 다음과 같은 변형을 줄때 마방진이 형성됩니다. 
    for j in [6,5,1,0,4,8,7,3,2]:     
        temp2.append(i[j])
    if is_mabang(temp1):
        result04.append(temp1)
    elif is_mabang(temp2):
        result04.append(temp2)
    else:
        print(i)

'''
In [1]: result04
Out[1]:
[[17, 89, 71, 113, 59, 5, 47, 29, 101],         
 [41, 89, 83, 113, 71, 29, 59, 53, 101],        
 [103, 79, 37, 7, 73, 139, 109, 67, 43],        
 [29, 131, 107, 167, 89, 11, 71, 47, 149],      
 [139, 127, 43, 7, 103, 199, 163, 79, 67],      
 [37, 151, 139, 211, 109, 7, 79, 67, 181],      
 [157, 151, 73, 43, 127, 211, 181, 103, 97],    
 [43, 181, 157, 241, 127, 13, 97, 73, 211],     
 [173, 149, 71, 29, 131, 233, 191, 113, 89],    
 [47, 191, 173, 263, 137, 11, 101, 83, 227],    
 [199, 151, 67, 7, 139, 271, 211, 127, 79],     
 [59, 197, 191, 281, 149, 17, 107, 101, 239],   
 [61, 199, 193, 283, 151, 19, 109, 103, 241],   
 [37, 271, 163, 283, 157, 31, 151, 43, 277],    
 [71, 233, 197, 293, 167, 41, 137, 101, 263],   
 [53, 251, 197, 311, 167, 23, 137, 83, 281],    
 [233, 197, 89, 29, 173, 317, 257, 149, 113],   
 [257, 191, 89, 11, 179, 347, 269, 167, 101],   
 [251, 233, 89, 29, 191, 353, 293, 149, 131],   
 [71, 269, 233, 353, 191, 29, 149, 113, 311],   
 [73, 277, 229, 349, 193, 37, 157, 109, 313],   
 [43, 349, 241, 409, 211, 13, 181, 73, 379],    
 [137, 281, 263, 353, 227, 101, 191, 173, 317], 
 [317, 263, 101, 11, 227, 443, 353, 191, 137],  
 [89, 347, 281, 431, 239, 47, 197, 131, 389],   
 [73, 379, 271, 439, 241, 43, 211, 103, 409],   
 [331, 283, 109, 19, 241, 463, 373, 199, 151],  
 [61, 379, 283, 463, 241, 19, 199, 103, 421],   
 [71, 419, 263, 443, 251, 59, 239, 83, 431],    
 [41, 443, 269, 479, 251, 23, 233, 59, 461],    
 [47, 443, 281, 491, 257, 23, 233, 71, 467],    
 [173, 347, 269, 359, 263, 167, 257, 179, 353], 
 [137, 359, 293, 419, 263, 107, 233, 167, 389], 
 [59, 467, 281, 491, 269, 47, 257, 71, 479],    
 [149, 347, 311, 431, 269, 107, 227, 191, 389], 
 [347, 311, 149, 71, 269, 467, 389, 227, 191],  
 [359, 311, 137, 47, 269, 491, 401, 227, 179],  
 [103, 379, 331, 499, 271, 43, 211, 163, 439],  
 [101, 431, 311, 491, 281, 71, 251, 131, 461],  
 [167, 359, 353, 479, 293, 107, 233, 227, 419]] 
In [2]: len(result04)
Out[2]: 40

홀수 마방진의 경우 center 는 무조건 값 list 의 중간값이어야 합니다. 회전했을때 (즉, 90,180,270도를 돌리거나 vertical & horizontal flip) 을 시켰을때도 하나로 친다는 얘기는 한개의 중간값에 대한 경우 하나만 찾으면 나머지는 안찾아도 된다는 얘깁니다. 조건에 맞족하는 경우 하나만 찾으면 됩니다.

a = np.array().reshape(3,3) 마방진에서

d1 = abs(a[0][0] - a[1][1]) = abs(a[2][2] - a[1][1])

d2 = abs(a[0][1] - a[1][1]) = abs(a[2][1] - a[1][1])

d3 = abs(a[0][2] - a[1][1]) = abs(a[2][0] - a[1][1])

d4 = abs(a[1][0] - a[1][1]) = abs(a[1][2] - a[1][1])

입니다.

따라서 0 < d2 < d4 인 범위를 정할 경우 중간값 이상의 소수에 대하여 combinations (x,y) 을 구하면 d2 = x-mid, d4 = y-mid 이 경우 abs(d1) = (d4+d2)/2 이며 abs(d3) = (d4-d2)/2 이고

나머지 6숫자는, r = [mid-d1, mid+d3, mid-d4, mid-d3, mid-d2, mid+d1]

all(x in 소수리스트 for x in r) 이 True 인 경우만 찾으면 됩니다.

검증코드는 1x9 array 를 받아 3x3 으로 reshape 한 후 가로줄, 세로기둥, 대각선의 각 합들을 구하고 그것들이 모두 같은지 봅니다.

알고리즘 및 판정

from math import sqrt
from itertools import combinations

# 범위내 소수찾기
def findPrimes(m):
    prime_nums = []
    for n in range(2,m):
        find_factor = list((i, n//i) for i in range(1, int(sqrt(n))+1) if n%i == 0)
        if len(find_factor) == 1:
            prime_nums.append(find_factor[0][1])
    return prime_nums

# 리스트내에서 마방진찾기    
def findMgcSqr(p):
    result = []
    for i in range(4,len(p)-4):
        mid = p[i]
        for x,y in combinations(p[i+1:],2):
            d2, d4 = x - mid, y - mid
            d1, d3 = int((d4+d2)/2), int((d4-d2)/2)
            r = [mid-d1, mid+d3, mid-d4, mid-d3, mid-d2, mid+d1]
            if all(value in p for value in r):
                r = [r[0],x,r[1],y,mid] + r[2:]
                result.append(r)
    return result

P = findPrimes(500) # 500 이하의 소수, 95개
S = findMgcSqr(P) # 가능한 마방진, 중간값기준
print(*iter(S), sep = "\n") # 하나씩 출력
print("Available MagicSquares by prime numbers: ", len(S)) # 모두 몇개?

검증코드

import numpy as np


def isMagicSquare(s):
    square = np.array(s).reshape(3,3) # 1x9 -> 3x3 matrix
    rows = iter(square) # 한개씩 반복하능하게
    columns = zip(*square) # square 의 가로줄 각각에 대하여 iter 하여 값을 순서대로 zip 함.
    diagonals = zip(*[(square[i][i], square[-i-1][i]) for i in range(len(square))]) # 대각선방향으로 각각 두개씩 tuple 로 묶어서 그걸 다시 iter 로 순서대로 zip 하면 대각선 두개 묶음나옴

    def sums(*args): # iterable 받아서 각각 sum 을 리스트에 모으기
        nums = []
        for a in args:
            nums.append(sum(a))
        return set(nums)

    if sums(*rows) == sums(*columns) == sums(*diagonals): return True # 합계 리스트 내에 sum 들을 각각 비교할때 모두 같으면 True
    else: return False # 하나라도 틀리면 마방진 아님.


if all(isMagicSquare(X) for X in S): print("All MagicSquares")
else: print("No all MagicSquares")

결과

>>>
[17, 89, 71, 113, 59, 5, 47, 29, 101]
[41, 89, 83, 113, 71, 29, 59, 53, 101]
[37, 79, 103, 139, 73, 7, 43, 67, 109]
[29, 131, 107, 167, 89, 11, 71, 47, 149]
[43, 127, 139, 199, 103, 7, 67, 79, 163]
[37, 151, 139, 211, 109, 7, 79, 67, 181]
[73, 151, 157, 211, 127, 43, 97, 103, 181]
[43, 181, 157, 241, 127, 13, 97, 73, 211]
[71, 149, 173, 233, 131, 29, 89, 113, 191]
[47, 191, 173, 263, 137, 11, 101, 83, 227]
[67, 151, 199, 271, 139, 7, 79, 127, 211]
[59, 197, 191, 281, 149, 17, 107, 101, 239]
[61, 199, 193, 283, 151, 19, 109, 103, 241]
[37, 271, 163, 283, 157, 31, 151, 43, 277]
[71, 233, 197, 293, 167, 41, 137, 101, 263]
[53, 251, 197, 311, 167, 23, 137, 83, 281]
[89, 197, 233, 317, 173, 29, 113, 149, 257]
[89, 191, 257, 347, 179, 11, 101, 167, 269]
[89, 233, 251, 353, 191, 29, 131, 149, 293]
[71, 269, 233, 353, 191, 29, 149, 113, 311]
[73, 277, 229, 349, 193, 37, 157, 109, 313]
[43, 349, 241, 409, 211, 13, 181, 73, 379]
[101, 263, 317, 443, 227, 11, 137, 191, 353]
[137, 281, 263, 353, 227, 101, 191, 173, 317]
[89, 347, 281, 431, 239, 47, 197, 131, 389]
[109, 283, 331, 463, 241, 19, 151, 199, 373]
[73, 379, 271, 439, 241, 43, 211, 103, 409]
[61, 379, 283, 463, 241, 19, 199, 103, 421]
[71, 419, 263, 443, 251, 59, 239, 83, 431]
[41, 443, 269, 479, 251, 23, 233, 59, 461]
[47, 443, 281, 491, 257, 23, 233, 71, 467]
[173, 347, 269, 359, 263, 167, 257, 179, 353]
[137, 359, 293, 419, 263, 107, 233, 167, 389]
[149, 311, 347, 467, 269, 71, 191, 227, 389]
[137, 311, 359, 491, 269, 47, 179, 227, 401]
[149, 347, 311, 431, 269, 107, 227, 191, 389]
[59, 467, 281, 491, 269, 47, 257, 71, 479]
[103, 379, 331, 499, 271, 43, 211, 163, 439]
[101, 431, 311, 491, 281, 71, 251, 131, 461]
[167, 359, 353, 479, 293, 107, 233, 227, 419]
Available MagicSquares by prime numbers:  40
All MagicSquares

먼저 sieve와 primes를 초기화한 다음 찾아봅니다.

#define N 500

char sieve[N] = {0,};
int primes[N] = {0,};
int len = 0;

int main(int argc, const char * argv[]) {
    // sieve
    sieve[0] = -1;
    sieve[1] = -1;
    for (int i=2; i<N; i++) {
        if (sieve[i] == -1)
            continue;
        sieve[i] = 1;
        primes[len++] = i;
        for (int p=i*2; p<N; p += i) {
            sieve[p] = -1;
        }
    }

    int count = count_magic_squares(primes, len);
    printf("%d\n", count);

    return 0;
}

마방진 채우는 건 백트래킹을 사용하면 될것 같습니다.

// assume that nums[] are sorted
int count_magic_squares(int nums[], int len) {
    int count = 0;
    int square[9] = {0,};
    bt(square, 0, nums, len, &count);
    return count;
}

bt의 기본뼈대는 금방 만들어집니다. 한글자씩 채우면서 마방진 룰을 어기면 더이상 진행하지 않고 백트랙합니다.

void bt(int square[], int i, int nums[], int len, int* count) {
    if (i == 9) {
        ++(*count);
        return;
    }

    for (int k=0; k<len; k++) {
        if (nums[k]) { // not used
            square[i] = nums[k];
            nums[k] = 0; // use

            if (magic(square, i))
                bt(square, i+1, nums, len, count);

            nums[k] = square[i]; // recover
            square[i] = 0;
        }
    }
}

magic()은 채우진 숫자들로 마방진의 룰을 어기는지 확인합니다. (가로/세로/대각선 합이 같다) 여기서 조금이라도 빨리 룰 체크가 가능하도록 채우는 순서를 조금 바꿔봤습니다.

// 0 1 2
// 3 5 6
// 4 8 7

int magic(int s[], int i) {
    if (i < 3) return 1;
    int sum = s[0] + s[1] + s[2];

    if (i < 4) return 1;
    if (s[0] + s[3] + s[4] != sum) return 0;

    if (i < 5) return 1;
    if (s[2] + s[5] + s[4] != sum) return 0;

    if (i < 6) return 1;
    if (s[3] + s[5] + s[6] != sum) return 0;

    if (i < 7) return 1;
    if (s[0] + s[5] + s[7] != sum) return 0;
    if (s[2] + s[6] + s[7] != sum) return 0;

    if (i < 8) return 1;
    if (s[4] + s[8] + s[7] != sum) return 0;
    if (s[1] + s[5] + s[8] != sum) return 0;
    return 1;
}

이러면 .. 오래 걸려도 답은 나옵니다. 그런데.. 그 답은 90도씩 4번 회전 가능하고 거기다 거울상이 있으니 총 8번 중복되어 계산됩니다. 그냥 마지막에 8로 나눠줘도 되고 쉽게 기준 점을 정해서 0번이 2번과 4번보다 커야하고(회전 방지), 2번이 4번보다 커야 하고(거울상 방지)의 규칙을 추가하면 됩니다.

오래 걸리는 문제는.. 잘 모르겠네요. 일단 넣어보는 소수가 정렬되었다는 점을 이용해서 우선은 더 큰 값이 들어오는 걸 방지하는 것도 될 것 같아요. 그래서 룰을 어겼을 때 0을 반환하는 대신 더 큰 값은 시도하지 마라는 의미로 -1을 반환한다면 .. 조금 더 단축될 것 같습니다.

int magic(int s[], int i) {
    if (i < 2) return 1;
    if (s[0] < s[2]) return -1;    // 중복제거 및 s[2] 자리에 더 큰 값은 시도하지 마라 (-1)

    if (i < 3) return 1;
    int sum = s[0] + s[1] + s[2];  // 합을 구했음

    int s4 = sum - s[0] - s[3];    // s[4] 자리에 와야 하는 값이 정해짐
    if (s4 < 2) return -1;         // 만약 2보다 작다면 더 큰 값은 시도해볼 필요도 없음
    if (s4 > N || sieve[s4] != 1) return 0;  // 그 값이 소수 아니어도 더 볼 필요없음
    ...
}

이런 식으로 룰 체크를 조금 강화하니 조금 빨라지긴 하네요.

참, -1 반환 때문에 bt()도 바뀌었습니다.

void bt(int square[], int i, int nums[], int len, int* count) {
    if (i == 9) {
        //printsquare(square);
        ++(*count);
        return;
    }

    for (int k=0; k<len; k++) {
        if (nums[k]) { // not used
            square[i] = nums[k];
            nums[k] = 0; // use

            int result = magic(square, i);
            if (result == 1)
                bt(square, i+1, nums, len, count);

            nums[k] = square[i]; // recover
            square[i] = 0;

            if (result == -1) {  // no further search
                break;
            }
        }
    }
}

탐색 범위를 줄이기 위해 소수를 합이 같은 녀석들끼리 따로 모아놓고 찾습니다.

찾을 때는 중앙값과 대각선까지 5개 숫자를 정해 놓고 나머지는 계산으로 뽑아서 검사합니다.

중요한 부분은 다른 분들 풀이를 참고했습니다. 중앙값만 계산해도 대부분 필터링이 되네요.

답이 39개 나오는데 1개는 어디 갔을까요...

def primes_under(n):
    N = set(range(2, n))
    primes = []
    while N:
        m = min(N)
        primes.append(m)
        N -= set(range(m, n, m))
    return primes

# key: 두 소수의 합, value: 합이 key인 소수쌍들의 집합
def pairs_w_sum(primes):
    pairs = dict()
    for i in range(len(primes)):
        for j in range(i+1, len(primes)):
            p = primes[i]
            q = primes[j]            
            if p + q not in pairs: pairs[p + q] = set()
            pairs[p + q] ^= {p, q}

    return pairs

# 4쌍 이상이고, primes에 중앙값이 있는 집합만 정렬해서 리턴
def filter_pairs(pairs, primes):
    pairs2 = dict()
    for key in pairs:
        if len(pairs[key]) >= 8 and (key // 2) in primes:
            pairs2[key] = list(sorted(pairs[key]))
            #print(pairs2[key])

    return pairs2

def count(UB):
    primes = primes_under(UB)
    pairs = filter_pairs(pairs_w_sum(primes), primes)
    cnt = 0

    for key in pairs:
        lst = pairs[key]
        C = key // 2
        l = len(lst)
        for i in range(l // 2):
            UL, DR, = lst[i], lst[l - i - 1]    # \ 대각선
            for j in range(i + 1, l // 2 + 1):
                UR, DL = lst[j], lst[l - j - 1] # / 대각선
                sum = UL + C + DR

                # 외곽 수직, 수평 4줄
                U, D, L, R = sum-(UL+UR), sum-(DL+DR), sum-(UL+DL), sum-(UR+DR)
                used = {C, U, D, L, R, UL, UR, DL, DR}

                # 가운데 수직, 수평 2줄
                if (U in lst and D in lst and L in lst and R in lst) and \
                    U+C+D == L+C+R == sum and \
                    len(used) == 9:
                    cnt += 1
                    #print([UL, U, UR, L, C, R, DL, D, DR])
                    break

    return cnt

print(count(500)) #39

from itertools import combinations

def prime(n):
  ret = set(range(2, n+1))
  while 1:
    tmp = min(ret)
    if tmp > n**0.5: break
    else           : yield tmp
    ret -= set(range(tmp, n+1, tmp))
  yield from ret

def search_magic(v):
  ret = []
  for i in combinations(v, 3):
    c = min(i)
    t1, t2 = (j-c for j in i if j-c)
    if (t1+t2)%2 or (t1-t2)%2: continue
    t3, t4 = (t1+t2)/2, (t1-t2)/2
    tmp = [c+i for i in (t3, t4, -t1, -t2, -t3, -t4) if c+i in v]
    if len(tmp) == 6: yield tmp+[*i]

print(len([*search_magic([*prime(500)])]))  # 40