Peter Luschny

## PERFECT AND OPTIMAL RULERS

### Introduction

Rulers, which we can buy at the stationery shop, are edged with a regular sequence of marks and look somewhat like this

'-'-'-'-'-'-'

Rather dull. Recently, however, the manufacturers of rulers found out that they can reduce production costs by deleting some of the marks from the rulers without reducing their functionality. For example with the ruler

'-'---'--'

you can also measure all the distances you can measure with the standard ruler of equal length.

How many marks can be deleted, if we want to be able to measure all distances up to some length? In essence, it is this question which we are going to study here and which turns a monotonous ruler in an interesting object of combinatorics.

First of all, we want to be sure that we can measure all length of equal or shorter size then the length of the ruler. Clearly this is a sine qua non for a reasonable ruler. A ruler which has this feature will be called a complete ruler, because it has enough marks to fulfill this basic task. Rulers which can not (incomplete rulers) will not be considered here.

Next we focus on those complete rulers which have as few marks as possible. They will be called perfect rulers. They are complete rulers, but the deletion of a single mark would make them incomplete. There is a certain aesthetic impression associated with perfect rulers: they accomplish a task with a minimum of complexity.

Interestingly, perfect rulers with the same number of marks still differ in their usefulness. Let us consider the following three perfect rulers. All of them have 4 marks. (The numbers correspond to the positions of the marks on the ruler.)

[0, 1, 3, 4],   [0, 1, 3, 5],   [0, 1, 4, 6]

But with the third one you can measure 33% more lengths compared to the first one. So you would certainly choose the third one if you would like to buy one of these rulers.

But the decisive thing about the third ruler is that there does not exist any other ruler with 4 marks, which is complete and can measure more lengths. Thus the marks cover a maximal metering range. Rulers with this quality are called optimal rulers. Optimal rulers thus combine the aesthetics of a minimum of complexity with a maximum of usefulness.

Rulers with this feature are quite rare creatures. For example there are 52012 perfect rulers with 14 marks, but only 4 out of them are optimal rulers. The enumeration of complete, perfect and optimal rulers will be one of our major tasks.

Optimal rulers, which exist for every number of marks, have a very peculiar length. In fact their lengths are, in dependency on the number of marks,

0, 1, 3, 6, 9, 13, 17, 23, ...

But these numbers also count the maximal number of edges in a graceful graph on n nodes, linking the study of perfect rulers to graph theory, combinatorics and additive number theory.

Clearly we would like to know how to construct such optimal rulers. At first sight this does not seem to be an easy task, but there is an exciting, still unproved conjecture, that there is a blueprint for optimal rulers, which would allow a very simple construction. It says that a ruler with more then 14 marks is optimal only if it is a Wichmann-ruler, named after B. Wichmann, who introduced them in 1962 in A note on restricted difference bases.

### Rulers in a Nutshell

A ruler is a strict increasing finite sequence of marks, which are understood as nonnegative numbers. By convention the first mark is set to 0.

Rulers are counted with regard to their length L and the number of segments S (a segment is the space between two adjacent marks). A ruler with M marks has S = M - 1 segments.

A ruler is:

• COMPLETE, if all distances d with 1 ≤ d ≤ L can be measured with the ruler.

• PERFECT, if there is no complete ruler with the same length which possesses fewer marks.

• OPTIMAL, if there is no perfect ruler with the same number of marks which has a greater length.

## Rulers

In [3]:
%load_ext sage

In [4]:
def Partsum(T) :
return [add([T[j] for j in range(i)]) for i in (0..len(T))]


Given a list T the function partsum returns the list of partial sums of T. For example [1,2,3,4,5] maps on [0,1,3,6,10,15] (the triangular numbers A000217). Note that a '0' is prepended as the first sum is the empty sum.

In [5]:
Partsum([1,2,3,4,5])

Out[5]:
[0, 1, 3, 6, 10, 15]


Recall that a composition of a positive integer n is defined as a list of positive integers which sum to n. If we apply the summation of parts to compositions we get new objects. In the case n = 4 this is displayed in the plot below.

This is all we need to give the basic definition: A ruler is the list of partial sums of an integer-composition.

In [6]:
def Ruler(L, S) :
return map(Partsum, Compositions(L, length=S))

In [7]:
Ruler(6, 3)

Out[7]:
[[0, 4, 5, 6],
[0, 3, 5, 6],
[0, 3, 4, 6],
[0, 2, 5, 6],
[0, 2, 4, 6],
[0, 2, 3, 6],
[0, 1, 5, 6],
[0, 1, 4, 6],
[0, 1, 3, 6],
[0, 1, 2, 6]]


Let us clarify the terminology: a member of the ruler is called a mark, the length of a ruler L is the difference between the last and the first mark; so all the rulers in the last example have length 6. The second important quantity for classification is the number of segments S. A segment is the space between two marks. So the ruler has [0,  2,  5, 6] has length 6 and 3 segments. The number of marks is M = S + 1.

We denote the set of all rulers which have length L and S segments by R(L, S), where L >= 0 and S >= 0. The empty ruler is [] and the trivial ruler [0, 1, 2, ..., n] is characterized by L = S.

## Complete Rulers

A ruler is complete if all distances between 1 and the length of the ruler can be measured.

In [8]:
def isComplete(R) :
S = Set([])
L = len(R)-1
for i in range(L,0,-1) :
for j in (1..i) :
S = S.union(Set([R[i]-R[i-j]]))
return len(S) == R[L]

In [9]:
def CompleteRuler(L, S) :
return filter(isComplete, Ruler(L, S))


We give some examples.

In [10]:
list(CompleteRuler(6, 3))

Out[10]:
[[0, 2, 5, 6], [0, 1, 4, 6]]

In [11]:
for i in (0..5) :
for c in CompleteRuler(5, i) : print c

[0, 3, 4, 5]
[0, 2, 4, 5]
[0, 1, 3, 5]
[0, 1, 2, 5]
[0, 2, 3, 4, 5]
[0, 1, 3, 4, 5]
[0, 1, 2, 4, 5]
[0, 1, 2, 3, 5]
[0, 1, 2, 3, 4, 5]



If we put all complete rulers of the same length L into a bag we get the CompleteRulers(L).

In [13]:
def CompleteRulers(L) :
return sum([CompleteRuler(L, i) for i in (0..L)],[])

In [14]:
CompleteRulers(5)

Out[14]:
[[0, 3, 4, 5],
[0, 2, 4, 5],
[0, 1, 3, 5],
[0, 1, 2, 5],
[0, 2, 3, 4, 5],
[0, 1, 3, 4, 5],
[0, 1, 2, 4, 5],
[0, 1, 2, 3, 5],
[0, 1, 2, 3, 4, 5]]


Counting complete rulers:

In [18]:
for n in (0..11):
print [len(CompleteRuler(n,k)) for k in (0..n)]

[1]
[0, 1]
[0, 0, 1]
[0, 0, 2, 1]
[0, 0, 0, 3, 1]
[0, 0, 0, 4, 4, 1]
[0, 0, 0, 0, 8, 27, 20, 7, 1]
[0, 0, 0, 0, 4, 40, 48, 27, 8, 1]
[0, 0, 0, 0, 0, 38, 90, 75, 35, 9, 1]
[0, 0, 0, 0, 0, 30, 134, 166, 110, 44, 10, 1]



CompleteRulers(L, S)

 L\S 0 1 2 3 4 5 6 7 8 9 10 11 sum 0 1 1 1 0 1 1 2 0 0 1 1 3 0 0 2 1 3 4 0 0 0 3 1 4 5 0 0 0 4 4 1 9 6 0 0 0 2 9 5 1 17 7 0 0 0 0 12 14 6 1 33 8 0 0 0 0 8 27 20 7 1 63 9 0 0 0 0 4 40 48 27 8 1 128 10 0 0 0 0 0 38 90 75 35 9 1 248 11 0 0 0 0 0 30 134 166 110 44 10 1 495

This is sequence A103294. In the triangle the length of the rulers increase from top to bottom and the number of segments increase from left to right, entries above the main diagonal are 0. The number of complete rulers with length n are the row sums A103295. The most important sequence is the irregular diagonal to the right of the zero entries: the number of perfect rulers with length n A103300.

## Perfect Rulers

A ruler r is a perfect ruler if it is complete and no complete ruler with the same length possesses less marks.

In [19]:
def PerfectRulers(L) :
for i in (0..L) :
R = CompleteRuler(L, i)
if len(R) > 0 : return R
return []

In [20]:
PerfectRulers(5)

Out[20]:
[[0, 3, 4, 5], [0, 2, 4, 5], [0, 1, 3, 5], [0, 1, 2, 5]]

In [23]:
for i in (0..6) :
print [c for c in PerfectRulers(i)]

[[0]]
[[0, 1]]
[[0, 1, 2]]
[[0, 2, 3], [0, 1, 3]]
[[0, 2, 3, 4], [0, 1, 3, 4], [0, 1, 2, 4]]
[[0, 3, 4, 5], [0, 2, 4, 5], [0, 1, 3, 5], [0, 1, 2, 5]]
[[0, 2, 5, 6], [0, 1, 4, 6]]



A more appropriate name for len(x) would be cardinality(x).

In [24]:
len(PerfectRulers(10))

Out[24]:
38

In [25]:
len(CompleteRulers(10))

Out[25]:
248


## Representations of rulers

First we look at a symbolic representation which captures the intuitive meaning of a ruler.

In [26]:
def Ruler_AsRuler(R) :
l = 0
S = ""
for i in range(0, len(R)) :
while l < R[i] - i :
S = S + '-';
l = l + 1
S = S + '|';
return S

In [27]:
Ruler_AsRuler([0,1,4,10,16,18,21,23])

Out[27]:
'||--|-----|-----|-|--|-|'


If we substitute the symbols by zeros and ones we might get the people from the CS department on board.

In [28]:
def Ruler_AsBinaryString(R) :
l = 0
S = ""
for i in range(0, len(R)) :
while l < R[i] - i :
S = S + '0';
l = l + 1
S = S + '1';
return S

In [29]:
for c in PerfectRulers(5) :
print Ruler_AsBinaryString(c)

100111
101011
110101
111001


In [30]:
Ruler_AsBinaryString([0, 2, 5, 6])

Out[30]:
'1010011'

In [31]:
def Ruler_AsInteger(R) :
return int(Ruler_AsBinaryString(R),2)

In [32]:
Ruler_AsInteger([0, 2, 5, 6])

Out[32]:
83


Finally we look at the difference representation of a ruler.

In [33]:
def DiffRep(R):
L = []
for i in (1..len(R)-1) :
L.append(R[i] - R[i-1])
return L

In [34]:
DiffRep([0, 4, 5, 9])

Out[34]:
[4, 1, 4]

In [35]:
Partsum([4,1,4])

Out[35]:
[0, 4, 5, 9]


Clearly the difference representation of a ruler is nothing but the composition to which we referred in the definition of a ruler.

In [36]:
def Composition(R):
return DiffRep(R)


## Ordering rulers

In [38]:
for i in (0..5) :
print sorted([Ruler_AsInteger(c) for c in CompleteRulers(i)])

[1]
[3]
[7]
[11, 13, 15]
[23, 27, 29, 31]
[39, 43, 47, 53, 55, 57, 59, 61, 63]


In [39]:
for i in (0..7) :
print sorted([Ruler_AsInteger(c) for c in PerfectRulers(i)])

[1]
[3]
[7]
[11, 13]
[23, 27, 29]
[39, 43, 53, 57]
[83, 101]
[143, 151, 167, 171, 179, 203, 205, 211, 213, 229, 233, 241]


In [41]:
SPR = sorted(sum(([Ruler_AsInteger(c) for c in PerfectRulers(i)] for i in (0..7)),[]))
print SPR

[1, 3, 7, 11, 13, 23, 27, 29, 39, 43, 53, 57, 83, 101, 143, 151, 167, 171, 179, 203, 205, 211, 213, 229, 233, 241]



If we speak of the natural order of rulers we mean the order given by this procedure.

The output below shows the perfect rulers as binary strings in their natural order.

In [44]:
for c in SPR : print bin(c).lstrip("0b")

1
11
111
1011
1101
10111
11011
11101
100111
101011
110101
111001
1010011
1100101
10001111
10010111
10100111
10101011
10110011
11001011
11001101
11010011
11010101
11100101
11101001
11110001



Let's have a quick look at the Hemming distances of the rulers with equal length in this sequence.

In [46]:
def HammingDistance(s, t) :
assert len(s) == len(t)
return sum(a != b for a, b in zip(s, t))

In [47]:
T = sorted([Ruler_AsInteger(c) for c in PerfectRulers(12)])
b = "0000000000000"
L = []
for c in T :
a = bin(c).lstrip("0b")
L.append(HammingDistance(a, b))
b = a
L

Out[47]:
[6, 4, 6, 2, 2, 6, 4, 6, 6, 2, 6, 2, 6, 4]


## Optimal rulers

A ruler is optimal if it is perfect and no perfect ruler with the same number of segments has a greater length.

Each row and each column of the matrix CompleteRuler(L, S) has only finite many entries different from 0. Moreover, if S <= L and CompleteRuler(L, S) = 0 then CompleteRuler(K, S) = 0 for all L <= K and CompleteRuler(L, N) = 0 for all N <= S.

Thus we have the characterization of a nonempty ruler R:

R is perfect iff R is in CompleteRuler(L, S) and CompleteRuler(L, S-1) = 0;
R is optimal iff R is in CompleteRuler(L, S) and CompleteRuler(L+1, S) = 0.

The list below displays only one out of the pair [r, reverse(r)], i. e. disregards the mirror-symmetric rulers.

In [48]:
OptimalRulers = [
[0],
[0, 1],
[0, 1, 3],
[0, 1, 4, 6],
[0, 1, 4, 7, 9],
[0, 1, 2, 6, 9],
[0, 1, 6, 9, 11, 13],
[0, 1, 4, 5, 11, 13],
[0, 1, 2, 6, 10, 13],
[0, 1, 8, 11, 13, 15, 17],
[0, 1, 4, 10, 12, 15, 17],
[0, 1, 2,  8, 12, 15, 17],
[0, 1, 2,  8, 12, 14, 17],
[0, 1, 2,  6, 10, 14, 17],
[0, 1, 2,  3,  8, 13, 17],
[0, 1, 4, 10, 16, 18, 21, 23],
[0, 1, 2, 11, 15, 18, 21, 23],
[0, 1, 4, 10, 16, 22, 24, 27, 29],
[0, 1, 3,  6, 13, 20, 24, 28, 29],
[0, 1, 2, 14, 18, 21, 24, 27, 29],
[0, 1, 3, 6, 13, 20, 27, 31, 35, 36],
[0, 1, 3, 6, 13, 20, 27, 34, 38, 42, 43],
[0, 1, 3, 6, 13, 20, 27, 34, 41, 45, 49, 50],
[0, 1, 2, 3, 23, 28, 32, 36, 40, 44, 47, 50],
[0, 1, 5, 8, 12, 21, 30, 39, 48, 53, 54, 56, 58],
[0, 1, 3, 6, 17, 24, 27, 38, 45, 49, 53, 57, 58],
[0, 1, 3, 6, 17, 20, 27, 35, 45, 49, 53, 57, 58],
[0, 1, 2, 8, 15, 16, 26, 36, 46, 49, 53, 55, 58],
[0, 1, 2, 6,  8, 17, 26, 35, 44, 47, 54, 57, 58],
[0, 1, 2, 3, 27, 32, 36, 40, 44, 48, 52, 55, 58],
[0, 1, 2, 8, 15, 16, 26, 36, 46, 56, 59, 63, 65, 68],
[0, 1, 2, 5, 10, 15, 26, 37, 48, 54, 60, 66, 67, 68],
[0, 1, 2, 5, 10, 15, 26, 37, 48, 59, 65, 71, 77, 78, 79],
[0, 1, 2, 5, 10, 15, 26, 37, 48, 59, 70, 76, 82, 88, 89, 90],
[0, 1, 2, 5, 10, 15, 26, 37, 48, 59, 70, 81, 87, 93, 99, 100, 101]
]


There are 155050 perfect rulers with length in [0, 123], but only 74 optimal rulers in this range.

Optimal rulers are the most interesting creatures in our setup. Unfortunately no blueprint for these rulers is known. Therefore we turn to the Wichmann rulers which might provide an easy to compute substitute for optimal rulers. However, for convenience, we provide a function which returns for some small arguments optimal rulers based on known data.

In [49]:
OR = [[0,0],[1,1],[3,2],[6,3],[9,4],[13,5],[17,6],[23,7],[29,8],[36,9]]

def OptimalRulers(S):
assert(S in (0..9))
return CompleteRuler(OR[S][0],OR[S][1])

In [50]:
for i in (0..6) :
for c in OptimalRulers(i) :
print Ruler_AsBinaryString(c)

1
11
1011
1101
1010011
1100101
1001000111
1010010011
1100100101
1110001001
10010001000111
10100000110011
10101001000011
11000010010101
11001100000101
11100010001001
100010000100001111
100100010001000111
100101000100000111
101001000100000111
101001010000010011
101010100100000011
110000001001010101
110010000010100101
111000001000100101
111000001000101001
111000100010001001
111100001000010001



## Wichmann rulers

B. Wichmann. A note on restricted difference bases.
J. London Math. Soc. 38, 1962, 465-466

Let us look at two examples first.

In [51]:
R = [0, 1, 2, 8, 15, 16, 26, 36, 46, 56, 59, 63, 65, 68]
Composition(R)

Out[51]:
[1, 1, 6, 7, 1, 10, 10, 10, 10, 3, 4, 2, 3]

In [52]:
S = [0, 1, 2, 5, 10, 15, 26, 37, 48, 54, 60, 66, 67, 68]
Composition(S)

Out[52]:
[1, 1, 3, 5, 5, 11, 11, 11, 6, 6, 6, 1, 1]


Both rulers R and S are optimal. However the second one has a special structure which can be read of from his associated composition. To see this let us look at the type of the composition of S.

type(R) = [1*2, 6*1, 7*1, 1*1, 10*4, 3*1, 4*1, 2*1, 3*1]

type(S) = [1*2, 3*1, 5*2, 11*3, 6*3, 1*2]

The symbolic notation 'p*m' means 'p occures m times in immediate succession'. 'p' denotes the part and 'm' the multiplicity.

Def. The type of a ruler is the type of its associated composition.

Def. A ruler is of Wichmann type if its type has the form, for r >= 0, s >= 0,

W(r, s) = [1*r, r+1, (2r+1)*r, (4r+3)*s, (2r+2)*(r+1), 1*r]

We see that the second ruler is of Wichmann type: type(S) = W(2,3).

If a ruler R is of Wichmann type then the number of segments of R is S = 4r+s+2 and the length of R is L = 4r(r+s+2) + 3(s+1).

In [53]:
def flatten(I) :
""" Utility function unrelated to rulers """
L = []
for i in I :
if hasattr(i, "__iter__") :
L.extend(flatten(i))
else :
L.append(i)
return L

def Wichmann(r,s) :
C = [[1]*r,r+1,[2*r+1]*r,[4*r+3]*s,[2*r+2]*(r+1),[1]*r]
return Partsum(flatten(C))

In [54]:
Wichmann(2,3)

Out[54]:
[0, 1, 2, 5, 10, 15, 26, 37, 48, 54, 60, 66, 67, 68]


Or as a binary strings:

In [55]:
Ruler_AsBinaryString(Wichmann(2,3))

Out[55]:
'111001000010000100000000001000000000010000000000100000100000100000111'


Since the 1s are the markers of the ruler and the number of markers is M = S + 1 the binary string representation of a Wichmann ruler of type (r,s) contains exactly 4r+s+3 times the 1.

## Are optimal rulers of Wichmann type? The optimal ruler conjecture

Not every optimal ruler is a ruler of Wichmann type. In the last section we saw an example for such a ruler.

However it was conjectured:

All optimal rulers with more than 13 segments are either Wichmann rulers or the mirror images of Wichmann rulers.

If the conjecture is true than the sequence A004137 continues 168, 183, 198, 213, 232, 251, 270, 289, 308, 327,...

A more cautious conjecture is: In the set of optimal rulers with S segments where S > 12 there exists at least one ruler of Wichmann type.

In view of the small basis of evidence at the current state this conjecture should preferably be regarded as a challenge for research. In any case it is an interesting observation and an investigation of the relation between optimal rulers and Wichmann rulers appears to be a worthwhile endeavor.

It should also be noted that rulers with different Wichmann types can fall into the same class of optimal rulers as the example below shows.

In [56]:
print Wichmann(2,8)
print Wichmann(3,4)

[0, 1, 2, 5, 10, 15, 26, 37, 48, 59, 70, 81, 92, 103, 109, 115, 121, 122,123]
[0, 1, 2, 3, 7, 14, 21, 28, 43, 58, 73, 88, 96, 104, 112, 120, 121, 122, 123]



## Counting Rulers

Note that we use the number of segments S rather than the number of marks M as our primary id.  A larger table can be found here.

 Seg-ment Length Perfect Optimal Sum 0 0 1 1 1 1 1 1 2 2 1 3 3 2 3 4 3 9 5 4 6 2 4 7 12 24 8 8 9 4 5 10 38 88 11 30 12 14 13 6 OEIS A103297 A004137 A103300 A103299 A103301

The challenges now are:

• Extend this table!
• Generate perfect rulers fast!
• (Dis-)Prove the optimal ruler conjecture.

## Appendix

### Prime number rulers

A prime number ruler is a ruler such that the interior marks of the ruler are on prime number positions.

Prime number rulers were studied by Naoyuki Tamura. He proved that there exist only finite many prime number rulers and gave a complete list of the minimal prime number rulers which is shown below.

 segments length perfect optimal ruler 2 3 ★ / ★ [0, 2, 3] 3 4 ★ [0, 2, 3, 4] 3 6 ★ / ★ [0, 2, 5, 6] 4 8 ★ [0, 2, 3, 7, 8] 6 12 [0, 2, 3, 5, 7, 11, 12] 6 14 ★ [0, 2, 3, 5, 7, 13, 14] 6 14 ★ [0, 2, 3, 7, 11, 13, 14] 7 18 ★ [0, 2, 3, 5, 7, 11, 17, 18] 7 18 ★ [0, 2, 3, 5, 11, 13, 17, 18] 7 20 ★ [0, 2, 3, 7, 11, 17, 19, 20] 8 20 [0, 2, 3, 5, 7, 11, 13, 19, 20] 8 20 [0, 2, 3, 5, 11, 13, 17, 19, 20] 8 24 ★ [0, 2, 3, 5, 7, 11, 17, 23, 24] 9 24 [0, 2, 3, 5, 11, 13, 17, 19, 23, 24] 9 30 ★ [0, 2, 3, 5, 7, 11, 17, 23, 29, 30] 9 32 ★ [0, 2, 3, 7, 13, 17, 23, 29, 31, 32] 10 32 [0, 2, 3, 7, 11, 17, 19, 23, 29, 31, 32] 11 38 [0, 2, 3, 5, 7, 13, 19, 23, 29, 31, 37, 38] 12 44 [0, 2, 3, 5, 7, 11, 19, 23, 29, 31, 37, 43, 44] 12 44 [0, 2, 3, 5, 11, 13, 19, 23, 29, 37, 41, 43, 44] 12 44 [0, 2, 3, 5, 11, 17, 19, 23, 31, 37, 41, 43, 44] 12 44 [0, 2, 3, 5, 11, 19, 23, 29, 31, 37, 41, 43, 44] 12 44 [0, 2, 3, 7, 11, 13, 19, 23, 29, 37, 41, 43, 44] 12 44 [0, 2, 3, 7, 11, 17, 19, 23, 31, 37, 41, 43, 44] 12 44 [0, 2, 3, 7, 11, 19, 23, 29, 31, 37,41, 43, 44] 15 62 [0, 2, 3, 7, 13, 17, 23, 29, 31, 37, 43, 47, 53, 59, 61, 62] 15 62 [0, 2, 3, 7, 13, 19, 23, 29, 31, 37, 41, 47, 53, 59, 61, 62] 15 62 [0, 2, 3, 7, 13, 19, 23, 29, 31, 37, 43, 47, 53, 59, 61, 62]

From this table we infer that the only possible lengths of minimal prime number rulers are:

3, 4, 6, 8, 12, 14, 18, 20, 24, 30, 32, 38, 44, 62.


Compare sequence A227956 in the Online Encyclopedia of Integer Sequences. The only possible lengths of perfect prime number rulers are:

3, 4, 6, 8, 14, 18, 20, 24, 30, 32.


The only possible lengths of optimal prime number rulers are:

3, 6.


There are 102 prime number rulers in total, 28 of which are minimal prime number rulers, 12 perfect prime number rulers and 2 optimal prime number rulers.

A simple fact: The length of a prime number ruler is p + 1 for some prime number p.

Reference: Naoyuki Tamura, Complete List of Prime Number Rulers, Information Science and Technology Center, Kobe University, 2013.

A Sage implementation can be found at Prime Rulers.