Copyright (C) 2017 Ben Asselstine

 A set of tools for finding and managing 3x3 magic squares of squares.

Finding a 3x3 square of squares is an ongoing exploration for amateur mathemeticians. As of the time of this writing, nobody has found one with 8 perfect squares. Maybe you can find it! For more information, see an introduction to the subject, the pages on MultiMagie, and the Numberphile video that inspired me to write this software.

This project is for generating and managing 3x3 magic squares that contain perfect squares. It uses libgmp (GNU Multiprecision Arithmetic Library) for arbitrary-precision arithmetic on very large integers. Fituvalu is written to be easily parallelizable with GNU Parallel.

What does the name mean?
The word Fituvalu is a compound word of the Samoan words for "seven" (fitu) and "eight" (valu).

There are top level programs that we want to run primarily, and then there are helper programs. These programs are designed in the UNIX fashion; they are command-line programs with the output of one program being passed to the input of another via pipes.

Note: Similarly named-programs with the small- prefix means they use 64-bit integers, and not arbitrary-precision integer types. They are much faster (but limited), and are included to verify the correctness of their arbitrary-precision counterparts.

There are three approaches used to create 3x3 magic squares.

There are two kinds of progressions that we are searching for:
The "step" progression, and the "fulcrum" progression.

Step Progression:

               D1       E1

The step progression is defined by two distances, D1 and E1.

Fulcrum Progression:

              D1    E1   F1

The fulcrum progression is defined by three distances, where the D1 distance adds up to the sum of E1 and F1.

These progressions can be trivially transformed into 3x3 magic squares.

Common Options


When using arbitrary-precision integers it is faster to deal with their raw formats. Many of these programs have an option for reading or writing the raw GMP formats rather than textual representations of numbers. When reading and writing billions of records, the time savings can be considerable.


Many of these programs iterate over perfect squares. Instead of iterating to the next square, we can go NUM squares forward using this option.

Top-level Programs

step-progression and fulcrum-progression programs iterate through squares to find 3x3 magic squares that have at least 5 perfect squares.

program: step-progression-1
Usage: step-progression-1 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX

Here are the perfect squares we're searching for in this program:

              ^  ^  ^       ^                           
              1  2  3       4                           

All of the progression programs will provide a visual depiction of the squares they're searching for in their --help. The --squares option causes the program to take the starting perfect square from the standard input so you can have more control over what progressions get generated.

You run it like this:

        $ step-progression-1 1 1000 | head -n1
        1, 25, 49, 121, 145, 169, 241, 265, 289
program: step-progression-2
Usage: step-progression-2 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
                 ^  ^       ^  ^                        
                 1  2       3  4                        
program: step-progression-3
Usage: step-progression-3 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
               ^     ^       ^     ^
               1     2       3     4
program: reverse-step-progression-1
Usage: reverse-step-progression-1 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
                                  ^       ^  ^  ^       
                                  4       3  2  1       
program: fulcrum-progression-1
Usage: fulcrum-progression-1 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
          ^        ^  ^     ^                           
          1        2  3     4                           
program: fulcrum-progression-2
Usage: fulcrum-progression-2 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
                   ^  ^     ^  ^                        
                   1  2     3  4                        
program: fulcrum-progression-3
Usage: fulcrum-progression-3 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
                            ^  ^  ^     ^               
                            1  2  3     4               
program: fulcrum-progression-4
Usage: fulcrum-progression-4 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
           ^        ^        ^     ^
           1        2        3     4
program: fulcrum-progression-5
Usage: fulcrum-progression-5 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
                    ^ ^       ^   ^
                    1 2       3   4
program: reverse-fulcrum-progression-1
Usage: reverse-fulcrum-progression-1 [--squares=FILE] [--in-binary] [--out-binary] [--increment=NUM] MIN MAX
                      ^     ^  ^  ^
                      4     3  2  1
program: morgenstern-search
Usage: morgenstern-search [--filter=NUM] [--out-binary] [--mem] -12345678 FILE

Whereas the prior programs generate 3x3 magic squares with 4 perfect squares and then hope that a 5th, 6th, or 7th perfect square shows up, the morgenstern-search programs generate 3x3 magic squares with 5 perfect squares and then hope that a 6th or 7th perfect square shows up.

The approach is create two arithmetic progressions of 3 perfect squares with a single square in common. The results are much larger, and less varied in terms of magic square types.

The standard input provides the parametric MN values -- two values per record to assist in the transformation. The seq-morgenstern-mn program provides these numbers.

The --filter=NUM option displays magic squares that have NUM or more perfect squares.

The --mem option reads FILE into memory and treats the standard input differently. The idea here is that the stdin contains the row numbers that we're going to check. Simply put, the row number an index into FILE which is M,N, and then we get P,R as every record from 0 to the record number. This option avoids the generation of duplicate magic squares, and saves a lot of time.

When -i used with --mem, the file that is passed on the standard input needs to be encoded with convert-binary-gmp-records-to-text -i --ull.

To generate Morgenstern type 1 squares, use the -1 option. For type 2, use the -2 option and so on. To see what the configuration looks like you can use the --show=TYPE option.

Morgenstern type 1 squares have 5 perfect squares in this configuration:

        |  A^2  |       |  C^2  |
        |       |  E^2  |       |
        |  G^2  |       |  I^2  |

In this case, A^2, E^2, I^2 is one arithmetic progression of three perfect squares, and C^2, E^2, G^2 is the other. They share E^2 in common.

Morgenstern type 2 squares have 5 perfect squares in this configuration:

        |       |  B^2  |       |
        |  D^2  |  E^2  |  F^2  |
        |       |  H^2  |       |

Morgenstern type 3 squares have 5 perfect squares in this configuration:

        |  A^2  |  B^2  |       |
        |       |  E^2  |       |
        |       |  H^2  |  I^2  |

Morgenstern type 4 squares have 5 perfect squares in this configuration:

        |       |  B^2  |       |
        |  D^2  |       |  F^2  |
        |  G^2  |       |  I^2  |

Morgenstern type 5 squares have 5 perfect squares in this configuration:

        |  A^2  |       |  C^2  |
        |  D^2  |       |  F^2  |
        |       |  H^2  |       |

Morgenstern type 6 squares have 5 perfect squares in this configuration:

        |  A^2  |  B^2  |       |
        |       |       |  F^2  |
        |  G^2  |  H^2  |       |

Morgenstern type 7 squares have 5 perfect squares in this configuration:

        |       |  B^2  |  C^2  |
        |  D^2  |       |       |
        |       |  H^2  |  I^2  |

Morgenstern type 8 squares have 5 perfect squares in this configuration:

        |       |  B^2  |  C^2  |
        |       |  E^2  |       |
        |  G^2  |  H^2  |       |
program: 3sq-pair-search
Usage: 3sq-pair-search [--in-binary] [--out-binary] N1, N2, N3,
  or:  3sq-pair-search [OPTION...] FILE

The idea here is to make 3x3 magic squares from two arithmetic progressions of three perfect squares. The progressions require a single number in common; either the first number (the lowest number) or the center number. When they don't have a number in common, this program tries to scale the progressions forward so they do.

These progressions are generated by running find-3sq-progressions, find-3sq-progressions-mn, find-3sq-progressions-gutierrez, find-3sq-progressions-dist and mine-3sq-progressions.

This program accepts a progression on the standard input, and check it against a set of progressions by passing in FILE. Alternatively we can check it against a single progression by passing in N1..N3 as arguments on the command-line.

program: scissor-square
Usage: scissor-square [--in-binary] [--out-binary] N1, N2, N3,
  or:  scissor-square [--in-binary] [--out-binary] FILE

Try to create a 3x3 of magic square from two progressions of three squares that have the same difference between the squares, or the differences must be divisible. This program tries create a magic square with the following layouts:

        +------+------+------+     +------+------+------+
        |      |  X3  |  N2  |     |  X2  |      |  N3  |
        +------+------+------+     +------+------+------+
        |  N3  |      |  X1  | and |      |  N2  |  X1  |
        +------+------+------+     +------+------+------+
        |  X2  |  N1  |      |     |  N1  |  X3  |      |
        +------+------+------+     +------+------+------+

The first progression is given as the arguments N1, N2, N3. The second progression is passed in on the standard input. Alternatively you can try the progressions on the standard input with a set of progressions by passing in FILE.

The three values must be perfect squares, separated by a comma and terminated by a newline, and must be in ascending order. This program also generates 3x3 nearly-magic squares.

Helper Programs

program: permute-square
Usage: permute-square [--display-offset] [--show-all]

Try all permutations of a given progression on the standard input until it becomes a magic square. With the --display-offset option, this program generates transformations for use in the transform-square program. With the --show-all option, this command will show all possible permutations of the square.

program: auto-transform-square
Usage: auto-transform-square [-n]

Like permute-square but much faster because it uses transform-square. It takes the transformations in the first 100 squares and uses them to transform the rest. The 100 can be changed with the -n option.

program: transform-square
Usage: transform-square [OPTION...] TRANSFORM-FILE

Transform a set of progressions into magic squares given a data file of transformations. When a progression can't be transformed it is displayed on the standard error.

program: count-squares
Usage: count-squares [--filter=NUM-SQUARES]

Count how many perfect squares a magic square has, or filter squares that only have a given number of squares using the --filter option.

program: reduce-square
Usage: reduce-square [OPTION...]

Determine what the lowest value is for the given magic square and have it still retain all of the perfect squares.

program: rotate-square
Usage: rotate-square [OPTION...]

Given a magic square, show the rotations and reflections of it.

program: uniq-squares
Usage: uniq-squares [OPTION...]

Given a set of magic squares, only show the unique ones. E.g. take into account all rotations, reflections and reductions.

program: type-square
Usage: type-square [--filter=NUM_SQUARES:TYPE] [---no-show-squares]

Given a magic square determine which of the configurations it is. Using the --filter option, this program filters a given set of magic squares such that they only have the given type. It displays records of the form: 5:12 (e.g. this magic square has 5 perfect squares, laid out in type 12.)

program: sort-square
Usage: sort-square [--prefix=NAME]

Given a magic square, put it in a file called fives, sixes, sevens etc. according to how many perfect squares it contains.

program: odd-square
Usage: odd-square [--even] [--mixed] [--doubly-even]

Given a magic square, show it if it's comprised of odd numbers. The --even option shows even instead of odd, --doubly-even shows the squares where all the numbers are divisible by 4, and --mixed has both odd numbers and even numbers.

program: number-square
Usage: number-square [--min] [--max] [--filter]

Given a 3x3 magic square on the standard input, output the magic number. Only do so if the square is a magic square. Use --min and --max to display the smallest and largest values in the square. Use the --filter=NUM option to show only squares that have NUM as a magic number.

program: check-square
Usage: check-square [OPTION...]

Given a set of 9 numbers on the standard input, check if it's a magic square. Only display it if it is.

program: lucas-square
Usage: lucas-square [--show-abc] [--in-binary] [--out-binary] [--inverse] [--from-abc]

Given a 3x3 magic square on the standard input, only show it if it is in the Lucas family of squares. The Lucas family has the following structure:

        |    c - b    | c + (a + b) |    c - a    |
        | c - (a - b) |      c      | c + (a - b) |
        |    c + a    | c - (a + b) |    c + b    |

The --inverse option only shows squares that aren't in the Lucas family. The --show-abc option displays the a, b, and c values that generate the magic square. The --from-abc takes in the a, b, and c values from the standard input, and generates a 3x3 magic square.

program: progression-type-square
Usage: progression-type-square [--filter=TYPE] [--show-def]

Given a magic square on the standard input, determine which of the two progression types it is. 1 is to "step progression", and 2 is "fulcrum progression". Using the --filter option, this program filters magic squares based on progression type. The --show-def option shows the distances between the numbers. For "step progression" squares, it shows 2 distances, and for "fulcrum progression" squares" it shows 3 distances.

program: display-square
Usage: display-square [--number-line] [--no-magic-number] [--show-sums] [--simple]

Pretty-print a 3x3 square and put the values inside the cells. This program will also display the sums of the rows and columns and diagonals when the square is not magic.

program: sq-seq [--out-binary]
Usage: sq-seq [OPTION...] FIRST
  or:  sq-seq [OPTION...] FIRST LAST

Generate a sequence of perfect squares. This program works like seq from GNU Coreutils, but loops over squares. It can be given an increment to loop over every Nth square in the sequence.

program: 4sq
Usage: 4sq [--filter=NUM] [--type=NAME] [--in-binary] [--out-binary] [--increment=NUM] START FINISH
  or:  4sq [--filter=NUM] [--type=NAME] [--in-binary] [--out-binary] [--increment=NUM] FINISH

Generate a progression of 4 perfect squares. 4 squares in an arithmetic progression is impossible (e.g. where the distance between them is the same), but if the distance between one or more of the squares is varied, it can still form the basis of a 3x3 magic square. 4sq has many different strategies for generating 4 square progressions which is controlled by the --type option. It can iterate over squares between a START and a FINISH, or it can take the starting squares from the standard input.

program: complete-4sq-progression
Usage: complete-4sq-progression [--filter=NUM] [--type=NAME] [--in-binary] [--out-binary]

Generate a nine number progression based on four squares from the standard input, based on a progression strategy defined by the --type option. The --filter=NUM option only shows progressions with NUM squares or more. The idea is that it accepts input from the 4sq program.

program: gen-progression
Usage: gen-progression [--filter=NUM] [--type=NAME] [--in-binary] [--out-binary] [--increment=NUM] START FINISH
  or:  gen-progression [--filter=NUM] [--type=NAME] [--in-binary] [--out-binary] [--increment=NUM] FINISH

Generate a progression of 9 numbers that have at least 4 perfect squares. It can iterate over squares between a START and a FINISH, or it can take the starting squares from the standard input. This progression can be easily transformed into a 3x3 magic square using the auto-transform-square command. The --filter=NUM option only shows progressions with NUM squares or more.

This program is equivalent to using 4sq and piping it into complete-4sq-progression. It's faster because it's not importing/exporting numbers to and from strings, and it doesn't involve a pipe.

The following are equivalent:

        $ step-progression-1 1 50000  >/tmp/dump1
        $ gen_progression --type=step-progression-1 --filter 5 1 50000  >/tmp/dump2
        $ diff /tmp/dump1 /tmp/dump2
program: seq-morgenstern-mn
Usage: seq-morgenstern-mn MAX

Generates a sequence of numbers that has the form: M > N > 0, where M and N are coprime, and with one being even and the other one odd. The output of this program is suitable for input into the morgenstern-search- programs, as well as the find-3sq-progressions program.

program: mine-3sq-progressions
Usage: mine-3sq-progressions [--show-diff] [--in-binary] [--show-root] [--show-sum]

Accept 3x3 magic squares on the standard input and output any arithmetic progressions that are formed by three perfect squares. The idea here is that we are mining the results of our previously generated magic squares for progressions that can be used to make new magic squares. This program generates useful input for 3sq-pair-search.

program: find-3sq-progressions
Usage: find-3sq-progressions [--in-binary] [--increment=NUM] [--no-root] [--out-binary] NUM

Given a perfect square on the standard input, generate three equidistant perfect squares. NUM is how many times to advance the second square forward, take the diff and check if the third number is a square or not. sq-seq generates squares for this program to use. This program generates useful input for complete-3sq-progression, and 3sq-pair-search.

program: find-3sq-progressions-mn
Usage: find-3sq-progressions-mn [--in-binary] [--no-root] [--out-binary]

Given a pair of MN parametric values on the standard input, generate three equidistant perfect squares. It generates records of four numbers separated by commas, with the final number being the square root of the third square. This program generates useful input for complete-3sq-progression, and 3sq-pair-search.

program: find-3sq-progressions-gutierrez
Usage: find-3sq-progressions-gutierrez [--show-diff] [--in-binary] [--no-root] [--out-binary] MAX [C [E]]

An implementation of an algorithm described by Eddie N. Gutierrez on www.oddwheel.com, for generating progressions of perfect squares. These progressions are suitable for the diagonal line of a 3x3 magic square, as seen in Andrew Bremner's magic square of 7 squares. It seems that B and C are always prime numbers, and E is always even, but that's all I know. MAX is how many times the program will iterate its sequence of numbers forward.

This is one of the few programs in this suite of tools that takes negative numbers as arguments on the command-line. You will need to put -- before a negative number, to signal to the program that your negative numbers are not command-line options. For example:

$ find-3sq-progressions-gutierrez -n -- 4 1 -3
4, 4, 4, 
1, 1, 1, 
4, 4, 4, 
49, 25, 1, 
Here are the values I've tested:
  MAX   B   C  E
  --- --- --- --
  100   1  -3  2
  100   1 -13 14
  100   1 -17  6
  100   1 -19 10
  100   1 -79 20
  100   1   3  2
  100   1  11 10
  100   1  17  4
  100   1  19  6
  100   1  57 28
  100  -3   1  4
  100  -7   1  8
  100 -15   1  8
  100 -17   1 18
  100   3   1  2
  100   9   1  4
  100  33   1  8
  100  65   1 16

program: find-3sq-progressions-dist
Usage: find-3sq-progressions-dist [--show-diff] [--in-binary] [--no-root] [--out-binary] [--show-sum] [--start=NUM] [--tries=NUM] [--range=NUM] [--no-second] [--no-third] [DISTANCE]

Generate arithmetic progressions of three squares, that are DISTANCE apart. If DISTANCE is not provided as an argument, it is read from the standard input. By default this program checks 5 million squares for a progression with the given distance. It searches every perfect square starting at 1, for two other perfect squares DISTANCE apart.

The --range=NUM options stops the program from checking when we get NUM away from the starting point.

The --no-second and --no-third options generate progressions where the 2nd or 3rd numbers are not perfect squares.

program: complete-3sq-progression
Usage: complete-3sq-progression NUM [--increment=NUM]

Given an arithmetic progression of three perfect squares and the root of the third square on the standard input, generate a progression of four squares. The fourth square iterated forward NUM times. This program generates useful input for complete-4sq-progression.

program: convert-binary-gmp-records-to-text
Usage: convert-binary-gmp-records-to-text [--num-columns=NUM] [--inverse] [--ull]

The progression programs have an option to dump records in raw format (--out-binary), so we don't spend time converting the numbers to their textual representation. This program converts them from the GMP raw format to that textual representation. It reads from the standard-input and outputs the records on the standard output. The --inverse option converts from text to binary. The --ull option converts unsigned long longs instead of GMP numbers.

How do you use all of this?

You do something like:

$ step-progression-1 1 36893488147419103232 > data &

(The big number is 2^65. It's okay that it's not a perfect square.)

Then while it's running you can check it by doing this:

$ cat data | auto-transform-square | type-square | sort | uniq -c
   1565 5:12
   1667 5:21
      4 5:8
      6 5:9
      2 6:14

(e.g. Take the progressions, turn them into magic squares, show the type of square along with the number of squares it has, sort them, and then count them all up.)

We can see there are two six-square magic squares with configuration type 14. But are they the same, or is it a duplicate?

$ cat data | auto-transform-square | count-squares -f 6 | uniq-squares | type-square | sort | uniq -c
      1 6:14

Ah we can see there's only one unique six found.

(e.g. Take the progressions, turn them into magic squares, filter out all but the six square magic squares, throw away the duplicate squares, show the type of square along with the number of squares it has, sort them, and then count them all up.)

Did you find a seven? You can look at it like this:

$ cat data | auto-transform-square | count-squares -f7 | uniq-squares | display-square 
|  373^2 |  289^2 | 565^2  |
| 360721 |  425^2 |  23^2  |    magic number = 541875
|  205^2 | 527^2  | 222121 |

Do you want to take a look at a particular types of square?

$ cat data | auto-transform-square | type-square -f 5:8 -f 5:9 | uniq-squares | wc -l

(e.g. Pull the magic squares that have 5 perfect squares, and are in configurations 8 and 9. Get rid of the duplicates and count up the uniques ones.)

How about those morgenstern programs?

$ seq-morgenstern-mn 100 | morgenstern-search-type-8 100  -f 6 | check-square >data

The morgenstern-search- programs do not require any transformations -- they do however, need to be checked for correctness because the algorithm generates squares with non-distinct values.

How can I mine the magic squares I've already found for other magic squares?

$ cat mysquares | mine-3sq-progressions | uniq > data
$ cat data | sort -g | uniq > data.sorted
$ cat data.sorted | 3sq-pair-search data.sorted > moresquares

Instead of using mine-3sq-progressions, you can also use:

$ seq-morgenstern-mn 10000 | find-3sq-progressions-mn | uniq > data


$ sq-seq 1000000000 | find-3sq-progressions 100000 | uniq > data

Running out of cpu cycles? Parallelize it with GNU Parallel.

If you have a 32 core system you can do something like this:

$ sq-seq 1 1237940039285380274899124224 | parallel -j32 'echo {} | step-progression-1 --squares=- 1 1237940039285380274899124224' > data

(The big number is 2^90)

With the right configuration, you can even parallelize over remote machines.

$ sq-seq 1 1237940039285380274899124224 | parallel -S 32/: -S 4/ -S 4/ 'echo {} | step-progression-1 --squares=- 1 1237940039285380274899124224' > data

This does 32 cores locally and 4 remotely at servers, and (ssh must be configured to be able to login without passwords)

This runs the step-progression-1 program once for every line that sq-seq 1 1237940039285380274899124224 generates.

How can I make it go faster? A performance increase can be gained by:

$ seq-morgenstern-mn 3500  > /tmp/indata
$ split -d -n l/6 --suffix-length=3 /tmp/indata indata-
$ cat indata-000 | convert-binary-gmp-records-to-text -i -n2 > indata-000.bin
$ cat indata-001 | convert-binary-gmp-records-to-text -i -n2 > indata-001.bin
$ cat indata-002 | convert-binary-gmp-records-to-text -i -n2 > indata-002.bin
$ cat indata-003 | convert-binary-gmp-records-to-text -i -n2 > indata-003.bin
$ cat indata-004 | convert-binary-gmp-records-to-text -i -n2 > indata-004.bin
$ cat indata-005 | convert-binary-gmp-records-to-text -i -n2 > indata-005.bin
$ mv indata-*.bin /tmp/
$ scp /tmp/indata-*.bin
$ seq -f "%03g" 0 5 | parallel -S 3/ -S 3/:  'cat /tmp/indata-{}.bin | morgenstern-search-type-2 /tmp/indata-{}.bin -f 6 -i' > data

We generate /tmp/indata, and then the split command generates: indata-000, indata-001 and so on. Next we convert them from text to binary gmp records. After manually copying the files into place, we run seq so that we get one number for each of our jobs, and it's the same number of the indata-???.bin file it will operate on.

This method avoids the overhead of regularly setting up and tearing down ssh connections, as well as starting up and tearing down program executions. It has a drawback of the output data not being flushed to data immediately.

It also has the unfortunate drawback of checking the same combination of M,N and P,R twice.

For example, it will check M=2584, N=1597 with P=3448, R=2173, and then later on it will check M=3448, M=2173 with P=2584, R=1597. This means half of our calculations are redundant which is not good at all.

Can we do better than this? Yes! We can use the --mem option.

$ seq-morgenstern-mn 3500  > /tmp/indata
$ seq 0 `cat /tmp/indata | wc -l` > /tmp/records
$ split -d -n r/6 --suffix-length=3 /tmp/records /tmp/records-
$ scp /tmp/indata /tmp/records-*
$ seq -f "%03g" 0 5 | parallel -S 3/ -S 3/:  'cat /tmp/records-{} | morgenstern-search-type-2 /tmp/indata -f 6 --mem' > data

In the first step, we make our data file. And in the second step we generate a sequence of row numbers -- it's a bit fancy with the "`" backticks, but it is just finding the number of lines in the file. In the third step we split the records into 6 files just like before, but now it says r/6 instead of l/6. This causes the row numbers to be round-robined, so that high and low record numbers (indices) are evenly distributed in our files. In the fourth step we copy our files into place. Now, in the fifth step there's a lot going on. We're running seq like in the previous example, where it's one for each job, which is one for each of our record datafiles (e.g. /tmp/records-000, /tmp/records-001 and so on). Parallel is running 3 jobs on this machine (-S 3/:), and 3 on a remote machine. You can see our data file /tmp/indata s an argument to morgenstern-search-type-2. We have passed in the whole data file here because we want to process portions of it, but we're not quite sure yet which portions.

The --mem option reads the whole data file into memory and this is what signals the program to read in row numbers (instead of M,N values) from the standard input, which controls the processing of rows of the data file /tmp/indata.

Recall that in the previous example we were checking combinations twice (e.g. where M,N,P,R was the same as P,R,M,N.), and we're trying to do better than that. With the --mem option, the program gets its first M,N pair as the row number as an index into /tmp/indata, and then we get our P,R pair as every other M,N pair up to the row number. For example, if the record number is 5, we will check the M,N pair in 5 with M,N pairs in 0 through 4. This has the side effect of low record numbers taking much less time than higher record numbers, but that's okay because we round-robined our record numbers in the split command.

The disadvantage of this method is that our whole data file has to be in memory.

Do you want to extend your search? Here's how! In this case we'll be extending our search from 3500 to 5000 and not using the --mem option.

$ seq-morgenstern-mn 3500 > /tmp/3500
$ seq-morgenstern-mn 5000 > /tmp/5000
$ diff --new-line-format="" --unchanged-line-format="" /tmp/5000 /tmp/3500 > /tmp/5000.diff

Now /tmp/5000.diff contains the new values we need to search with. Split and convert, and you're searching for new magic squares.

Extending a search with the --mem option is a little different. We generate the data file like so:

$ seq-morgenstern-mn 5000 > /tmp/5000

But the difference here is that we want our row numbers to start after a certain number of rows (after the rows that belong to 3500):

$ seq-morgenstern-mn 3500 | wc -l
$ cat /tmp/5000 | wc -l
$ seq 2481760 5065076 > /tmp/records

And then we split them like before.

Here's another example of parallelization.

$ seq 1 5000000000000000000 25578614457001433088040  |  parallel --tmpdir=. -S 3/:  'find-3sq-progressions-dist --start={} --range=5000000000000000000' 377850167358576840913440 > data

In this case we know there are several three square progressions that have a delta of 377850167358576840913440. We got this number by repeatedly multiplying 3360 by 13. We use find-3sq-progressions-dist to find progressions with a specific distance, but it would take a long time to find the answer using only a single core, so we employ parallellization.

seq gives us the starting point, which gets passed to the --start option. On the first core, GNU parallel will run this:

find-3sq-progressions-dist --start=1 --range=5000000000000000000' 377850167358576840913440

And on the second core, it runs this:

find-3sq-progressions-dist --start=5000000000000000001 --range=5000000000000000000' 377850167358576840913440

And so on. The results get dumped to data, and we just manually stop when it gets filled with 3 records.

We picked the ending point of seq (which is the 25578614457001433088040 number) by estimating based on prior runs. It is the first value of the first progression from the prior run, multiplied by ten. Maybe it's too big, or maybe it's too small to find the progressions we're looking for. If it's too small, we'll know where we haven't searched and if it's too big, we can stop it early.

One more example to illustrate where a beginner might go wrong:

$ step-progression-1 1 1000 | head -n1 | display-square 
| 1^2 | 5^2 | 7^2 | 75
| 11^2| 145 | 13^2| 435
| 241 | 265 | 17^2| 795
  363   435   507   435

step-progression-1 does not output magic squares! They have to passed through the auto-transform-square program. Also, display-square will show the row and column sums when the square isn't a magic square.

$ step-progression-1 1 1000 | head -n1 | auto-transform-square | display-square 
| 5^2| 241|13^2|
|17^2| 145| 1^2|        magic number = 435
|11^2| 7^2| 265|

How long do these programs take to run?

(times are in seconds)

                                 1 mil    10 mil
step-progression-1               0.344     3.865
step-progression-2              61.806  2004.184
step-progression-3               9.934   317.211
step-progression-4              48.073  1605.358
reverse-step-progression-1       0.358     3.844
fulcrum-progression-1            0.385     3.973
fulcrum-progression-2           62.021  1973.596
fulcrum-progression-3            0.558     6.364
fulcrum-progression-4            0.324     3.281
fulcrum-progression-5           47.827  1612.199
reverse-fulcrum-progression-1    0.384     4.179

Hey, you forgot to mention anything about these sdt programs.

You can see a description of the square diff triangle in the introduction.

The square diff triangle is a method to find arithemetic prorgressions of three squares by subtracting every perfect square from every other perfect square. It is an exhaustive search.

The sdt programs are for doing this exhaustive search in a parallel way.

It's nice because it turns the problem of finding magic squares into a disk-bound sorting problem. It's bad because it doesn't really find a lot of magic squares.

Here is how the "sdt" programs fit together:

$ seq-sdt 1000 > step-1
$ cat step-1 | sdt -o > step-2
$ rm step-1
$ sort-sdt-results -i -o -z 8 -r step-2 > step-3
$ cat step-3 | filter-sdt -i -o -z 8 --stdout > step-4
$ rm step-3
$ cat step-4 | find-3sq-progressions-sdt -i > step-5
$ rm step-4
$ cat step-5 | scissor-square step-5 | check-square > results2
$ cat step-5 | 3sq-pair-search step-5 | count-squares -f6 -f7 -f8 -f9 > results
$ rm step-5
$ cat results results2 | unique-squares | wc -l

This finds 104 unique 3x3 magic squares with 6 perfect squares in a few seconds.

And here's what's going on at each step:

The seq-sdt program generates the spine of the square diff triangle; and the 1000 means we want 1000 rows.

If you take a look at the final row of that output, it shows that it is working with the 10000th square, which has a diff of 999999 when subtracted against 1, and the minimum diff is 1999 when subtracted against the 999th square. Because of the way the square diff triangle works, there are no differences less than 1999 beyond this point in the triangle. If you use a value of 2^16, the largest value in the output will be 2^32.

The sdt program generates the columns for each row. This is what does the subtracting of every square against the square that the row represents.

This can generate a lot of data very quickly. When seq-sdt is run with an argument of 50000, it makes a 24 GB file.

Now comes the sorting. The sort-sdt-results program takes the file and puts it in order, by splitting it up into small chunks, sorting each chunk and then merging the files in a parllel fashion. This is the step that takes a lot of time.

The next step is filter-sdt, and it doesn't do a whole lot but it gets rid of records that can't possibly form a progression of three squrares. This step can also split up file, so if you want to parallelize the next step, it's easy to do that.

The next step is finding the progressions of three squares in the remaining data. Because the data is sorted by the difference of two squares, this program can quickly see if three of them fit together into a progression.

Lastly we run the progressions through the magic square generation programs to find magic squares from the progressions.

This method finds fewer and fewer unique magic squares as it goes forward. If you want to use this method to find 5000 magic squares, it will take more disk space than most people have. I stopped at 1600 unique magic squares, with N=50000, and the curve was starting to level off (approximately 200 new unique magic squares per 10000 rows of sdt). With large values of N, the triangle gets ridiculously big, with a whole lot of repeated magic squares.

How do I build these programs?

You'll need a GNU/Linux system, with the normal sort of tools installed: GNU Bash, GNU GCC, GNU make, GNU automake, GNU autoconf, GNU libtool, and GNU Coreutils.

Also install libgmp (GNU Multiprecision Arithmetic Library.)

All of these packages are available for installation in the package management system of your GNU/Linux distribution.

If you are using Fedora you would type something like this:

    $ dnf install automake autoconf make gcc libtool gmp-devel
(Bash and Coreutils are already installed by default)

    $ git clone http://sv.gnu.org/p/fituvalu/fituvalu.git
    $ cd fituvalu
    $ ./configure
    $ make
    $ su -c 'make install'

Final Words

I wish you good luck in your search for the 3x3 magic square of squares! If you use this software to find anything interesting, please let me know! You can use the Magic Square Checker if you are unsure.

Or if you intend on running this software on a large number of machines in a serious computational effort, also please let me know. I'd love to be in the loop on your exciting endeavor.

If you have any questions about this software, you can contact me by opening a bug on the project page for Fituvalu.

This page is licensed under the terms of the Creative Commons Attribution ShareAlike 4.0 Intl License.