Magic Square Generator

magsq.c

/*
 * https://gist.github.com/bradenbest/af3c83230e11cf1b5aaeb05598398df6
 * magsq v1.0 by Braden Best, 2018
 * Public Domain, no warranty.
 *
 * A simple program that calculates a magic square of order n.
 *
 * The cool part is that I did all of this with a flat "1D" array using the
 * formula y * w + x to serialize a coordinate pair (x, y) as an array index,
 * which technically makes this a hash table.
 *
 * I wouldn't recommend trying to extend this unless you know what you are doing
 * This program is heavy with macros and tight operations that may not be
 * intuitive to a novice.
 *
 * This program works with pretty much any input, although there are things you
 * should know about certain landmarks.
 *
 * - Above 45
 *       the square takes up my entire screen and then some.
 * - Above 99
 *       the alignment starts getting thrown off by 5-digit numbers resulting in
 *       a jagged printout.
 * - In the 1000s
 *       the program starts to use MEGABYTES of memory. This is not bad memory
 *       management, it's just the fact that an int is 4 bytes and you're asking
 *       for a square of a million+ of them.
 * - Above 1625
 *       n(n^2+1)/2 (the magic sum for n) exceeds 2^31-1. Assuming ints are
 *       32-bit on your system (they usually are), this means we start getting
 *       our first negative sums.
 * - Above 3772
 *       the magic square checker breaks, as it starts reporting the
 *       magic squares as invalid. It's highly unlikely that all three magic
 *       square algorithms break at this point, especially since n^2 is still
 *       within the limits
 * - In the 10000s
 *       the program starts taking a really long time to execute. Again, this is
 *       related to just how big a number 10000^2 is and the fact that we're
 *       working with an algorithm that is O(n^2) at best.
 *       Memory usage: >100MB
 * - Above 46340
 *       n^2 surpasses 2^31-1, which means we start getting our first negative
 *       numbers in the *square*. Technically, the squares are still
 *       mathematically sound if you work in {(mod 2^64) - 2^32}, but
 *       information will be lost. Also, the output will be so big at this point
 *       that copying it down would probably be a lost cause.
 *       Memory usage: >8GB
 *
 *       At this point, changing all of the types to long unsigned will afford
 *       you more inputs as the limit will then become 2^63-1, but ask yourself:
 *       do you really need a magic square of order > 46340? Are you really
 *       willing to wait however long it takes to generate it?  Do you have the
 *       system resources to allocate f(x) = 8(x^2 + 1) bytes? Before you answer
 *       that, f(46340) is 17,179,164,808, or ~17GB. That's more memory than
 *       *my* system has. And don't forget about singly even squares, which use
 *       f(x) + f(x/2) bytes.  If you choose n=46342, it will cost 21GB.
 *
 * changelog:
 * 1.0
 *   ! release
 *   + implemented magsq_even (all inputs > 2 are now supported)
 * 0.5
 *   + added 4x4 party trick
 *   + added -p switch
 * 0.4
 *   + YXLOOP macro
 *   + SWAPA macro
 *   + implemented magsq_deven (inputs of 4, 8, 12, etc are now supported)
 * 0.3
 *   + added magsq_from_user()
 *   + added -t switch
 * 0.2
 *   + added MAGSQ_EXSUM macro
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

/* wraps x or y value in either direction to be between 0 and width */
#define WRAP(xory, width) \
    ( ((width) + (xory)) % (width) )

/* serializes 2D coordinate pair to 1D array index */
#define D2_1(x, y, width) \
    ( WRAP((y), (width)) * (width) + WRAP((x), (width)) )

/* wrapper for D2_1 using an array (coord pair) */
#define D2_1A(A, width) \
    D2_1( (A)[0], (A)[1], (width) )

/* adds the given x and y increments to an array (coord pair) */
#define ADDA(A, xi, yi) \
    ( ((A)[0] += (xi)), ((A)[1] += (yi)) )

/* generates a generic 2-dimensional Y/X loop header
 * [!] cannot be used as an expression
 */
#define YXLOOP(yi, yc, xi, xc) \
    for(int y = (yi); y < (yc); ++y) \
        for(int x = (xi); x < (xc); ++x)

/* swaps the array contents of index i and j */
#define SWAPA(A, i, j) \
    ( ((A)[(i)] = (A)[(i)] + (A)[(j)]), \
      ((A)[(j)] = (A)[(i)] - (A)[(j)]), \
      ((A)[(i)] = (A)[(i)] - (A)[(j)]) )

/* computes the expected sum for n (OEIS A006003) */
#define MAGSQ_EXSUM(n) \
    ( (n) * ((n) * (n) + 1) / 2 )

/* magsq_linesum wrapper macros */
#define MAGSQ_XSUM(ms, y) \
    magsq_linesum((ms), (y), magsq_xcoord)

#define MAGSQ_YSUM(ms, x) \
    magsq_linesum((ms), (x), magsq_ycoord)

#define MAGSQ_DSUM(ms) \
    magsq_linesum((ms), 0, magsq_dcoord)

#define MAGSQ_DISUM(ms) \
    magsq_linesum((ms), 0, magsq_dicoord)

/* generates an int(int,int,int) function for use with magsq_linesum
 * [!] cannot be used as an expression
 * [!] the first four arguments must be identifiers
 */
#define LAMBDA_GET_COORD(name, arg_seed, arg_var, arg_width, formula) \
    int name(int arg_seed, int arg_var, int arg_width){ return (formula); }

typedef int *magsq_t;
/* magsq_t: { width, [array of size width^2] } */

static struct {
    int use_lpmagsq;            // -l
    int use_report_magsq_check; // -c
    int test_mode;              // -t
    int party_trick;            // -p
} settings;

/* [mag]ic [sq]uare [struct]ure
 * n:   the order n
 * out: magic square
 *
 * generates and prepares a new magic square structure
 * [!] allocates from the heap
 */
magsq_t
magsq_struct(int n)
{
    size_t outsz = sizeof (int) * (1 + n * n);
    int *out = malloc(outsz);

    memset(out, 0, outsz);
    *out = n;
    return out;
}

LAMBDA_GET_COORD(magsq_xcoord, y, x, width, D2_1(x, y, width))
LAMBDA_GET_COORD(magsq_ycoord, x, y, width, D2_1(x, y, width))
LAMBDA_GET_COORD(magsq_dcoord, unused, d, width, D2_1(d, d, width))
LAMBDA_GET_COORD(magsq_dicoord, unused, d, width, D2_1((width-1)-d, d, width))

/* [mag]ic [sq]uare [linesum]
 * ms:        a magic square
 * seed:      the seed, used to select a specific row or column
 * get_coord: a function of type int(int,int,int) that computes a magic square
 *            index given the seed, the current iteration, and the width of the
 *            magic square. Applicable functions are defined with the macro
 *            LAMBDA_GET_COORD
 * out:       the sum of the given row, column or diagonal
 */
int
magsq_linesum(magsq_t ms, int seed, int (*get_coord)(int,int,int))
{
    int *map = ms + 1;
    int width = *ms;
    int total = 0;

    for(int var = 0; var < width; ++var)
        total += map[get_coord(seed, var, width)];

    return total;
}

/* [mag]ic [sq]uare [check]
 * ms:  a magic square
 * out: true (1) or false (0)
 *
 * checks if a magic square is valid.
 * This will be useful in implementing magsq_even and magsq_deven
 */
int
magsq_check(magsq_t ms)
{
    int *map = ms + 1;
    int width = *ms;
    int sum = MAGSQ_EXSUM(width);

    for(int y = 0; y < width; ++y)
        if(MAGSQ_XSUM(ms, y) != sum)
            return 0;

    for(int x = 0; x < width; ++x)
        if(MAGSQ_YSUM(ms, x) != sum)
            return 0;

    return MAGSQ_DSUM(ms) == sum && MAGSQ_DISUM(ms) == sum;
}

/* [mag]ic [sq]uare [odd]
 * n:   the order n
 * out: magic square
 *
 * generates an odd magic square
 */
magsq_t
magsq_odd(int n)
{
    magsq_t out = magsq_struct(n);
    int *map = out + 1;
    int coord[2] = {n / 2, 0};

    map[D2_1A(coord, n)] = 1;

    for(int i = 1; i < n*n; ++i){
        ADDA(coord, +1, -1);

        while(map[D2_1A(coord, n)] > 0 || (coord[0] == -1 && coord[1] == n))
            ADDA(coord, -1, +2);

        coord[0] = WRAP(coord[0], n);
        coord[1] = WRAP(coord[1], n);
        map[D2_1A(coord, n)] = i + 1;
    }

    return out;
}

/* [mag]ic [sq]uare [d]oubly [even]
 * n:   the order n
 * out: magic square
 *
 * generates a doubly even magic square
 */
magsq_t
magsq_deven(int n)
{
    magsq_t out = magsq_struct(n);
    int *map = out + 1;

    for(int i = 0; i < n*n; ++i)
        map[i] = i + 1;

    // swap top left and bottom right corners
    YXLOOP(0, n/4, 0, n/4)
        SWAPA(map, D2_1(x, y, n), D2_1((n-1)-x, (n-1)-y, n));

    // swap top right and bottom left corners
    YXLOOP(0, n/4, 3 * n/4, n)
        SWAPA(map, D2_1(x, y, n), D2_1((n-1)-x, (n-1)-y, n));

    // rotate center (has to be one-sided or else it will undo itself)
    YXLOOP(n/4, 2 * n/4, n/4, 3 * n/4)
        SWAPA(map, D2_1(x, y, n), D2_1((n-1)-x, (n-1)-y, n));

    return out;
}

/* [mag]ic [sq]uare [even]
 * n:   the order n
 * out: magic square
 *
 * generates a singly even magic square
 * algorithm based on this article: http://www.1728.org/magicsq3.htm
 * this is the most complicated, so there are comments explaining the logic
 * to help you, the reader, understand what the function is doing. "How x works"
 * comments are generally bad practice, so this is the only time I will do this.
 */
magsq_t
magsq_even(int n)
{
    magsq_t out = magsq_struct(n);
    magsq_t quadrant = magsq_odd(n / 2);
    int *map = out + 1;
    int const qn = n/2;
    int const qsq = qn * qn;
    int const xl = (qn + 1) / 2;
    int const xr = xl - 2;
    int const ym = xl - 1;

    /* Copy the mini magic square to all four quadrants
     * adding m(n/2)^2 where m = ...
     * 0 2
     * 3 1
     */
    YXLOOP(0, qn, 0, qn){
        map[D2_1(x, y, n)] = (quadrant + 1)[D2_1(x, y, qn)];
        map[D2_1(x + qn, y + qn, n)] = (quadrant + 1)[D2_1(x, y, qn)] + qsq;
        map[D2_1(x + qn, y, n)] = (quadrant + 1)[D2_1(x, y, qn)] + 2 * qsq;
        map[D2_1(x, y + qn, n)] = (quadrant + 1)[D2_1(x, y, qn)] + 3 * qsq;
    }

    free(quadrant);

    /* Swap the top left and bottom left quadrants in this pattern:
     * # # . . .
     * # # . . .
     * . # # . .
     * # # . . .
     * # # . . .
     */
    YXLOOP(0, qn, 0, xl)
        if((x == 0 && y != ym) || (x == xl - 1 && y == ym) || (x > 0 && x < xl - 1))
            SWAPA(map, D2_1(x, y, n), D2_1(x, y + qn, n));

    /* Swap the top right and bottom right quadrants in this pattern:
     * . . . . #
     * . . . . #
     * . . . . #
     * . . . . #
     * . . . . #
     */
    YXLOOP(0, qn, n - xr, n)
        SWAPA(map, D2_1(x, y, n), D2_1(x, y + qn, n));

    return out;
}

/* [mag]ic [sq]uare
 * n:   the order n
 * out: magic square
 *
 * generates a magic square of order n.
 */
magsq_t
magsq(int n)
{
    if(n == 2 || n < 1)
        printf("Warning: %i is outside of the domain of magsq(n).\n", n);

    if(n % 4 == 0)
        return magsq_deven(n);

    if(n % 2 == 0)
        return magsq_even(n);

    return magsq_odd(n);
}

/* [p]rint [mag]ic [sq]uare
 * ms: a magic square
 *
 * pretty-prints a magic square. Monospace font assumed
 * If n > 99, the spacing will be thrown out of alignment
 */
void
pmagsq(magsq_t ms)
{
    int width = *ms;
    int *map = ms + 1;

    for(int y = 0; y < width; ++y){
        for(int x = 0; x < width; ++x)
            printf("%4i  ", map[D2_1(x, y, width)]);

        puts("\n\n");
    }
}

/* [l]inear [p]rint [mag]ic [sq]uare
 * ms: a magic square
 *
 * prints a magic square "linearly"; i.e. in order.
 */
void
lpmagsq(magsq_t ms)
{
    int *map = ms + 1;
    int length = *ms * *ms;

    for(int i = 0; i < length; ++i)
        printf("%i%c ", map[i], i < length - 1 ? ',' : ' ');

    puts("");
}

/* [report] [mag]ic [sq]uare [check]
 * ms: magic square
 *
 * reports whether or not ms is a magic square
 */
void
report_magsq_check(magsq_t ms)
{
    static char *boolstr[] = {"is not", "is"};

    printf("Sum(%i): %i\n", *ms, MAGSQ_EXSUM(*ms));
    printf("This %s a magic square.\n", boolstr[magsq_check(ms)]);
}

/* get int
 * out: number
 *
 * reads stdin for a number, terminated by any character outside of [0-9]
 */
int
get_int(void)
{
    static char buf[1024];
    char *bufp = buf;
    int ch;

    while((ch = getchar()) != EOF && isdigit(ch))
        *bufp++ = ch;

    *bufp++ = '\0';

    if(ch == EOF){
        puts("magsq: reached end of input. abort.");
        exit(1);
    }

    return atoi(buf);
}

/* [get] [pos]itive [int]
 * out: number
 *
 * calls get_int and enforces that n > 0
 */
int
get_pos_int(void)
{
    int n;

    while((n = get_int()) <= 0)
        ;

    return n;
}

/* [mag]ic [sq]uare [from user]
 * out: magic square
 *
 * reads a magic square from stdin
 */
magsq_t
magsq_from_user(void)
{
    magsq_t ms;
    int *map;

    printf("Enter a value for n >");
    ms = magsq_struct(get_pos_int());
    map = ms + 1;
    printf("Enter the magic square in sequence (e.g. 2 7 6 9 5 1 4 3 8)\n");
    printf("Expected: %i numbers\n>", *ms * *ms);

    while(map - ms < *ms * *ms + 1)
        *map++ = get_pos_int();

    return ms;
}

/* party trick
 * out: magic square
 *
 * reads a number N from stdin, and generates a 4x4 magic square equal to N
 * It's only impressive the first couple times.
 */
magsq_t
party_trick(void)
{
    magsq_t out = magsq_struct(4);
    static int preset[16] = {
         0,  1, 12, 7,
        11,  8,  0, 2,
         5, 10,  3, 0,
         4,  0,  6, 9
    };
    int *map = out + 1;
    int n;

    memcpy(map, preset, sizeof (int) * 16);
    printf("Enter a number between 22 and 99 >");

    while((n = get_pos_int()) < 22 || n > 99)
        ;

    map[0] = n - 20;
    map[13] = n - 19;
    map[11] = n - 18;
    map[6] = n - 21;

    return out;
}

void
usage(void)
{
    puts("magsq - generates a magic square");
    puts("usage: magsq [options] <arguments>");
    puts("arguments:");
    puts("  n   (required) the order (width) of the magic square to be generated");
    puts("options:");
    puts("  -c  [c]heck if magic square is a magic square");
    puts("  -l  print magic square's contents [l]inearly");
    puts("  -t  [t]est a user-provided magic square (implies -c)");
    puts("  -p  [p]arty trick");
    exit(1);
}

int
check_switch(char *arg)
{
    struct {
        char *flag;
        int *cfg;
        int end;
    } switches[] = {
        {"-c", &(settings.use_report_magsq_check)},
        {"-l", &(settings.use_lpmagsq)},
        {"-t", &(settings.test_mode)},
        {"-p", &(settings.party_trick)},
        {.end = 1}
    };

    for(int i = 0; !switches[i].end; ++i)
        if(strncmp(arg, switches[i].flag, strlen(switches[i].flag)) == 0)
            return *(switches[i].cfg) = 1;

    return 0;
}

magsq_t
parse_args(char **args)
{
    char *arg_n = NULL;

    while(*++args)
        if(!check_switch(*args))
            arg_n = *args;

    if(settings.party_trick)
        return party_trick();

    if(settings.test_mode){
        settings.use_report_magsq_check = 1;
        return magsq_from_user();
    }

    if(arg_n == NULL)
        usage();

    return magsq(atoi(arg_n));
}

int
main(int argc, char **argv)
{
    magsq_t ms = parse_args(argv);

    pmagsq(ms);

    if(settings.use_lpmagsq)
        lpmagsq(ms);

    if(settings.use_report_magsq_check)
        report_magsq_check(ms);

    free(ms);
    return 0;
}

If you think this is complicated, you should see it after the preprocessor unwraps the macros.

before:

// swap top left and bottom right corners
YXLOOP(0, n/4, 0, n/4)
    SWAPA(map, D2_1(x, y, n), D2_1((n-1)-x, (n-1)-y, n));

after:

for(int y = (0); y < (n/4); ++y) for(int x = (0); x < (n/4); ++x)
    ( ((map)[(( ( (((n)) + ((y))) % ((n)) ) * (n) + ( (((n)) + ((x))) % ((n)) ) ))] = (map)[(( ( (((n)) + ((y))) % ((n)) ) * (n) +
        ( (((n)) + ((x))) % ((n)) ) ))] + (map)[(( ( (((n)) + (((n-1)-y))) % ((n)) ) * (n) + ( (((n)) + (((n-1)-x))) % ((n)) ) ))]),
        ((map)[(( ( (((n)) + (((n-1)-y))) % ((n)) ) * (n) + ( (((n)) + (((n-1)-x))) % ((n)) ) ))] = (map)[(( ( (((n)) + ((y))) % ((n)) ) *
        (n) + ( (((n)) + ((x))) % ((n)) ) ))] - (map)[(( ( (((n)) + (((n-1)-y))) % ((n)) ) * (n) + ( (((n)) + (((n-1)-x))) % ((n)) ) ))]),
        ((map)[(( ( (((n)) + ((y))) % ((n)) ) * (n) + ( (((n)) + ((x))) % ((n)) ) ))] = (map)[(( ( (((n)) + ((y))) % ((n)) ) * (n) + ( (((n)) +
        ((x))) % ((n)) ) ))] - (map)[(( ( (((n)) + (((n-1)-y))) % ((n)) ) * (n) + ( (((n)) + (((n-1)-x))) % ((n)) ) ))]) );

macros:

#define YXLOOP(yi, yc, xi, xc) \
    for(int y = (yi); y < (yc); ++y) \
        for(int x = (xi); x < (xc); ++x)

#define SWAPA(A, i, j) \
    ( ((A)[(i)] = (A)[(i)] + (A)[(j)]), \
      ((A)[(j)] = (A)[(i)] - (A)[(j)]), \
      ((A)[(i)] = (A)[(i)] - (A)[(j)]) )

#define D2_1(x, y, width) \
    ( WRAP((y), (width)) * (width) + WRAP((x), (width)) )

#define WRAP(xory, width) \
    ( ((width) + (xory)) % (width) )

They really reduce complexity. Any C programmer worth their salt ought to be proficient with them.