Skip to content

Algorithms

Villanova Summer REU

Tags: math, research, website

So it was a fine summer morning. Our research team was about to start our project “Counting Arithmetical Structures on Cycles with a Doubled Edge”. From the description of the project it is immediately obvious that this kind of counting can be done by brute force, and me being the only computer scientist in my group of mathematicians I start coding up a couple algorithms for this.

Given the fact this article’s thesis isn’t arithmetical structures or combinatorics, you can read more about those here:

We’ll start simple.

What is a graph?

Here we see a couple graphs, some of them have shapes that are so recognizable that they’re given names.

In this case, we were looking at doubled edge cycles:

It’s literally just a cycle graph and you add an extra edge between two vertices. Additionally, the set of cycle graphs is usually denoted as (C(n)) where (n) is the amount of vertices. Then, the set of doubled edge cycles is denoted as (C^\sim(n)).

What is an arithmetical structure?

  1. You take a graph and you put integers on each vertex.
  2. If the number at the vertex divides the sum of its neigboring vertex, then it’s an arithmetical structure.

Technically, you also need an extra condition but this is the gist.

Let’s see an example.

We have a graph with four vertices: 1, 1, 2, 3. Let’s verify if each vertex divides the sum of its neighbors:

We checked the divisibility condition for all the vertices and we have confirmed:

This is an arithmetical structure

Remarks

Armed with this knowledge, it is pretty obvious that we can just brute force all the possible combination of numbers on each vertex and just check if the divisibility conditions are met.

Let’s code!

For example, for (C^\sim(3)) we have a weird pizza shaped thing:

If we label the top two vertices “a” and “b” and the bottom vertex “c” we can say that a valid arithmetical structure meets the following criteria (in Python):

# Their `gcd` must be 1.
assert gcd([a, b, c]) == 1
# The sum of `b + b + c` or `2 * b + c` must be divisible by `a`.
assert (2 * b + c) % a == 0
# The sum of `a + a + c` or `2 * a + c` must be divisible by `b`.
assert (2 * a + c) % b == 0
# The sum of `a + b` must be divisible by `c`.
assert (a + b) % c == 0

Three of these properties have to do with divisibility, so let’s group them and pop them into a function:

def is_valid(a: int, b: int, c: int) -> bool:
    divisible = (2 * b + c) % a == 0 and (2 * a + c) % b == 0 and (a + b) % c == 0
    return gcd([a, b, c]) == 1 and divisible

and iterate over all possible values:

for a in range(1, 1000):
    for b in range(1, 1000):
        for c in range(1, 1000):
            if is_valid(a, b, c):
                print(a, b, c)

This will print all of the valid arithmetical structures on (C^\sim(3)) with vertex values under 1000.

Now if we wanted to do this for (C^\sim(6)) . . .

def is_valid(a: int, b: int, c: int, d: int, e: int, f: int) -> bool:
	divisible = (2 * b + f) % a == 0
	divisible &= (2 * a + c) % b == 0
	divisible &= (b + d) % c == 0
	divisible &= (c + e) % d == 0
	divisible &= (d + f) % e == 0
	divisible &= (e + a) % f == 0
	return gcd([a, b, c, d, e, f]) and divisible

for a in range(1, 1000):
	for b in range(1, 1000):
		for c in range(1, 1000):
			for d in range(1, 1000):
				for e in range(1, 1000):
					for f in range(1, 1000):
						if is_valid(a, b, c, d, e, f):
							print(a, b, c, d, e, f)

Yup, this is it . . . Looks great doesn’t it ?

This doesn’t actually run on the compute server we host our jupyter notebooks on so I had to transpile this to C++ to even be able to compute anything.

auto is_valid(std::vector<uint64_t> arr) -> bool {
    for (uint64_t i = 1; i < arr.size() - 2; ++i) {
        if ((arr.at(i) + arr.at(i + 2)) % arr.at(i + 1) == 0)
            return false;
    }

    return (arr.at(arr.size() - 2) + arr.at(0)) % arr.at(arr.size() - 1) == 0
        && (2 * arr.at(1) + arr.at(arr.size() - 1)) % arr.at(0) == 0
        && (2 * arr.at(0) + arr.at(2)) % arr.at(1) == 0;
}

auto search() {
    for (uint64_t a = 1; i < 1000; ++i)
    for (uint64_t b = 1; i < 1000; ++i)
    for (uint64_t c = 1; i < 1000; ++i)
    for (uint64_t d = 1; i < 1000; ++i)
    for (uint64_t e = 1; i < 1000; ++i)
    for (uint64_t f = 1; i < 1000; ++i) {
        auto const v = std::vector{a, b, c, d, e, f};
        if (is_valid(v) and
            std::gcd(a, std::gcd(b, std::gcd(c, std::gcd(d, std::gcd(e, f))))) == 1)
            std::print("{}\n", v);
    }
}

It has a complexity of O(n^m) where n is the upper bound and m is the degree of the doubled edge cycle. Additionally, I will save you the horrors of looking at the C++ version of all of this . . . (for now)

The C++ version took 5 days to output all of the structures for (C^\sim(8)) . . .

Optimizations

I won’t get into the details, but it turns out that we only need to vary the two vertices that share the doubled edge. The rest of the vertices on the graph can be constructed from those two. For simplicity let’s just say that we have a function F that will return the rest of the doubled edge cycle.

The function goes a little like this (in C++):

auto F(uint64_t x, uint64_t y) -> std::vector<uint64_t> {
    auto lst = std::vector{x, y};
    lst.reserve(x + y);
    while (lst.back() > 1)
        lst.push_back(mod(-lst.at(lst.size() - 2u), lst.back()));

    return lst;
}

C++ modulus doesn’t compute negative modulus correctly so I had to implement that myself with:

auto mod(auto a, auto m) {
    return ((a % m) + m) % m;
}

This means the search space went from n^m to just n^2! Now we need to iterate over every pair of two numbers from 1 to the upper bound and check if we can merge them at the start.

auto search(uint64_t upper_bound) -> std::unordered_set<std::vector<uint64_t>> {
    std::unordered_set<std::vector<uint64_t>> ct;

    for (uint64_t a = 1; a <= upper_bound; ++a) {
        for (uint64_t b = 1; b <= upper_bound; ++b) {
            if (std::gcd(a, b) != 1) continue;

            auto forwards = F(2 * a, b);
            auto backwards = F(2 * b, a);
        }
    }

}

First of all, we want to find the total count for