The University of Southern Mississippi

The Aquila Digital Community

Mathematics Undergraduate Theses School of Mathematics and Natural Sciences

12-2014

Magic Surfaces
Brandi Crystal Moore
University of Southern Mississippi

Follow this and additional works at: https://aquila.usm.edu/math_undergraduate_theses

Recommended Citation
Moore, Brandi Crystal, "Magic Surfaces" (2014). Mathematics Undergraduate Theses. 1.
https://aquila.usm.edu/math_undergraduate_theses/1

This Other is brought to you for free and open access by the School of Mathematics and Natural Sciences at The
Aquila Digital Community. It has been accepted for inclusion in Mathematics Undergraduate Theses by an
authorized administrator of The Aquila Digital Community. For more information, please contact
aquilastaff@usm.edu.
 The University of Southern Mississippi

MAGIC SURFACES

by

Brandi Crystal Moore

Abstract of a Thesis
Submitted to the Department of Mathematics
of The University of Southern Mississippi
in Partial Fulfillment of the Requirements
for the Degree of Bachelor of Science

December 2014
 ABSTRACT

MAGIC SURFACES

by Brandi Crystal Moore

December 2014

Magic figures are discrete, two-dimensional (2-D) objects. We translated the definitions
of magic squares and magic circles in an attempt to find three-dimensional (3-D) surfaces
that adhere to the same definitions. This is what we defined as a magic surface.
First, we translated these definitions by thinking of the discrete integer entries in the
magic squares and magic circles as volumes under the magic surfaces; these were evaluated
by volume integrals. We then translated these definitions in a similar way to apply line
integral constraints.
We found polynomial functions that satisfied thse re-definitions of conditions for a
magic surface. In the case over magic squares using volume integrals, we were able to form
conjectures about the polynomial solutions and the systems of equations that formed them.

iii
 C OPYRIGHT BY

B RANDI C RYSTAL M OORE

2014
 The University of Southern Mississippi

MAGIC SURFACES

by

Brandi Crystal Moore

A Thesis
Submitted to the Department of Mathematics
of The University of Southern Mississippi
in Partial Fulfillment of the Requirements
for the Degree of Bachelor of Science

Approved:

Director

Chair, Department

December 2014
 ACKNOWLEDGEMENTS

Thank you Dr. John Harris for giving me the concept of this research in the first place.
And for poking your head in with funny stories and interesting opinions on my color choices.
Thank you Dr. Jeremy Lyle for making my life more difficult, for about a month. Your
idea on making the function orthogonal to give us the uniqueness we were hoping to achieve
was well-intended. Unfortunately, Sage and Mathematica were unable to solve the non-linear
equations. Maybe we can explore it more in the future.
Thank you Mother for not letting your eyes glaze over during all my attempts to explain
this project to you. And thank you for not yelling at me when the dishes did not get done
because I was working on it.
Thank you Dawn, my sister, for letting me vent my frustrations about this project and
paper while boosting my confidence.
Finally, thank you Dr. Karen Kohl, the other half of the "we", for holding my hand
through this entire project. Thank you for teaching me parts of calculus my education lacked.
Thank you for explaining the same concept three times before I finally understood it, and
not getting visibly frustrated with me. Thank you for being available to communicate and
meet numerous times throughout this process. And thank you for talking about our pets and
puzzles and everything in-between.

v
 TABLE OF CONTENTS

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii

ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Magic Figures 1
1.2 Magic Surfaces 4

2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1 Overlay 5
2.2 Storage 5
2.3 Algorithm 6

3 Magic Squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.1 Volume 7
3.2 Line Integrals 11

4 Magic Circles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.1 Volume 15
4.2 Line Integrals 18

5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5.1 Results 20
5.2 Future Work 20

APPENDIX

A Sage Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
A.1 Common Functions 21
A.2 Functions for volume over a square 23
A.3 Functions for line integrals over a square 25
A.4 Functions for volume over a circle 27
A.5 Functions for line integrals over a circle 30

B Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
B.1 Solutions for volume over a square 33
B.2 Solutions for line integrals over a square 45

vi
B.3 Solution for volume over a circle 47
B.4 Solution for line integrals over a circle 50

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

vii
 1

Chapter 1

INTRODUCTION

Calculus has been an integral part of the emergence of the modern world. Before tangent
lines and areas under a curve interested mathematicians, recreational mathematics had been
intriguing both math professionals and those with no mathematical training at all.
We will study discrete, 2-D magic figures: both magic squares and magic circles.
Throughout the course of this research, we will bridge integral calculus and recreational
mathematics.

1.1 Magic Figures

Magic figures will be the portion of recreational mathematics that we will study. These are
two-dimensional (2-D), discrete structures that have certain properties that must be met that
involve the placement of integers into various shapes. All magic figures have conditions on
the placement of integers subject to the sums of those integers.
We will extend these properties and definitions to three-dimensional (3-D) continuous
surfaces in all respects.

1.1.1 Magic Squares

Magic squares are just one of these magic figures. Anderson [1] explains magic squares’
long and diverse history beginning around 2200 B.C. Magic squares have been used for
astronomy, divination, alchemy, interpreting philosophy, and medicine in many cultures
including Chinese, Islamic, Indian, West African, and European. The Arabic term for magic
squares translates to "harmonious distribution of the numbers."
In general, they are an n × n square that meets certain conditions:

1. Entries are consecutive, unique integers that range from 1 to n2 .

2. The sum of the entries in each row and column and diagonal are all equal.

This is an example of a 3 × 3 magic square:

 2

This figure shows the summing conditions: the entries in each row and each column and
each diagonal sum to the same constant, 15.
Definition 1. The magic constant is the sum of each entry in each row, column, and diagonal,
given by [2]
2
1 n 1 n2 (n2 + 1) 1 2
M2 (n) = ∑ k = = n(n + 1) (1.1)
n k=1 n 2 2
We can use this constant to give a more concrete definition of those summing conditions.
If Mx,y is the xyth entry in the n × n magic square.
Let Hi be the ith horizontal row:

Hi = Σnj=1 Mi j (1.2)

And let V j be the jth vertical column:

V j = Σni=1 Mi, j (1.3)

Definition 2. For an n × n square to be semi- magic both conditions listed must be satisfied:

1. Hi = M2 (n) ∀1 ≤ i ≤ n

2. V j = M2 (n) ∀1 ≤ j ≤ n

Definition 3. For the same square to be magic, it must be semi-magic and

1. Σi Mi,i = M2 (n)

2. Σi Mi,n−i+1 = M2 (n)

for the main diagonal and opposite diagonal, respectively.

 3

Magic squares have been studied for over 4 millennia, so we know quite a bit about them.
For instance, the smallest, non-trivial magic square is the 3 × 3 as there is no 2 × 2. There
are different types of magic squares aside from the normal ones we are using: panmagic,
bimagic, multimagic, alphamagic. There are some algorithms that can generate the magic
square by hand, such as the LUX method, lozenge method, and the Siamese method.
Molony [4] has written computer code that can generate these different types of magic
squares. His code is broken into cases based on the sizes of magic squares: odd, singly even,
and doubly even.
Definition 4. A doubly even number is a number divisible by 2 and 4.
Definition 5. A singly even number is a number divisible by 2 and not 4.
The number of distinct magic squares (those not generated through rotation and re-
flection) has been calculated for a few sized magic squares: 1, 880, 275305224 for 3 × 3,
4 × 4, and 5 × 5, respectively. The number of all 6 × 6 magic squares is unknown, though
it is estimated to be around 1.7745 ×1019 . Weisstein [7] also notes that "it is an unsolved
problem to determine the number of magic squares of an arbitrary order." These open
questions about magic squares make them all the more interesting.

1.1.2 Magic Circles

Another type of magic figure is the magic circle[6]. Although there are two types of magic
circles, only one is interesting to our work. This type of magic circle has entries placed
around a center in rings along 8 radial lines. Therefore, the number of rings must divide the
number of entries less one by 8.
Definition 6. The magic circle conditions are:

1. Entries are consecutive, unique integers ranging from 1 to n

2. The sum of the entries around a ring sums to the magic circle constant

3. The sum of the entries across diameter lines sums to the magic circle constant,
excluding the entry placed in the center

We denote the number of rings in a magic circle by k = n−1 8 .

For example, this figure is a 4-ring magic circle with entries 1, 2, 3, ..., 33.
 4

33−1
• k= 8 =4

• The innermost ring sums to 3+24+10+23+16+20+27+15=138

• The vertical diameter line sums to 30+7+22+10+27+28+8+6=138

There are many magic circles of the same size; different entries in the center can lead to
many more permuatations of the entries. Just as with magic squares, no formula is known
for the number of magic circles of a given size.

1.2 Magic Surfaces

Using these definitions of magic squares and circles, we will extend these two-dimensional,
discrete objects into three-dimensions. We will adhere to the constraints about integer
placement and summing of those entries by evaluating both volume and line integrals. We
will search for continuous and smooth functions of 2 variables that satisfy a redefinition of
the magic properties.
To do so, we will view f (x, y) as a surface over a magic square placed in the first quadrant
with the bottom left corner at the origin. In one case, we will consider each entry to be a
quantity of the volume integral over a 1 × 1 subsquare of the xy-plane. Alternatively, we will
view the entries in such subsquares as the evaluation of the line integral in both directions.
We will also study functions over magic circles, with the center placed at the origin of the
plane and using polar coordinates, by evaluating the volume and line integrals over sections
set equal to the corresponding entries.
 5

Chapter 2

Methodology

Throughout the course of this research the computer algebra system Sage was used to do
the searching for solutions and graphing. The use of Mathematica’s interface through Sage
was considered briefly, but, in the end, Sage was able to solve the systems of equations
using its own processes. In addition, the Sage server proved to be more user-friendly. After
getting the solutions to the equations, we were able to substitute those solutions back into
the equation and graph.

2.1 Overlay

To adhere to the constraint of unique integers being represented in a magic figure, we will
view the surface over a magic figure in the xy-plane. For the square, we will place the
function over the first quadrant of the Cartesian plane. In the case of the magic circles, we
imagined the center of the circle at the origin.

2.2 Storage

We stored the magic squares and circles as matrices in Sage.

 
8 1 6
3 5 7 (2.1)
4 9 2

is a 3 × 3 magic square, and

 
9 0 0 0 0 0 0 0
 3 24 10 23 16 20 27 15
 
11 25 22 13 1 33 28 5  (2.2)
 
26 18 7 19 31 12 8 17
29 2 30 14 21 4 6 32

is a 4-ring magic circle.

Here, note that Sage starts indexing matrices at 0. Therefore, M[2, 1] for the magic
square’s matrix is 7.
 6

2.3 Algorithm

We used a magic square generator written by Molony [4] that created a magic square of a
given degree represented as a matrix using different algorithms for the different cases of
sized squares. No algorithm was used for magic circle creation; they were input by hand.
To search for a magic surface function, we developed the following algorithm:

1. generate a polynomial of a given degree, p, with 2 specified variables and unknown

coefficients using create_poly, a function we wrote

2. Set up equations by setting each integral evaluated over a certain region equal to the
corresponding entry from that figure

3. Attempt to solve for the unknown coefficients in create_poly

4. If no solution was found, increase the degree of the polynomial, p, and repeat the
algorithm

After the algorithm found a solution, we produced a graph of the solution.

Here is some sample output from create_poly:

k3*x^3 + k6*x^2*y + k9*x*y^2 + k12*y^3 + k2*x^2 + k5*x*y

+ k8*y^2 + k1*x + k4*y + k0

[k0, k4, k8, k12, k1, k5, k9, k2, k6, k3]

This includes the polynomial produced and the list of unknown coefficients.
 7

Chapter 3

Magic Squares

3.1 Volume

To make the surface generated by the polynomial in x and y of degree p magic over a square,
it must satisfy certain conditions:

1. The volume integral over the ith horizontal row equals the magic constant: Hi = M2 (n)

2. The volume integral over the jth vertical column equals the magic constant: V j =
M2 (n)

3. The volumes over subsquares are consecutive, unique integers

To accomplish this, we will start with a known n × n magic square. We imagine the
integer entries as being the volume under the magic surface in that subsquare’s specific
region. To carry over the uniqueness of the integer entries in the magic square to the 3-D
figure we hope to find, we will evaluate the volume integral over a subsquare and set it equal
to the integer entry overlaid.

The constraint equations in step 2 of the algorithm will be of the form

Z n−i Z j+1
f (x, y)dxdy = M[i, j] (3.1)
n−i−1 j
 8

where n is the dimension of the matrix that represents the square, and 0 ≤ i, j < n. Due
to the additive interval (region) property of integration [3], this ensures the volume over a
horizontal or vertical strip is equal to the magic constant.
There will be n2 such equations will be of this form since we are evaluating the vol-
ume integrals over each subsquare and there are n2 subsquares. For example, we want
R3R1
2 f (x, y)dxdy = 8 = M[0, 0]. This is the top left subsquare of the 3 × 3 magic square.
0
We evaluate this region first because we need it set equal to the M[0, 0] entry in the matrix:
the first entry of the matrix lies in the top lefthand corner of the grid. We go through the rest
of the subsquares and set them equal to their corresponding matrix entry to solve for the
unknown coefficients.

3.1.1 3 × 3 Example

For example, we represent the magic square shown as the matrix below. Note that
Sage/Python indexes rows and columns starting with 0; this determined the limits of
integration in the constraint equation.

 
8 1 6
3 5 7 (3.2)
4 9 2
We will search for a polynomial, f (x, y), that when placed over the figure above, satisfies
the volume integral constraint.
For the evaluation of this 3 × 3 magic square’s volume,

• We start by assuming f (x, y) is a degree 3 polynomial in variables x and y. The

automatically generated polynomial is f (x, y) = k3 x3 + k6 x2 y + k9 xy2 + k12 y3 + k2 x2 +
k5 xy+k8 y2 +k1 x+k4 y+k0 , with unknown coefficients [k0 , k4 , k8 , k12 , k1 , k5 , k9 , k2 , k6 , k3 ].

• The system of 9 equations is set up by evaluating the volume integrals of the polyno-
mial by looping through values of i and j.
1
• The solution of the system gave f (x, y) = r1 x3 + 11 (2r2 +45)y3 − 92 (r1 +1)x2 +3x2 y−
9 2 2 1 3 3 247
22 (2r2 + 23)y − 6xy + 2 (11r1 + 9)x + r2 y + 9xy − 2 r1 − 11 r2 + 44 where r1 and r2
are free variables.

• After choosing to set r1 = 0 and r2 = 0, we have the solution given by f (x, y) =

45 3 9 2 2 207 2 2 9 247
11 y − 2 x + 3x y − 22 y − 6xy + 2 x + 9xy + 44 .
 9

This is the graph of the solution:

The red dots on the surface are at the height of the entry of the subsquare in the center of
the subsquare being represented. This shows how the height of the surface varies throughout
the entire region being evaluated, with the average height of the subsquare being equal to
the entry.

3.1.2 Results

Using the search_mtx_sq_vol function (Appendix A) for magic squares 3 × 3 through

12 × 12, (Appendix B) we were able to notice some patterns. We assume the pattern holds
for all magic squares of the same size. We used the magic square test function by Molony
[4] to test all 3 × 3 magic squares; it would take an exceptionally long amount of time to try
to generate the 4 × 4 magic squares, so we generalized from the 3 × 3 results. We confirmed
that all magic squares of size 3 × 3 have the same results shown in the table below.
 10

Square Size Equations Degree Unknowns Free Variables Unknowns – Free Vars.
(p+1)(p+2)
n n2 p 2
3 9 3 10 2 8
4 16 5 21 6 15
5 25 7 36 12 24
6 36 10 66 30 36
7 49 11 78 30 48
8 64 13 105 42 63
9 81 15 136 56 80
10 100 18 190 90 100
11 121 19 210 90 120
12 144 21 253 110 243

From the table of results, we formed several conjectures. Note that Conjecture 1 implies
Conjecture 2.
Conjecture 1. The degrees of polynomials needed to solve a system of equations when
evaluating the volume integrals over a n × n magic square increase in a pattern of 2, 2, 3, 1,
2, 2, 3, 1, ... starting at the degree of the 3 × 3 magic square, 3.
Conjecture 2. The degrees of polynomials are odd unless the size of the n × n magic square
is singly even.
Conjecture 3. There only need be one less non-free unknown than equations in order to
solve the system, unless the size of the n × n square is singly even.
These singly even magic squares are interesting, not only because they change these
patterns, but because they also are built differently than the odd and doubly even magic
squares. The most common, if not only, way to generate a singly even magic square is first
by sectioning it half vertically and horizontally, then filling in these sections by the Siam
method [5]. Finally, since this is not a magic square yet, elements in almost every column
must be exchanged with another until the magic square is achieved.
We noticed that the graphs of solutions are exaggerated more than expected. Some fall
below the z-axis much further away from the xy-plane than we had thought since we had
pictured the surfaces above the magic square. The increase in the entries made an increase
in the "waves" seen on the surface up until the 7 × 7; once the squares reached this size, they
became what looked to be relatively flat with a sharp angle near one of the borders of the
square.
 11

3.1.3 Non-Magic Evaluation

After noticing these patterns, we tested some non-magic squares for the same conditions.
These squares merely contained consecutive, unique integers; the summing conditions of
the entries were not met.
Square Size Equations Degree Unknowns Free Variables Unk.-Free Vars.
n n2 2(n − 1) n(2n − 1) n(n − 1)
3 9 4 15 6 9
4 16 6 28 12 16
5 25 8 45 20 25
6 36 10 66 30 36
As shown in the table, these non-magic squares with the same sizes as the magic squares
we evaluated do not have the same conjectured properties. The number of non-free unknowns
is equal to the number of equations in each of these few cases, whereas the magic squares
evaluated had one less non-free unknown than equations unless the size of the magic square
was singly even. Another difference between the magic and non-magic squares evaluated
is that all of the magic squares were given an odd degree of polynomial (unless, again, if
the square’s size was singly even) for the surface, but the non-magic squares’ surfaces were
produced with even degrees, regardless of size of square.
From these observations on a small sample size, we conclude that the surfaces produced
and the patterns noted truly make for a magic surface in the volume integral constraint over
a magic square. A more thorough study of surfaces over non-magic squares needs to be
carried out.

3.2 Line Integrals

The other way we evaluated whether a surface has magic properties is by testing the
line integrals, with respect to arc length, down the center of the row/column. The magic
 12

conditions for the line integral case are:

1. The kth horizontal line integral across a row equals the magic constant: Hk = M2 (n)

2. The l th vertical line integral down a column equals the magic constant: Vl = M2 (n)

3. The diagonal and anti-diagonal equal the magic constant: D1 = D2 = M2 (n)

4. The line integrals over the subsquares are consecutive, unique integers

Just as with the volume integral case, we overlaid the surface above a magic square and set
the evaluated integrals over a region equal to the entry in that subsquare.

Recall from [3] that the line integral of f (x, y) with respect to arc length along a path C
is deonted by Z
f (x, y) ds, (3.3)
C
where x(t) and y(t) parameterize C, and
s 2  2
dx dy
ds = + dt (3.4)
dt dt
The horizontal line integrals over subsquares are given by
s 
dx 2
Z b
ds = + 0 dt ds = dx f (x, y)dx (3.5)
dt a

This produced n2 equations of the form

Z j+1
f (x, n − i − 1/2)dx = M[i, j], 0 ≤ i, j ≤ n (3.6)
j

The limits of integration are the changing x-values as the y-value stays the same; the y-value
is written as n − i − 12 because we looped through the i’s until all y-values down the middle
1
(where the 2 comes from). These limits are as such since the first entry in the matrix is in
the top left subsquare of the magic square.
The vertical line integrals over subsquares are given by
s  2 Z b
dy
ds = 0 + dt ds = dy f (x, y)dy (3.7)
dt a

This again produced n2 equations of this form

Z n−i
f ( j + 1/2, y)dy = M[i, j], 0 ≤ i, j ≤ n (3.8)
n−i−1
 13

Here, the limits of integration are the changing y-values. The x-values are looped through
going down the center of the subsquares, j + 12 . Again, because the matrix representing the
magic square begins in the top.
Unlike the volume integrals, we adhere more strictly to the definition of a magic square
by also evaluating the diagonals of each region and setting it equal to the entry in that
subsquare.
We parameterize C for the main diagonal represented by the line y = x by

x(t) = t, y(t) = x = t (3.9)

We parameterize C for the anti-diagonal represented by the line y = −x + n by

x(t) = t, y(t) = −t + n (3.10)

So here, ds for the diagonals is

s   
dt 2 dt 2 √
q
ds = + dt = (1)2 + (±1)2 dt ds = 2dt (3.11)
dt dt

For the diagonal, the n equations are of the form

Z i+1 √
f (t, −t + N) 2dt = M[i, i]. (3.12)
i

Similarly, the n equations for the anti-diagonal are of the form

Z n−i √
f (t,t) 2dt = M[i, n − i − 1]. (3.13)
n−i−1

3.2.1 3 × 3 Example

We evaluated the same square for the line integral constraint as we did the volume integral
constraint,  
8 1 6
3 5 7 (3.14)
4 9 2
A 6th degree polynomial was needed to solve these 24 equations; this output can be
found in Appendix B.
This is the graph of the solution.
 14

Here, the red lines drawn on the surface represent the vertical line integrals, and the green
lines represent the horizontal integrals.

3.2.2 Results

Square Size Equations Degree Unknowns Free Variables Unknowns-Free Variables

n 2n(n + 1) 2n (2n + 1)(n + 1) n+1
3 24 6 28 4 24
4 40 8 45 5 40
5 60 10 66 6 60
6 84 12 91 7 84

Unlike the volume integrals, it seems that the patterns from the volume case do not carry
over. Here, the increase in degree is always by 2. Also, the number of unknowns minus the
number of free variables is equal to the number of equations.
 15

Chapter 4

Magic Circles

4.1 Volume

We will redraw magic circles so that the entries are now representing volumes of sections of
the circle:

To make the surface in polar coordinates magic over a circle, it must satisfy certain
conditions:

1. The volume integral of a wedge over a diameter line equals the magic constant

2. The volume integral of a washer over a ring equals the magic constant

3. The volumes over sections are consecutive, unique integers

As with the squares, we imagine the circle’s entries as being the volume under the magic
surface in that section’s specific region. We stored these entries as a matrix with the center
being the only entry in the top row; the rest of the matrix was filled in starting at the first
ring on angle 0 and moving counter clockwise.
 
9 0 0 0 0 0 0 0
 3 24 10 23 16 20 27 15
 
11 25 22 13 1 33 28 5  (4.1)
 
26 18 7 19 31 12 8 17
29 2 30 14 21 4 6 32
 16

The volume integral constraints over a section represented by an entry in the matrix have
this form: iπ π
4 +8
Z Z j+1/2
iπ π
f (r, θ )rdrdθ = M[ j, i] (4.2)
4 −8 j−1/2

We decided to keep the center for the volume case though we know there is a point of
singularity there. We searched through many different types of polynomials, f (r, θ )

• with terms km ri sink hθ

• with terms km ri cos j hθ

• with terms produced by (r + θ ) p × (sin θ + cos θ )

• with terms km ri θ j

It was finally decided to use regular polynomials in r and θ , km ri θ j , since we were unable
to produce solutions with the trigonometric functions included.
We were able to get solutions, but when they were graphed there was a discontinuity at
θ = 0 = 2π. To remedy this, we set the limit of the polynomial as θ → 0 equal to the limit
of the polynomial as θ → 2π

lim f (r, θ ) = lim f (r, θ ) (4.3)

θ →0+ θ →2π −

Then we looped through the derivatives of the polynomial and took the limits of them to
make the function smooth.

lim f (d) (r, θ ) = lim f (d) (r, θ ) (4.4)

θ →0+ θ →2π −

where 1 ≤ d ≤ p, and p is the degree of the polynomial.

We are matching the coefficients in these equations and appending these new equations
into the system of equations of the volume integrals.

4.1.1 One Result

This is the result of evaluating volume integrals over sections set equal to the entries of the
circle in
 17

The solution required the polynomial f (r, θ ) to be of degree 11. The number of equations,
99, far exceeded the number of unknowns, 78, needed to solve the system of equations.
There are so many equations because not only are the volume integrals included in the
system, but also all of the equations created when matching the coefficients in the limits
of the derivatives. Since the degree of the polynomial is 11, we can take the derivative 11
times; this gives us 12 new equations (the 11 derivatives and the original function). There
were 35 free variables, so the number of unknowns actually needed to solve the equation
is 78-35=33; this happens to be the number of entries, n. The full solution can be found in
Appendix B.
This is the graph of the solution:

Even after including the new constraints to make the function smooth and continuous
at θ = 0 = 2π, there is still an odd dip in the graph. There is a discontiniuity at the origin,
though we evaluated the volume integral in that region and set it equal to the central entry.
 18

4.2 Line Integrals

This is our reimagining of a magic circle while taking the line integrals:

For example, if we were to look at the paths over 18, we would take the line integral
from θ = π8 to θ = 3π
8 (between green dotted lines) along r = 3 (the blue arc) and set it
equal to 18 then take the line integral from r = 2.5 to r = 3.5 (between red dotted lines)
π
along the radial line θ = 4 (the purple line).
Recall that the line integral is given by
Z
f (r, θ )ds (4.5)
C

where C is either a path along a radius or a path along an arc. In polar coordinates, we have
p
ds = dr2 + r2 dθ 2 (4.6)

and we will simplify this in the two cases of paths.

When r is constant, θ varies:
p
ds = 0 + r2 dθ 2 ds = rdθ (4.7)

so the line integral constraint along an arc segment of a ring is

Z (2 j+1) π
8
f (i, θ )rdθ = M[ j, i]. (4.8)
(2 j−1) π8

When θ is constant, r varies:

p
ds = dr2 + 0 ds = dr (4.9)
 19

so the line integral constraint along a radial line is

Z (2 j+1) 1
2 iπ
f (r, )dr = M[ j, i]. (4.10)
(2 j−1) 12 4

4.2.1 One Result

We were able to find a solution to the line integral constraint over the same magic circle we
evaluated the volume integral over. The degree of the polynomial needed was 12, which
produced 91 unknowns. The system contained 155 equations, but there were only 15 free
variables in the solution. This is the only case where there were more equations than non-free
unknowns. The solution can be found in Appendix B.
This is the graph of the solution:

There is the discontinuity at the origin, just as with the volume integral case. In this line
integral case, we did not evaluate over the central entry.
 20

Chapter 5

Conclusion

5.1 Results

We were able to translate the definition of these discrete magic figures to generalize magic
surfaces in these four cases. To do so, we wrote an algorithm to search for such surfaces
over a magic figure of arbitrary size.
We were also able to notice patterns in the magic square cases; these included the degrees
of the polynomials and the number of unknowns needed to solve the system dependent on
the size of the square. The non-magic squares we tested did not have the same properties as
the magic squares of the same sizes.

5.2 Future Work

The magic surfaces produced in this text have only been satisfying to some degree. There is
still more work to be done; a lifetime could be spent exploring what we have not been able
to and verifying what we have.
One priority would be to place additional constraints on the magic surfaces so that the
surfaces produced would be strictly above the xy-plane or rθ -plane and bounded above.
This would correct the appearance of "flat" areas in our graphs. We would want to be able to
look closer at the curvature of the surface to see what other interesting properties we can
identify.
We would also want to see if magic squares and circles exist that satisfy both the volume
and line integral cases. Another option would be to relocate the discrete squares in an
attempt to find surfaces that satisfy more than one than one magic square of different sizes.
Upon further research, we would want to test what different kinds of functions satisfy
these magic constraints, especially those that are non-polynomial.
A final goal would be to further generalize these magic surfaces by getting away from
these discrete magic figures while still retaining uniqueness in some other way.
 21

Appendix A

Sage Functions

A.1 Common Functions

d e f create_vars ( viter , varstr , num ) :

# c r e a t e t h e l i s t num o f v a r i a b l e s u s i n g t h e i t e r a t o r v i t e r
r e t u r n [ var ( varstr + s t r ( viter . n e x t ( ) ) ) f o r i i n range ( num ) ]

d e f create_poly ( deg , v1 , v2 ) :
# c r e a t e s a p o l y n o m i a l o f d e g r e e deg i n t h e 2 v a r i a b l e s s p e c i f i e d
num =(1+ deg ) ^2
var ( ' r T ' )
K=[]
citer = i t e r ( range ( 1 0 0 0 ) )
L=create_vars ( citer , ' k ' , num )
poly =0
m=0
f o r i i n range ( deg + 1 ) :
f o r k i n range ( deg + 1 ) :
i f ( i+k ) <=deg :
m+=1
t=L [ m ] * v1 * * ( i ) * v2 * * ( k )
poly=t+poly
K . append ( L [ m ] )

r e t u r n poly , K

d e f matcher ( f , deg ) :
# matches c o e f f i c i e n t s of degrees of r to b e t t e r s o l v e f o r k_i ' s
# t a k e s i n t h e f u n c t i o n and t h e h i g h e s t d e g r e e o f t h a t f u n c t i o n
f ( r , T ) =f
R=f ( r , T ) . limit ( T =0 , d i r = ' + ' )
# t a k e s t h e l i m i t o f t h e f u n c t i o n as t h e t a goes t o 0 from t h e r i g h t
L=f ( r , T ) . limit ( T=2 * pi , d i r = '− ' )
# t a k e s t h e l i m i t o f t h e f u n c t i o n a s t h e t a g o e s t o 2 * p i f r o m t h e ←-
left
# t h i s i s t o remove t h e d i s c o n t i n u i t y on t h e c i r c l e where 0=2 * p i
eq = L==R
 22

# l o o p t h r o u g h d e g r e e , match c o e f f i c i e n t s o f p o w e r s o f r t o ←-
p r o d u c e a l i s t o f new e q u a t i o n s
eqlist = [ ]
f o r i i n range ( deg + 1 ) :
#+1 b e c a u s e r a n g e g o e s f r o m 0 t o d−1
eqn=L . coefficient ( r , i ) ==R . coefficient ( r , i )
# s e t s t h e c o e f f i c i e n t s on e i t h e r s i d e o f t h e e q u a t i o n e q u a l t o ←-
each o t h e r
eqlist . append ( eqn ) # a d d s t h e new e q u a t i o n t o t h e l i s t

r e t u r n eqlist
# r e t u r n s t h e l i s t o f new e q u a t i o n s

d e f polar_plot3d ( f ) :
# p l o t s f u n c t i o n i n 3−d p o l a r c o o r d i n a t e s
f ( r , T ) =f
print f(r , T)
r e t u r n parametric_plot3d ( ( r * cos ( T ) , r * sin ( T ) , f ( r , T ) ) , ( T , 0 , 2 * pi ) , ( ←-
r , 0 , 4 ) , mesh=True , dots=True , plot_points = [ 1 0 0 , 1 0 0 ] )
 23

A.2 Functions for volume over a square

d e f search_mtx_sq_vol ( f , L , M , pr=False , solver= ' s a g e ' ) :

# s e a r c h e s f o r a f u n c t i o n o f g i v e n f o r m f o v e r a magic s q u a r e ←-
r e p r e s e n t e d as a m a t r i x M
#L = l i s t o f unknown p a r a m e t e r s
f ( x , y ) =f
N = M . dimensions ( ) [ 0 ]
EqnList = [ ]
f o r i i n range ( N ) :
f o r j i n range ( N ) :
h=integrate ( integrate ( f , x , j , j + 1 ) , y , N−i −1 ,N−i )
# c a l c u l a t e s subsquare volume
i f pr :
print i , " , " ,j , " : " , h , "=" , M [ i , j ]
EqnList . append ( h==M [ i , j ] )
# adds t h e volume c o n s t r a i n t t o t h e l i s t

# creates a dictionary of solutions

i f solver== ' s a g e ' :
sol = solve ( EqnList , L , solution_dict=True )

i f sol = = [ ] :
d = sol
r e t u r n "No s o l u t i o n "
else :
d=sol [ 0 ]
# s o l v e i n Sage
else :
d = solvemma ( EqnList , L )
# s o l v e in Mathematica

fsubs = f . subs_expr ( d )
# s u b s t i t u t e into f the assignments given in dictionary d

freevars = Set ( fsubs . variables ( ) ) . difference ( Set ( [ x , y ] ) )

# substitute that solution into f
f o r freevar i n l i s t ( freevars ) :
var ( freevar )

p r i n t " # e q u a t i o n s = " , l e n ( EqnList )

p r i n t " # unknowns = " , l e n ( L )
p r i n t " # f r e e v a r s = " , l e n ( freevars )
r e t u r n fsubs , freevars

d e f search_sq_vol_deg ( M ) :
# loops through degree of polynomial u n t i l s o l u t i o n i s found
deg =3
 24

sol= ' No s o l u t i o n '

w h i l e ( sol== ' No s o l u t i o n ' and deg <=25) :
poly=create_poly ( deg )

sol=search_mtx_sq_vol ( poly [ 0 ] , poly [ 1 ] , M )

print " "
deg+=1
p r i n t " deg = " , deg−1
r e t u r n sol

d e f conv_matrix_to_points ( M ) :
# converts the matrix M i n t o points in the center of the subsquare
N = M . dimensions ( ) [ 0 ]
L = []
f o r i i n range ( N ) :
f o r j i n range ( N ) :
x = j + 1/2
y = N−1−i + 1 / 2
z = M[i , j]
L . append ( ( x , y , z ) )
return L

d e f solve_and_graph ( N ) :
# searches f o r a s o l u t i o n then graphs i t
M=make_magic_square ( N )
print M
f=search_sq_vol_deg ( M )
print f [0]
# f u n c t i o n t h a t has a l r e a d y been s u b s t i t u t e d w i t h s o l u t i o n s
print f [1]
# l i s t of free variables
fvlist=f [ 1 ]
d = {}
# c r e a t e s an e m p t y d i c t i o n a r y
f o r r i n fvlist :
# i t e r a t e s through the l i s t of f r e e v a r ia b l e s
d [ r ]=0
# s e t s each f r e e v a r i a b l e i n d i c t i o n a r y equal t o z e r o
print d
func=f [ 0 ] . subs_expr ( d )
# s u b s t i t u t e s dictionary assignments into the function

pts=point3d ( conv_matrix_to_points ( M ) , size =10 , color= ' r e d ' )

graph=plot3d ( func , ( x , 0 , N ) , ( y , 0 , N ) )
show ( graph+pts )
 25

A.3 Functions for line integrals over a square

d e f search_mtx_sq_line ( f , L , M , pr=False , solver= ' s a g e ' ) :

# s e a r c h e s f o r a f u n c t i o n o f g i v e n f o r m f o v e r a magic s q u a r e ←-
r e p r e s e n t e d as a m a t r i x M
#line integrals
#L = l i s t o f unknown p a r a m e t e r s
f ( x , y ) =f
var ( ' t ' )
print f
N = M . dimensions ( ) [ 0 ]
EqnList = [ ]

f o r i i n range ( N ) :
f o r j i n range ( N ) :
h=integrate ( f ( x , N−i − 1 / 2 ) , x , j , j + 1 )
v=integrate ( f ( j + 1 / 2 , y ) , y , N−i −1 ,N−i )
i f pr :
print " Horiz : " ,i , " , " ,j , " : " , h , "=" , M [ i , j ]
print " Vert : " ,i , " , " ,j , " : " , v , "=" , M [ i , j ]
EqnList . append ( h==M [ i , j ] )
EqnList . append ( v==M [ i , j ] )

f o r i i n range ( N ) :
# main d i a g o n a l
d1 = integrate ( f ( t,−t+N ) * sqrt ( 2 ) , t , i , i + 1 )
EqnList . append ( d1==M [ i , i ] )

# a n t i −d i a g o n a l
d2 = integrate ( f ( t , t ) * sqrt ( 2 ) , t , N−1−i , N−i )
EqnList . append ( d2==M [ i , N−1−i ] )

p r i n t " # e q u a t i o n s = " , l e n ( EqnList )

p r i n t " # unknowns = " , l e n ( L )

# creates a dictionary of solutions

i f solver== ' s a g e ' :
sol = solve ( EqnList , L , solution_dict=True )
i f sol = = [ ] :
r e t u r n ' no s o l u t i o n '
else :
d=sol [ 0 ]
# s o l v e i n Sage
else :
d = solvemma ( EqnList , L )
# s o l v e in Mathematica
print "d =" , d

fsubs = f . subs_expr ( d )
# s u b s t i t u t e into f the assignments given in dictionary d
 26

freevars = Set ( fsubs . variables ( ) ) . difference ( Set ( [ x , y ] ) )

# substitute that solution into f
f o r freevar i n l i s t ( freevars ) :
var ( freevar )

r e t u r n fsubs , freevars

d e f plot_li ( f , d , N ) :
# draws t h e l i n e i n t e g r a l s on s u r f a c e
fsubs = f . subs_expr ( d )
# s u b s t i t u t e s dictionary values into f
fsubs ( x , y ) = fsubs
p r i n t " f s u b s = " , fsubs
g= plot3d ( fsubs ( x , y ) , ( x , 0 , N ) , ( y , 0 , N ) )

width = . 0 4
# w i d t h o f h i g h l i g h t on s u r f a c e

f o r c i n range ( N ) :
ctr =c + . 5
below = ctr − width / 2
above = ctr + width / 2
vert=plot3d ( fsubs ( ctr , y ) , ( x , below , above ) , ( y , 0 , N ) , rgbcolor←-
=(1 ,0 ,0) )
g=g+vert
horiz=plot3d ( fsubs ( x , ctr ) , ( x , 0 , N ) , ( y , below , above ) , rgbcolor←-
=(0 ,1 ,0) )
g=g+horiz
show ( g )
 27

A.4 Functions for volume over a circle

d e f search_mtx_cir_vol ( f , L , M , N , deg , solver= ' s a g e ' ) :

# s e a r c h e s f o r s o l u t i o n t o volume i n t e g r a l s over a c i r c l e
k = ( N−1) / 8
f ( r , T ) =f
rlines=pi / 4
# r a d i a l l i n e s are at every m u l t i p l e of pi /4
split=pi / 8
# s p l i t s t h e d i f f e r e n c e b e t w e e n e a c h r a d i a l l i n e t o g e t p r o p e r ←-
interval
dist = 1 / 2
#how f a r up and down t h e r a d i a l l i n e t o go on e i t h e r s i d e o f e n t r y ←-
for a complete i n t e r v a l of 1
EqnList = [ ]
f o r i i n range ( 8 ) :
f o r j i n range ( 1 , k + 1 ) :
v=integrate ( integrate ( f ( r , T ) * r , T , i * ( pi / 4 ) −(pi / 8 ) , i * ( pi / 4 ) ←-
+ ( pi / 8 ) ) , r , j −1/ 2 , j + 1 / 2 )
eq = v==M [ j , i ]
i f n o t eq :
# i f e q u a t i o n i s n o t " F a l s e " ( i n c o n s i s t e n t ) , t h e n append ←-
it
EqnList . append ( eq )
center=integrate ( integrate ( f ( r , T ) * r , r , 0 , 1 / 2 ) , T , 0 , 2 * pi ) # i was j u s t ←-
t e s t i n g t o s e e what would happen w i t h o u t t h e c e n t e r
# print " center = " , center
EqnList . append ( center==M [ 0 , 0 ] )

# makes f c o n t i n u o u s & s m o o t h a l o n g t h e t a =0 by m a t c h i n g ←-
c o e f f i c i e n t s o f f ( r , 0 ) t o f ( r , 2 * p i ) −−s e v e r a l e q u a t i o n s

# m a t c h e s d e r i v a t i v e s by l o o p i n g
f o r d i n range ( deg + 1 ) :
new_eqns = matcher ( diff ( f , r , d ) , deg−d )
EqnList += new_eqns

p r i n t " # e q u a t i o n s = " , l e n ( EqnList )

p r i n t " # unknowns = " , l e n ( L )

# create a dictionary of solutions

i f solver== ' s a g e ' :
sol = solve ( EqnList , L , solution_dict=True )

i f sol = = [ ] or sol== ' No s o l u t i o n ' :

d = sol
r e t u r n "No s o l u t i o n "
else :
d=sol [ 0 ]
# s o l v e i n Sage
else :
 28

d = solvemma ( EqnList , L )
# s o l v e in Mathematica
p r i n t "MMA d = " , d

i f s t r ( d ) . replace ( ' ' , ' ' ) == " {} " :

r e t u r n "No s o l u t i o n "

fsubs = f . subs_expr ( d )
# s u b s t i t u t e s into f the assignments given in dictionary d
p r i n t " f s u b s = " , fsubs

freevars = l i s t ( Set ( fsubs . variables ( ) ) . difference ( Set ( [ r , T ] ) ) )

# substitute that solution into f
f o r freevar i n freevars :
var ( freevar )
p r i n t " # f r e e v a r s = " , l e n ( freevars )

r e t u r n fsubs , freevars

d e f search_cir_vol_deg ( M , N , solver= ' s a g e ' , p = 8 ) :

# loops through degree of polynomial u n t i l s o l u t i o n i s found
deg =3
sol= ' No s o l u t i o n '
w h i l e ( sol== ' No s o l u t i o n ' and deg <=15) :
# loops through degrees of polynomials u n t i l a s o l u t i o n i s found
poly=create_poly ( deg , r , T )
# c r e a t e s the polynomial of degree in loop
p r i n t " deg = " , deg
sol=search_mtx_cir_vol ( poly [ 0 ] , poly [ 1 ] , M , N , deg , solver )
# searches for a solution using that polynomial
print " "
deg+=1
r e t u r n sol

d e f solve_and_graph_cv ( M ) :
# f i n d s the s o l u t i o n to the system of equations then graphs i t

rows , cols=M . dimensions ( )

# c o l s i s always 8
N = cols * ( rows −1)+1

f=search_cir_vol_deg ( M , N , solver= ' s a g e ' , p = 8 )

print f [0]
# f u n c t i o n t h a t has a l r e a d y been s u b s t i t u t e d w i t h s o l u t i o n s
print f [1]
# l i s t of free variables
fvlist=f [ 1 ]
 29

d = {}
# c r e a t e s an e m p t y d i c t i o n a r y
f o r r i n fvlist :
# i t e r a t e s through the l i s t of f r e e v a ri a b l e s
d [ r ]=0
# s e t s each f r e e v a r i a b l e i n d i c t i o n a r y equal t o z e r o
print d
func=f [ 0 ] . subs_expr ( d )
# s u b s t i t u t e s dictionary assignments into the function

graph=polar_plot3d ( func )
show ( graph )
 30

A.5 Functions for line integrals over a circle

d e f search_mtx_cir_li ( f , L , M , N , deg , solver= ' s a g e ' ) :

# searches for a s o l u t i o n for l i n e i n t e g r a l s over a c i r c l e
k = ( N−1) / 8
f ( r , T ) =f
rlines=pi / 4
# r a d i a l l i n e s are at every m u l t i p l e of pi /4
split=pi / 8
# s p l i t s t h e d i f f e r e n c e b e t w e e n e a c h r a d i a l l i n e t o g e t p r o p e r ←-
interval
dist = 1 / 2
#how f a r up and down t h e r a d i a l l i n e t o go on e i t h e r s i d e o f e n t r y ←-
for a complete i n t e r v a l of 1
EqnList = [ ]
f o r i i n range ( 8 ) :
f o r j i n range ( 1 , k + 1 ) :
line_T=integrate ( f ( i , T ) * r , T , ( 2 * j−1) * pi / 8 , ( 2 * j + 1 ) * pi / 8 )
# i n t e g r a l s along arc segments
eq = line_T==M [ j , i ]
i f n o t eq :
# i f e q u a t i o n i s n o t " F a l s e " ( i n c o n s i s t e n t ) , t h e n append ←-
it
EqnList . append ( eq )

f o r m i n range ( 8 ) :
f o r p i n range ( 1 , k + 1 ) :
line_r=integrate ( f ( r , m * pi / 4 ) , r , ( 2 * p−1) / 2 , ( 2 * p + 1 ) / 2 )
# i n t e g r a l s down r a d i a l l i n e s
eq1 = line_r==M [ p , m ]
i f n o t eq1 :
# i f e q u a t i o n i s n o t " F a l s e " ( i n c o n s i s t e n t ) , t h e n append ←-
it
EqnList . append ( eq1 )

# makes f c o n t i n u o u s & s m o o t h a l o n g t h e t a =0 by m a t c h i n g ←-
c o e f f i c i e n t s o f f ( r , 0 ) t o f ( r , 2 * p i ) −−s e v e r a l e q u a t i o n s

# m a t c h e s d e r i v a t i v e s by l o o p i n g
f o r d i n range ( deg + 1 ) :
new_eqns = matcher ( diff ( f , r , d ) , deg−d )
EqnList += new_eqns

p r i n t " # e q u a t i o n s = " , l e n ( EqnList )

p r i n t " # unknowns = " , l e n ( L )

# create a dictionary of solutions

i f solver== ' s a g e ' :
sol = solve ( EqnList , L , solution_dict=True )

i f sol = = [ ] or sol== ' No s o l u t i o n ' :

 31

d = sol
r e t u r n "No s o l u t i o n "
else :
d=sol [ 0 ]
# s o l v e i n Sage

fsubs = f . subs_expr ( d )
# s u b s t i t u t e s into f the assignments given in dictionary d
p r i n t " f s u b s = " , fsubs

freevars = l i s t ( Set ( fsubs . variables ( ) ) . difference ( Set ( [ r , T ] ) ) )

# substitutes that solution into f
f o r freevar i n freevars :
var ( freevar )
p r i n t " # f r e e v a r s = " , l e n ( freevars )

r e t u r n fsubs , freevars

d e f search_cir_li_deg ( M , N , solver= ' s a g e ' , p = 8 ) :

# loops through degrees of polynomials u n t i l s o l u t i o n i s found
deg =3
sol= ' No s o l u t i o n '
w h i l e ( sol== ' No s o l u t i o n ' and deg <=15) :
# l o o p s t h r o u g h d e g r e e s o f p o l y n o m i a l s u n t i l a s o l u t i o n i s ←-
found
poly=create_poly ( deg , r , T )
# c r e a t e s the polynomial of degree in loop
p r i n t " deg = " , deg
sol=search_mtx_cir_li ( poly [ 0 ] , poly [ 1 ] , M , N , deg )
# searches for a solution using that polynomial
print " "
deg+=1
r e t u r n sol

d e f solve_and_graph_cli ( M ) :
# searches f o r a s o l u t i o n then graphs i t

rows , cols=M . dimensions ( )

# c o l s i s always 8
N = cols * ( rows −1)+1

f=search_cir_li_deg ( M , N , solver= ' s a g e ' , p = 8 )

print f [0]
# f u n c t i o n t h a t has a l r e a d y been s u b s t i t u t e d w i t h s o l u t i o n s
print f [1]
# l i s t of free variables
fvlist=f [ 1 ]
 32

d = {}
# c r e a t e s an e m p t y d i c t i o n a r y
f o r r i n fvlist :
# i t e r a t e s through the l i s t of f r e e v a ri a b l e s
d [ r ]=0
# s e t s each f r e e v a r i a b l e i n d i c t i o n a r y equal t o z e r o
print d
func=f [ 0 ] . subs_expr ( d )
# s u b s t i t u t e s dictionary assignments into the function

graph=polar_plot3d ( func )
show ( graph )
 33

Appendix B

Output

B.1 Solutions for volume over a square

B.1.1 3 × 3

(x, y) |--> 3*x^2*y - 6*x*y^2 + 45/11*y^3 - 9/2*x^2 + 9*x*y -

207/22*y^2 + 9/2*x + 247/44

B.1.2 4 × 4

(x, y) |--> -x^3*y^2 + 4*x^2*y^3 - 2*x*y^4 + 4*x^3*y - 18*x^2*y^2 +

4/3*y^4 - 8/3*x^3 + 18*x^2*y + 87/2*x*y^2 - 4*x^2 - 86*x*y - 31*y^2 +
212/5*x + 196/3*y - 251/10
 34

B.1.3 5 × 5

(x, y) |--> -25/36*x^4*y^3 + 25/36*x^3*y^4 - 5/12*x^2*y^5 -

3343/234630*y^7 + 125/24*x^4*y^2 - 75/274*x*y^5 - 845/72*x^4*y -
115/4*x^3*y^2 + 85/18*x^2*y^3 + 75785/4932*x*y^4 + 5959/62568*y^5 +
275/36*x^4 + 175/2*x^3*y + 50*x^2*y^2 - 28225/274*x*y^3 - 2375/36*x^3 -
1910/9*x^2*y + 51935/274*x*y^2 + 259319/21330*y^3 + 550/3*x^2 -
409297/10428*y^2 - 1367375/9864*x + 3109409/62568
 35

B.1.4 6 × 6

(x, y) |--> -9/200*x^5*y^5 + 23/480*x^4*y^6 - 14873/36575*x^3*y^7 +

63/80*x^5*y^4 + 1074399/209000*x^3*y^6 + 45732/36575*x^2*y^7 -
21/4*x^5*y^3 - 829/96*x^4*y^4 - 3716337/209000*x^3*y^5 -
406827/26125*x^2*y^6 + 17330637021/21943846400*y^8 + 1287/80*x^5*y^2 +
593/8*x^4*y^3 + 5528529/104500*x^2*y^5 - 4433/2400*x*y^6 -
913658751/78370880*y^7 - 2031/100*x^5*y - 19697/80*x^4*y^2 +
182385237827/3134835200*y^6 + 279/50*x^5 + 5125/16*x^4*y +
12953538/26125*x^3*y^2 + 136819/480*x*y^4 - 31489222203/313483520*y^5 -
50517/560*x^4 - 1614474777/1463000*x^3*y - 144034443/104500*x^2*y^2 -
87799/40*x*y^3 + 108707391/292600*x^3 + 4420728393/1463000*x^2*y +
2579027/400*x*y^2 - 307366119/292600*x^2 - 601757/80*x*y +
175199678597/685745200*y^2 + 6027117/2800*x - 63680388517/685745200
 36

B.1.5 7 × 7

(x, y) |--> 1813/86400*x^6*y^5 - 49/3200*x^5*y^6 + 7/640*x^4*y^7 -

1152216777509/2002350733344*x^2*y^9 - 12691/34560*x^6*y^4 -
343/2880*x^5*y^5 - 2459/60480*x^3*y^7 +
424229067376793/33372512222400*x^2*y^8 +
148154691517/4768032501504*x*y^9 + 61397/25920*x^6*y^3 +
9751/1920*x^5*y^4 + 4165/6912*x^4*y^5 - 4291/5760*x^3*y^6 -
20369238101933/192533724360*x^2*y^7 - 29780121959/7946720835840*y^9 -
26411/3840*x^6*y^2 - 169099/4320*x^5*y^3 - 208201/6912*x^4*y^4 +
215861/17280*x^3*y^5 + 159308290039896599/400470146668800*x^2*y^6 -
31886405085913/2980020313440*x*y^7 + 1109381/129600*x^6*y +
589939/4800*x^5*y^2 + 26591803/103680*x^4*y^3 -
487257997230893/834312805560*x^2*y^5 +
237598679630501/1862512695900*x*y^6 + 1555574697943/993340104480*y^7 -
34349/10800*x^6 - 691439/4320*x^5*y - 218099/256*x^4*y^2 -
28076159/51840*x^3*y^3 - 4705041986147531/7946720835840*x*y^5 -
200191392376031/9933401044800*y^6 + 442813/7200*x^5 +
11930555/10368*x^4*y + 13910939/5760*x^3*y^2 +
3195109801524151/2980020313440*x*y^4 +
 37

1321227231308929/13244534726400*y^5 - 982499/2160*x^4 -
665470871/181440*x^3*y + 532115415438729829/400470146668800*x^2*y^2 -
11815380238207/62083756530*y^4 + 1335821/864*x^3 -
7660815562425389/4966700522400*x*y^2 -
35462255394170437/33372512222400*x^2 + 965065055275153/3311133681600*y^2
+ 733442813831413/827783420400*x - 20742857285153/137963903400

B.1.6 8 × 8

(x, y) |--> -1/8100*x^7*y^6 + 2/2025*x^6*y^7 - 1/1350*x^5*y^8 -

591893/2736228600*x^3*y^10 + 2/675*x^7*y^5 - 49/2025*x^6*y^6 +
4643/343980*x^4*y^8 + 96203/147987000*x^2*y^10 +
602752/1867167315*x*y^11 - 8254/10854718875*y^12 - 43/1620*x^7*y^4 +
151/675*x^6*y^5 + 1007/3240*x^5*y^6 - 11330/51597*x^4*y^7 +
949931/30402540*x^3*y^8 - 117184775749/21173677352100*x*y^10 +
44/405*x^7*y^3 - 391/405*x^6*y^4 - 21/5*x^5*y^5 - 5767/45675*x^2*y^8 +
176707/986792625*y^10 - 133/675*x^7*y^2 + 3929/2025*x^6*y^3 +
76229/3240*x^5*y^4 + 49115/2457*x^4*y^5 + 51519374791/130701712050*x*y^8
+ 248/2025*x^7*y - 1048/675*x^6*y^2 - 2866/45*x^5*y^3 -
4513919/29484*x^4*y^4 - 3052031/86175*x^3*y^5 +
100533763/7047000*x^2*y^6 - 98564/5221125*y^8 - 32/4725*x^7 +
584/2025*x^6*y + 64253/810*x^5*y^2 + 3535033/7371*x^4*y^3 +
 38

180837258499/547245720*x^3*y^4 - 40696/435*x^2*y^5 -
31723186349969/1008270350100*x*y^6 + 32/675*x^6 - 1492/45*x^5*y -
33415345/51597*x^4*y^2 - 135977698/120645*x^3*y^3 +
540328139/2959740*x^2*y^4 + 73231869440/373433463*x*y^5 +
234033161/140970375*y^6 - 9488/4725*x^5 + 13348000/51597*x^4*y +
333216251497/228019050*x^3*y^2 - 430192025563381/1058683867605*x*y^4 -
82176/7735*y^5 + 73024/1755*x^4 - 432722908/1809675*x^3*y +
350369137/2055375*x^2*y^2 + 4434531292/197358525*y^4 -
21720226096/59719275*x^3 - 1841944/1305*x^2*y +
412299974784368/588157704225*x*y^2 + 1572028256/1507275*x^2 -
6886639168/172297125*y^2 - 964919754032128/2156578248825*x +
2488659577/42664050

B.1.7 9 × 9

(x, y) |--> -227/1881600*x^8*y^7 + 199/1881600*x^7*y^8 -

199/2419200*x^6*y^9 + 9285205465/24872695981056*x^4*y^11 -
49072489175075/42404830518403584*x^3*y^12 + 681/179200*x^8*y^6 +
3/5600*x^7*y^7 + 7591/7257600*x^5*y^9 - 3723341587/358913361920*x^4*y^10
+ 10581236792557/446178772289600*x^3*y^11 +
142426942329895336919/16687311162160730685440*x^2*y^12 -
 39

8663/179200*x^8*y^5 - 32129/403200*x^7*y^6 - 133/57600*x^6*y^7 +

1835877693527/22611541800960*x^4*y^9 -
7358155626026948459201/25599852351042030028800*x^2*y^11 -
3624407839388419/2824795051414237440*x*y^12 + 1629/5120*x^8*y^4 +
14307/11200*x^7*y^5 + 10137/12800*x^6*y^6 - 517/1890*x^5*y^7 -
6986886834794099/2141658106990080*x^3*y^9 +
872567149522349606407/232198207265687347200*x^2*y^10 +
11365017920571191/412713562706625600*x*y^11 +
28851392356349/562120670186432000*y^12 - 617929/537600*x^8*y^3 -
298777/32256*x^7*y^4 - 1399861/96000*x^6*y^5 - 407/768*x^5*y^6 +
267264373182929267/10708290534950400*x^3*y^8 -
1416007114813076222516731/61439645642500872069120*x^2*y^9 -
15378840102551/14053016754660800*y^11 + 198243/89600*x^8*y^2 +
56461/1600*x^7*y^3 + 58073/512*x^6*y^4 + 36899263/576000*x^5*y^5 -
12940358817453/358913361920*x^4*y^6 +
439301735172330309566333/7585141437345786675200*x^2*y^8 -
5307988814555484233/1155597975578551680*x*y^9 - 3733367/1881600*x^8*y -
6196249/88200*x^7*y^2 - 1090885337/2419200*x^6*y^3 - 15443/24*x^5*y^4 +
27633744458949373/113057709004800*x^4*y^5 -
2368693798153612073/3441950529091200*x^3*y^6 +
144968208239455679/3079131296505600*x*y^8 +
714991301836277/4818177173026560*y^9 + 174471/313600*x^8 +
180749/2800*x^7*y + 1478349/1600*x^6*y^2 + 1033119847/362880*x^5*y^3 +
8838518388608775809/3569430178316800*x^3*y^5 -
17533497893710548509/103178390676656400*x*y^7 -
71410513934275113/70265083773304000*y^8 - 726387/39200*x^7 -
249932917/288000*x^6*y - 1584071/256*x^5*y^2 -
10546904785383203/2261154180096*x^4*y^3 -
112616887825535401031430479/102399409404168120115200*x^2*y^5 +
1141351/4480*x^6 + 47569292/7875*x^5*y +
1420240053149937/89728340480*x^4*y^2 -
36597753308989660409/2677072633737600*x^3*y^3 +
975636269939802108487/825427125413251200*x*y^5 - 4108957/2240*x^5 -
844685518536534691/44415528537600*x^4*y +
3365475335461191441241/176686793826681600*x^3*y^2 +
486081833516550206567746963/38399778526563045043200*x^2*y^3 +
 40

13579345850284010311/56212067018643200*y^5 +
4575299031428299/690908221696*x^4 -
159102058718947727386048049/7822177107262842508800*x^2*y^2 -
9491700599768581929541/962998312982126400*x*y^3 -
26993452693425414971/16060590576755200*y^4 -
765655727700947739/118373959995200*x^3 +
20788867269632604283217/1412397525707118720*x*y^2 +
1888507345344640837/430194390448800*y^3 +
2695770302196901835507913/325924046135951771200*x^2 -
4774313147302122125/1124241340372864*y^2 -
24127512673251150201/4374206173573600*x +
55207158100046568199/45672304452647600

B.1.8 10 × 10

(x, y) |--> 3275/1316818944*x^9*y^9 - 4301/2090188800*x^8*y^10 +

1674671/152613828864*x^7*y^11 - 17275/146313216*x^9*y^8 -
352429451/1046494826496*x^7*y^10 - 17906593/174415804416*x^6*y^11 +
129565/54867456*x^9*y^7 + 436391/162570240*x^8*y^8 +
11010559801/3662731892736*x^7*y^9 + 6622782529/2092989652992*x^6*y^10 +
355399567/639524616192*x^5*y^11 - 811375/31352832*x^9*y^6 -
 41

579881/8128512*x^8*y^7 - 29662151237/1046494826496*x^6*y^9 -
24050692789/1395326435328*x^5*y^10 -
1312920035239/248236047360000*x^4*y^11 +
10647819290969322369078665/3024023442003306992172417564672*y^15 +
3538387/20901888*x^9*y^5 + 3228997/3686400*x^8*y^6 +
108270069017/697663217664*x^5*y^9 +
449688275290507/2383066054656000*x^4*y^10 +
14997350741/5036256352512*x^3*y^11 -
363017177554495/551921264298657792*x^2*y^12 -
22702646006971801469569858753/256257986492724681410758940295168*y^14 -
42335675/62705664*x^9*y^4 - 23541179/3870720*x^8*y^5 -
986891102497/348831608832*x^7*y^6 -
171710293396513/79435535155200*x^4*y^9 -
4179467076509/43952782712832*x^3*y^10 +
10955969885841391/482931106261325568*x^2*y^11 +
264575345/164602368*x^9*y^3 + 525085691/20901888*x^8*y^4 +
217720130975093/6104553154560*x^7*y^5 +
2670133986569/99666173952*x^6*y^6 +
19239087132697/21976391356416*x^3*y^9 -
36024477018358565/131708483525816064*x^2*y^10 +
24787790595188217373034023687/1152008930286974092256159072256*y^12 -
29227375/13716864*x^9*y^2 - 141938303/2322432*x^8*y^3 -
102234006951785/523247413248*x^7*y^4 -
588270903813841/1744158044160*x^6*y^5 -
9801224327969/66444115968*x^5*y^6 +
7695387275990911/35462292480000*x^4*y^7 +
177761461600399633/150523981172361216*x^2*y^9 - 2006567/914457600*x*y^10
- 1145831322016708204222399220515/4032031256004409322896556752896*y^11 +
12242647/9144576*x^9*y + 100139840497/1219276800*x^8*y^2 +
511091341712363/915682973184*x^7*y^3 +
1927701262218415/1046494826496*x^6*y^4 +
2154119238077749/1162772029440*x^5*y^5 -
243622609899002017/113479335936000*x^4*y^6 +
376347808205122469116315897392731/241921875360264559373793405173760*y^10
- 150085/571536*x^9 - 705647837/13547520*x^8*y -
6020663189057/7267325184*x^7*y^2 - 2740442799091157/523247413248*x^6*y^3
 42

- 7022606093035135/697663217664*x^5*y^4 +
3971339384423017829/397177675776000*x^4*y^5 -
1774629388967329/2092989652992*x^3*y^6 + 590957/3386880*x*y^8 -
40335485192230355465946457285775/12096093768013227968689670258688*y^9 +
6017855/580608*x^8 + 20413120348289/36336625920*x^7*y +
64197985825711/8305514496*x^6*y^2 + 309537666100801/10900987776*x^5*y^3
- 1193853963572428007/47661321093120*x^4*y^4 +
384404807676492821/36627318927360*x^3*y^5 + 20076785/1016064*x*y^7 -
1029436509683/8478546048*x^7 - 108374733628361/20763786240*x^6*y -
114909377548043/2768504832*x^5*y^2 +
2711368195090581079/82745349120000*x^4*y^3 -
1217923552164376415/21976391356416*x^3*y^4 -
2289115685900587501717/1756113113677547520*x^2*y^5 -
18942231577/43545600*x*y^6 + 790280753981/692126208*x^6 +
1058940190154653/38066941440*x^5*y -
87448083339333127/4728305664000*x^4*y^2 +
821664155838501461/5494097839104*x^3*y^3 +
10028817592660310697805/1053667868206528512*x^2*y^4 +
1707823739/483840*x*y^5 - 55050768482191/8882286336*x^5 -
4499309168727277/21801975552*x^3*y^2 -
876213034624285879861/37630995293090304*x^2*y^3 -
1858076107/142884*x*y^4 +
61366652985458005310678320800905/252001953500275582681034797056*y^5 +
127236590835157/27581783040*x^4 + 315495409400235913/2398217310720*x^3*y
+ 848146842824167090735/45993438691554816*x^2*y^2 +
5510586071/290304*x*y^3 -
2638523024112239103307265318552755/4158032232754547114237074151424*y^4 -
18049237172069611/559584039168*x^3 -
368686702420475398601/114983596728887040*x^2*y +
21413081563/50803200*x*y^2 + 43700859290646910393/4529656840834944*x^2 -
26894770183/1693440*x*y +
622043084185718750288547468045563/1365010581459826072855605150720*y^2 -
882502211/620928*x +
1547404439098024312557762802261/1276795611909701258623440192
 43

B.1.9 11 × 11
 44

B.1.10 12 × 12
 45

B.2 Solutions for line integrals over a square

B.2.1 3 × 3

(x, y) |--> -860/157437*x^3*y^3*(244*sqrt(2) - 917) +

464/52479*x*y^5*(244*sqrt(2) - 917) + 11200/549*x^6*(sqrt(2) - 2) +
11200/549*y^6*(sqrt(2) - 2) + 2/17493*x^2*y^3*(65788*sqrt(2) - 243803) +
2/17493*x^3*y^2*(45796*sqrt(2) - 173831) - 4/5831*x*y^4*(24142*sqrt(2) -
90587) + 8/21*x^4*y*(2*sqrt(2) - 7) + 2/157437*x*y^3*(2570684*sqrt(2) -
9707593) - 2/3843*y^4*(1434602*sqrt(2) - 2867557) -
4/3843*x^4*(718948*sqrt(2) - 1439543) - 2/52479*x^3*y*(213940*sqrt(2) -
788543) - 9/5831*x^2*y^2*(19708*sqrt(2) - 73493) +
2/254065*x^3*(311636947*sqrt(2) - 624369942) +
2/254065*y^3*(310918306*sqrt(2) - 621473601) +
1/17493*x^2*y*(595172*sqrt(2) - 2097403) -
1/1524390*y^2*(3727375024*sqrt(2) - 7458339959) -
1/762195*x^2*(1866758069*sqrt(2) - 3739916929) -
2/5831*x*y*(94552*sqrt(2) - 316631) + 1/32012190*y*(13902375286*sqrt(2)
- 27830909141) + 1/10670730*x*(4644397109*sqrt(2) - 9212597044) +
5773446719/21341460*sqrt(2) - 833150771/1524390
 46

B.2.2 4 × 4

B.2.3 5 × 5

B.2.4 6 × 6
 47

B.3 Solution for volume over a circle

-725110537030224393188290344252529/682670175962816821923333980771520*r^4\
/pi -
1816344198150318412011493278392/20186548011797388851819339955*T*r^4/pi^2
+
784777562187011393923473581019758/723258363848099816044466406675*T^3*r^4\
/pi^4 -
8693749420287366160917192792936448/130909763856506066704048419608175*T^7\
*r^4/pi^8 +
120245502615175637802726789773296594851695/39378278073775300852828110866\
00081830656*r^3/pi +
46866480915604496356452058366172260586/485172982593103988437386371568248\
85*T*r^3/pi^2 - 8895890936711/2271257415*T^2*r^3/pi^3 +
4564976394637654835962001121672366796544/3476626289631247917145249303447\
61745*T^3*r^3/pi^4 - 1292250084832/21631023*T^4*r^3/pi^5 +
7962491015168/108155115*T^5*r^3/pi^6 -
15579428098048/324465345*T^6*r^3/pi^7 +
148650410203621850645737184799546710491136/89895622631893696143327160560\
57410835*T^7*r^3/pi^8 - 1581960331264/757085805*T^8*r^3/pi^9 -
252794010545435109302353936859759252885166022541/58146252369824891073347\
71993657776433991754240*r^2/pi -
83806663267582804818647948035429055658969826/358204980878365516015521970\
10939791609575*T*r^2/pi^2 +
282028709619621697/67075105011957*T^2*r^2/pi^3 -
863167217738955847198502086380600718109893733/51336117149087301024544971\
848827524549855*T^3*r^2/pi^4 +
1354854223626451472/9582157858851*T^4*r^2/pi^5 -
9835312974643392896/79851315490425*T^5*r^2/pi^6 +
107579239095808342214537136424444825439663581184/18583674407969602970885\
27980927556388704751*T^7*r^2/pi^8 -
738345193419931648/22358368337319*T^8*r^2/pi^9 +
621332324736827392/111791841686595*T^9*r^2/pi^10 -
8275028289473842532761107168168028052624433913/1938208412327496369111590\
664552592144663918080*r/pi +
9364964187711456996997210979761055208445706/1194016602927885053385073233\
 48

6979930536525*T*r/pi^2 + 22153982806023791/3194052619617*T^2*r/pi^3 +
360275714269500057283738153439787682219298524/17112039049695767008181657\
282942508183285*T^3*r/pi^4 - 807316029646392560/3194052619617*T^4*r/pi^5
+ 6384318233538350336/26617105163475*T^5*r/pi^6 -
54037861150955871093355121495760844913719967744/442468438284990546925839\
995458941997310655*T^7*r/pi^8 +
372241275299790848/5323421032695*T^8*r/pi^9 -
443785954896904192/37263947228865*T^9*r/pi^10 +
5717346601461109617441191721311744/9694258138588628534123134485*T^6/pi^7
-
61310065455980249771496234746078464/37699892761177999854923300775*T^7/pi\
^8 +
2025025581758006609163118597259264/1077139793176514281569237165*T^8/pi^9
-
5044845067290146275387176962490368/4523987131341359982590796093*T^9/pi^1\
0 +
3247074696180182124139484563898368/9694258138588628534123134485*T^10/pi^\
11 -
654037165692554899635184666148864/16157096897647714223538557475*T^11/pi^\
12
49
 50

B.4 Solution for line integrals over a circle

The solution was degree 12, and the output was too long to present here.
 51

BIBLIOGRAPHY

[1] Dawn L. Anderson. Magic squares. http://illuminations.nctm.org/Lesson.aspx?id=

655. [Online; accessed 26-October-2014].

[2] W. W. Ball and H. S. Coxeter. Mathematical Recreations and Essays. The Macmillan Company,
New York, NY, 1947.

[3] Ron Larson, Robert P. Hostetler, and Bruce H. Edwards. Calculus, 8e. Houghton Mifflin
Company, Massachusetts, 2006.

[4] Mike Molony. Panmagic (pandiagonal) and magic square generation, testing. http://blog.
dreamshire.com/panmagic-%pandiagonal-magic-square-generation-testing/,
November 26 2013.

[5] William H. Richardson. Magic squares of even order (4n + 2). http://www.math.wichita.
edu/~richardson/mathematics/magic%20squares/even-ordermagicsquares.html.
[Online; accessed 29-October-2014].

[6] Eric W. Weisstein. Magic circle. http://mathworld.wolfram.com/MagicCircles.html.

[Online; accessed 29-October-2014].

[7] Eric W. Weisstein. Magic square. http://mathworld.wolfram.com/MagicSquare.html.

[Online; accessed 29-October-2014].