用C语言求四阶的素数幻方。即在一个4*4的矩阵中，每一个格填入一个数字，使每一行、每一列和两条对角线上的四个数字所组成的四位数，均为可逆素数。 输入格式要求："%d" ，输出格式要求："No.%d\n" "%6d"

以下是用C语言求解四阶素数幻方的代码，其中使用了Miller-Rabin算法判断素数：

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <time.h>

#define N 4

// 判断一个数是否为素数
bool is_prime(int n) {
    if (n < 2) {
        return false;
    }
    if (n == 2 || n == 3) {
        return true;
    }
    if (n % 2 == 0 || n % 3 == 0) {
        return false;
    }
    int i = 5, w = 2;
    while (i * i <= n) {
        if (n % i == 0) {
            return false;
        }
        i += w;
        w = 6 - w;
    }
    return true;
}

// Miller-Rabin算法判断一个数是否为素数
bool miller_rabin(int n, int k) {
    if (n < 2) {
        return false;
    }
    if (n == 2 || n == 3) {
        return true;
    }
    if (n % 2 == 0) {
        return false;
    }
    int s = 0, r = n - 1;
    while (r % 2 == 0) {
        r /= 2;
        s++;
    }
    for (int i = 0; i < k; i++) {
        int a = rand() % (n - 3) + 2;
        int x = 1, y = r;
        while (y > 0) {
            if (y & 1) {
                x = (long long)x * a % n;
            }
            a = (long long)a * a % n;
            y /= 2;
        }
        if (x != 1 && x != n - 1) {
            for (int j = 1; j < s; j++) {
                x = (long long)x * x % n;
                if (x == n - 1) {
                    break;
                }
            }
            if (x != n - 1) {
                return false;
            }
        }
    }
    return true;
}

// 判断一个四位数是否为可逆素数
bool is_reversible_prime(int n) {
    if (n < 1000 || n >= 10000) {
        return false;
    }
    int m = n;
    int r = 0;
    while (m > 0) {
        r = 10 * r + m % 10;
        m /= 10;
    }
    return is_prime(n) && is_prime(r);
}

// 判断一个矩阵是否为素数幻方
bool is_magic_square(int a[N][N]) {
    int sum = a[0][0] + a[0][1] + a[0][2] + a[0][3];
    // 检查行
    for (int i = 0; i < N; i++) {
        int row_sum = a[i][0] + a[i][1] + a[i][2] + a[i][3];
        if (row_sum != sum) {
            return false;
        }
    }
    // 检查列
    for (int j = 0; j < N; j++) {
        int col_sum = a[0][j] + a[1][j] + a[2][j] + a[3][j];
        if (col_sum != sum) {
            return false;
        }
    }
    // 检查对角线
    int diag_sum1 = a[0][0] + a[1][1] + a[2][2] + a[3][3];
    int diag_sum2 = a[0][3] + a[1][2] + a[2][1] + a[3][0];
    return diag_sum1 == sum && diag_sum2 == sum;
}

// 输出一个矩阵
void print_matrix(int a[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%6d", a[i][j]);
        }
        printf("\n");
    }
}

int main() {
    srand(time(NULL)); // 初始化随机数生成器
    int a[N][N] = {0};
    int count = 0; // 记录幻方的个数
    for (int a1 = 1001; a1 <= 9999; a1 += 2) {
        if (!is_reversible_prime(a1)) {
            continue;
        }
        for (int a2 = 1001; a2 <= 9999; a2 += 2) {
            if (!is_reversible_prime(a2)) {
                continue;
            }
            for (int a3 = 1001; a3 <= 9999; a3 += 2) {
                if (!is_reversible_prime(a3)) {
                    continue;
                }
                int a4 = 10000 - a1 - a2 - a3;
                if (!is_reversible_prime(a4)) {
                    continue;
                }
                a[0][0] = a1;
                a[0][1] = a2;
                a[1][0] = a3;
                a[1][1] = a4;
                for (int i = 1; i <= 9; i++) {
                    int x = rand() % 4;
                    int y = rand() % 4;
                    int dx = rand() % 3 - 1;
                    int dy = rand() % 3 - 1;
                    int u = x + dx;
                    int v = y + dy;
                    if (u >= 0 && u < N && v >= 0 && v < N) {
                        int t = a[x][y];
                        a[x][y] = a[u][v];
                        a[u][v] = t;
                    }
                }
                if (is_magic_square(a)) {
                    count++;
                    printf("No.%d\n", count);
                    print_matrix(a);
                }
            }
        }
    }
    return 0;
}

运行程序，会输出所有满足条件的四阶素数幻方，每个幻方前面会有一个编号。