SurrealNumbers
Intro
This is a package implementing some parts of the Surreal Number system invented by John Horton Conway, and explained by Knuth in "Surreal Numbers: How Two ExStudents Turned on to Pure Mathematics and Found Total Happiness."
It isn't intended to be useful, so much as educational, and an interesting test of Julia itseld. It was educational for me to write, in terms of learning Julia. I wanted a task that would be painful to do in Matlab but "easy" to do in Julia. It might also be helpful, I hope, for someone trying to learn about Surreal numbers. I certainly did learn a lot about them that would have been easily brushed to the side if I had only gone through the theorems.
Background: the Surreal Numbers
Surreal numbers aren't numbers as we are taught, but they have many of the same properties. The tricky thing is that they are defined recursively from the very start. The nice part is that they use only set operations.
The definition is as follows: a surreal number <
all of the
members of the right set (technically no member of the left set is
There is a starting point  we can always use empty sets  as so the
first surreal number (usually denoted zero, because it will turn out
to be the additive identity) is
Then on the "first day" a new generation of surreals can be created in terms of the initial case. On the second day we create a new generation and so on. Each has a meaning corresponding to traditional numbers in order to place a consistent interpretation on standard mathematical operators defined on the surreals.
Just to reiterate, the tricky thing is that everything is
recursive. Even comparitors like <
, and hence, we can't even know
if something is a valid surreal in its own right, but only through
recursively investigating its component sets.
In any case, this is not a survey or tutorial on surreal numbers. There are many out there better than I can write, e.g.,
 https://en.wikipedia.org/wiki/Surreal_number
 https://conservancy.umn.edu/bitstream/handle/11299/174778/Hebert_umn_0130M_15912.pdf?sequence=1
 https://math.stackexchange.com/questions/816540/proofofconwayssimplicityruleforsurrealnumbers
 http://www.math.harvard.edu/~knill/teaching/mathe320_2015_fall/blog15/surreal1.pdf
 https://www.cuttheknot.org/WhatIs/Infinity/SurrealNumbers.shtml
 https://www.whitman.edu/Documents/Academics/Mathematics/Grimm.pdf
 https://www.tondering.dk/download/sur16.pdf
 https://en.wikibooks.org/wiki/Surreal_Numbers_and_Games
 https://www.encyclopediaofmath.org/index.php/Surreal_numbers
 https://www.scmp.com/article/241601/gamegoinspiresnewsurrealnumberstheory
 https://www.edge.org/responsedetail/27156
 https://ncatlab.org/nlab/show/surreal+number#conwaystyle_definitions
The purpose here is simply to introduce the key reasons for implementing this in Julia  it enables users to define highperformance types of their own, and those types can be recursively defined.
I know you can do this in other languages, and in some cases also achieve high performance. But it's so easy in Julia.
An example
Let's start off with some small examples of using the package. After
installing the package we can create basic surreals using (i) the
constructor, (ii) conversion from another real number, or (iii) a couple of
special functions, e.g., zero
and one
. There are two constructors,
one includes an extra string we'll call the shorthand for the
surreal. It's use in printing out numbers. The second constructor, and
many other operators leave this blank. The empty set is designated by
∅, which you can get in Julia by typing \emptyset TAB
.
julia> using SurrealNumbers
julia> z = zero(SurrealShort)
julia> x1 = SurrealShort("1", [z], ∅)
julia> x2 = SurrealShort("2", [x1], ∅)
julia> x_something = SurrealShort([z, x1], ∅)
julia> x_half = convert(SurrealShort, 1 // 2)
These commands create several surreals, starting at zero, then 1, and
2, showing the recursive construction. Then a surreal whose value we
might not know (to start with), and then we convert the rational value
1/2 into a surreal. The type SurrealShort
is named because the
subset of surreals that are implemented here (those with finite
representations) is called the short surreals. There are several
aliases to this in the code: SurrealFinite
(which was what I started
with, and so is the fundamental type), and SurrealDyadic
(for
reasons described below).
Note that the output of these varies: in the case where a shorthand
was defined it just outputs (in bold) the shorthand, but otherwise it
will show the bracket format of the components. To see the bracket
format even when there is a shorthand defined use the command pf
, e.g.
julia> println(" x_half = ", x_half, " = ")
julia> pf(x_half)
julia> println()
We have most of the standard arithmetic operators (division has some restrictions  see below), so you can do things such as
julia> x2 + x_half
which will produce quite a long sequence. To see what it is, convert back to a real number,
julia> float( x2 + x_half )
or do a picture of the recursion using DOT (which needs to be separately installed from GraphViz), e.g.,
julia> file = "test_dot.dot"
julia> FID = open(file, "w")
julia> surreal2dot(FID, x2 + x_half)
julia> close(FID)
julia> run(`dot Tsvg O $file`)
with the following result. Note that the red notations were added manually. Each box is a surreal number, designated by the number at the top of the box, and its left and right sets are in the corresponding boxes below, with links shown in red for left, and blue for right. The recursion for each is show below, down to the point where each recursion stops at zero.
That seems like enough to get you started, so now a little about the implementation.
Icky bits  implementation details
While surreals use sets, and Julia has a Set type, implementation wise, sorted arrays of surreals were a better bet because that exchanges a once off cost when constructing a surreal with a reduced cost for many operations throughout the lifetime of the surreal. Also Julia creates a whole suite of array functions automagically when you define scalars so using arrays was a low pain way to get things working. Also, although the surreals use sets, i.e., order of the elements is not important, almost all texts do write these sets in order.
The surreals include all real numbers (and infinity and epsilon and
others). However, many of these require infinitely large sets
So the type Surreal
is an abstract type (a subtype of Real
) with
at present only one useful subtype SurrealShort
, where finite
(here) means that the representation is finite, not that the number
itself is finite.
This is a pretty rich set by itself, but it doesn't cover even the entire set of Rationals; only the dyadic rationals. So a few words on them seem in order (see below).
The plan is to add a SurrealTrans
for the transfinites and other
surreals with infinite representations. But that is a little harder to
do, so I leave that for the moment.
The actual type structure (minus the constructor  see the code) is just
struct SurrealShort <: Surreal
shorthand::String
L::Array{SurrealShort,1}
R::Array{SurrealShort,1}
h::UInt64
end
Note the addition of a shorthand
string, which isn't necessary, but
carries a little bit of extra information to make pretty printing and
hence debugging a little easier.
The variable h
contains the hash value of the surreal, which is
calculated the first time it is used. This avoids incurring the cost
of this (recursive) calculation every time the value is needed, which
is often as the hash is used to speed up other operations.
Dyadic numbers
The dyadic rational numbers are those that have a denominator that is an exact power of 2, that is, numbers of the form (note seems that GitHub doesn't let me use Mathjax to interpret maths).
\[ x = \frac{ n }{ 2^k } \]
It turns out that every dyadic has a finite representation as a surreal, and every finite surreal represents a dyadic. The easiest way to represent this is as a DAG (Directed Acyclic Graph) of the constructions as shown below:
It might seem a little limiting to be restricted to this set, but remember that floating point numbers are dyadics. They are a (binary) integer (the significand) multiplied by 2 to the power of a (binary) exponent (just called the exponent). Thus all (finite) floating point numbers have a finite surreal representation (though it may be very, very long).
I did start writing a type for dyadics, but I'm not sure what it would be useful for, so it didn't get far (I realise the irony here  the surreals themselves aren't exactly useful).
Converting a number to a surreal
The description of a surreal given above generates a form. The forms
satisfy the rules of arithmetic (
I think of this loosely by analogy to the rationals, e.g., we can have two different forms that have the same value, and we usually call the same "number"
2 // 4 = 1 // 2
However, it seems important to distinguish equality (where two forms are identical) from equivalence (where they have the same value). In programming terms two "things" are equal when they are the same, not some airyfairy notion of equivalence, so equality and equivalence have different meanings and uses. Hence, here we have the relation ≅ defined to test equivalence separate from ==.
BTW, here we hit one of the weirdnesses of Julia; 99% of the time,
you can redefine operators and comparators to do whatever you like on
your new type. But you can't redefine
 https://docs.julialang.org/en/stable/devdocs/functions/#Builtins1
 https://discourse.julialang.org/t/overloadforcustomtype/4898
Anyway, I have defined congruence ≅
to do the same
thing, check for equivalence (equality of value). However, as there
are many possible surreal forms we could create to represent any given
real number, we have to chose one. Call that the canonical form. We
could define it in several ways, but the standard is

zero =>
$0 = {  }$ 
integers n =>
$n+1 = { n  }$ 
dyadic fraction => $ x = \frac{ n }{ 2^k } $ (for n odd) becomes
[ { x  1/2^{k}  x + 1/2^{k} } ]

negative number => use the identity that
$x = {X_R  X_L}$
This is the set of conversions implemented, and it includes all dyadics, and hence floats. Note the recursion inherent in the construction/conversion. So some care (particularly) with floats should be taken not to create a variable that exhausts the stack. The current code isn't very clever in checking for this another place some work is needed.
Examples:
julia> pf( convert(SurrealShort, 1//2) )
{ 0  1 }
julia> pf( convert(SurrealShort, 3//4) )
{ 1//2  1 }
julia> pf( convert(SurrealShort, 11//8) )
{ 3//2  5//4 }
The pf
function used here is a "printinfull", which prints the
left and right sets of the surreal, not the real equivalent
shorthand.
Converting a surreal back to a real
The conversion of canonical forms is relatively easy, but remember that noncanonical forms are possible, and can be quite counterintuitive.
For instance, naively, you might expect that the form ${ 3 
17 }$ could be mapped to the mean of the two sets, i.e., 10.
However,
Note, often in texts it is written
Conway defines a surreal
The secret of the conversion is again to use recursion, but that is not quite enough in this case. We use several tricks along the way:
 If
$x$ is equivalent to a known surreal such as 0 or 1, we convert directly  If
$x$ is negative, we use the identity that$x = { X_R  X_L }$  And, most importantly, if
$0 < x < 1$ we know$x$ will be the "simplest" dyadic rational number such that$x_L < x < x_R$ . In the interval$(0,1)$ simplest means having the denominator with the lowest power, i.e., in order of preference we would like the denomator to the$1,2,3, \ldots$ .
We can find the latter though a simple modification of the standard binary search a simplified version of which is shown below.
a = 0; b = 1
while true
d = (a + b) / 2
if xl < c < xr
return d
elseif c <= xl
a = d
elseif c >= xr
b = d
end
end
The first and last rules allow us to convert any number $x \in
[0,1]$. To convert numbers above this range, i.e.,
The result is not the world's most beautiful code  I'm sure it can be improved, but there are so many other inefficiency's here, I am not sure it warrants it.
Implementations of operators and standard functions
Most of the operators follow standard surreal definitions and defining them in Julia is easy. They are all recursive, as you might guess, and so very inefficient  I wouldn't want to do demanding computations this way, but they are easy to understand, for instance
+(x::SurrealShort, y::SurrealShort) = SurrealShort([x.X_L .+ y; x .+ y.X_L], [x.X_R .+ y; x .+ y.X_R] )
Notice that we are exploiting here Julia's natural extension of
operators from scalar to vectors (this is one of the reasons that
using arrays instead of sets for
x .+ [y_1, y_2]
once addition is defined on the scalar surreal type, without any additional definitions, and this is particular appealling here as the scalar addition operator is recursively defined in terms of the vector+scalar addition .+
.
Mutiplication was a bit of a bugbear to get right because I
misinterpreted the definition. The definition of multiplication of
X_L y + x Y_L  X_L Y_L
I assumed the +
operator in this definition was the standard
(surreal) addition of sets of surreals, but it isn't. Instead the
entire expression should be interpreted as "for all pairs of elements
from X_L and Y_L perform the above computation, and form a set from
the results." That's easy enough to do once you work out what you are
trying to do. There is another underlying problem which is that
canonical forms all have sets with 1 or 0 elements, and my broken
multiplication worked for them. It was only when I created products of
noncanonical forms that I saw problems, and then only when I had
cleaned out another misunderstanding from the code.
The one interesting thing to note is that even simple computations often generate noncanonical forms. Part of the aim of this package was to let people experiment and see such things. It is otherwise far to laborious to calculate anything but the very simplest cases (none of the above texts do any but the simplest).
One really useful example (for me to understand the surreals) is the following, which was obtained by multiplying the canonical form of 2 with itself. You might naively think that operators applied to canonical forms resulted in canonical forms, but this is not the case. The form is quite complex. Moreover, superficially it recurse into itself, i.e., the surreal number 4 appears twice in the tree. However, the form is defined in terms of the tree, so the dependent "4" is not the same surreal as the parent "4". That is, they are not equal, but they are equivalent. Thus the parent "4" is defined in terms of simpler forms (even though one of these is equivalent).
The other pieces of the toolkit are the standard things you expect to be able to do with numbers, e.g. round, sign, isinteger, ... I haven't implemented all of these, but hopefully enough that any others can be easily added.
There are two approaches: one is to use intrinsic surreal arithmetic,
e.g. sign
and abs
are implemented using native surreal arithmetic and operators. The result is that they looks almost exactly like it would for any other number.
sign(x::SurrealShort) = x<zero(x) ? one(x) : x>zero(x) ? one(x) : zero(x)
abs(x::SurrealShort) = x<zero(x) ? x : x
Actually, I don't even have to define abs
as I get this for free
because Julia has a similar operation defined abs(x::Real)
for all
real numbers. This is a great feature of Julia's type hierarchy and
multiple dispatch function selection.
The other approach to function definition is somewhat of a cheat. It
involves converting the number to a real, and then using the
appropriate operation on that field. I have tried to avoid that
approach when possible. At the moment, only the simplify
function
uses this approach. But if I am going to do more complex math
functions, e.g., logs or trig functions, I think I will have to take
this approach. The mathematical definitions of such in pure surreal
terms are obscure.
The cheat is used as part of one of the pieces of this that is hard to
implement, namely division. Division is welldefined on the surreals,
but, as even simple divisions such as 1/3 result in nondyadic rational
numbers, and hence have a infinite representations as a surreal, the
class of "finite" surreals defined here is not closed under
division. Eventually, this can be solved by outputting a nonfinite
surreal, but for the moment I have only implemented division when
dividing by
The show
command is designed to show the full structure unless there
is a "shorthand" defined for a surreal (most of the simple conversions
will set this up). This aids in viewing the surreals succintly, but
sometimes we want to see deeper. In the case where shorthand is
defined we can use the command pf
to see deeper, but it will stop at
the first layer below with a shorthand.
To see the full recursion I have implemented output of a surreal form
into the DOT syntax from GraphViz. The
function surreal2dot
can output a .dot
file, and then this can be
parsed (assuming you have GraphVis installed) by commands such as
using SurrealNumbers
s2 = convert(SurrealShort, 3//4)
file = "test_dot_s2.dot"
FID = open(file, "w")
surreal2dot(FID, s2)
close(FID)
run(`dot Tsvg O $file`)
x5 = SurrealShort( convert.(SurrealShort, [1, 1//2, 0]), [one(SurrealShort)] )
file = "test_dot_x5.dot"
FID = open(file, "w")
surreal2dot(FID, x5)
close(FID)
run(`dot Tsvg O $file`)
Which produces the figures like those below, illustrating the recursive definition of the two surreal numbers given.
The thing to note about these tree representations is that they are inefficient,
the same surreal forms are repeated. The code includes a function surreal2dag
that you use in exactly the same way to generate a DAG (Directed Acyclic Graph)
representation of the surreal as you see below.
See the code test/test_dot.jl
for other examples.
Uniquely surreal functions
There are two pieces that are unique to surreals:
Generation: the generation of a surreal is 1 + the generation of
the surreals used to construct it. Again this is easy to implement
recursively. Generation comes from Knuth's story where it's called
the "birth day" of the surreal. Generation can be thought of as a
function
\[ g(x) = \sup_{y \in X_L \cup X_R} \big[ g(y) + 1 \big] \]
Implicitly, generation is the tree depth, if viewed as the recursive tree shown above.
Generation is important as
Note that within an equivalence class of forms, not all forms have the same generation.
Canonicalise: convert a surreal into its equivalent canonical form. The easiest way to implement this was to use a similar cheat to that above, i.e., convert to a real, and then convert back to the equivalent surreal in canonical form.
Arrays of Surreals
We mentioned arrays of surreals earlier  they are used in the implementation instead of sets. One of the nice things about Julia us that you get many of the array operators and functions for free when you create scalar operators. So, for instance, you can write
convert.(SurrealShort, [1, 0, 1, 2] )
which will broadcasts the convert function across the array of integers to create an array containing the respective surreals. Likewise, we can use comprehensions to construct arrays, e.g.,
[ convert(SurrealShort, i) for i=1:2 ]
Or we could construct and iterator for the same thing (once floor
and some promotion rules are defined), e.g.., the iterator from 1 to 2 is
convert(SurrealShort, 1 ):convert(SurrealShort, 2)
However, in order to use arrays as sets we need, in the constructor
for a surreal to reduce the "set" to a sorted array containing unique
elements. Julia has nice sort and unique functions, but they rely on
hash functions, so we have to add a hash. These need to be recursive,
and work for arrays of surreals as well. The hash function help says
that we should implement such for new types such that isequal(x,y)
implies hash(x)==hash(y)
, with a second argument to be mixed in the
results. This is linked to the idea that we should separate the == and
the congruent comparisons  my hash is based on two forms being
equal, not being equivalent. In any case, Julia's syntax is again
simple and concise for specifying the hashes need (believe it or not
its another recursive function), but I need the hash to be really
fast. I use it to speed other pieces of code up, so I can't afford to
recurse every time. So instead, the hash value is calculated the first
time it is needed, and then stored in the surreal data structure
itself, so that we only need call the function once at most. We could
have defined the hash in the constructor, but it turns out that many
temporary surreals are created in multiplication (in particular), and
so avoiding the overhead when it isn't needed is just fine.
Other comments
Don't bother to tell me that this is horribly inefficient. It will never be anything but. Surreals were not created with numerical computing in mind. They are about as inefficient a way to do calculations as I can reasonably think of (excluding the addition of nullops everywhere).
I'm still learning Julia, and part of this is absorbing the style guide https://docs.julialang.org/en/stable/manual/styleguide/ so the code is still morphing to make it more efficient, and more stylish.
And note that this isn't really exploring an entire chunk of the surreals, i.e., the transfinite numbers that can be represented this way. I'll get to that one day, time gods willing.
There are other implementations of the surreals. For instance
 Coq https://dl.acm.org/citation.cfm?id=2150520
 Coq https://github.com/pirapira/surreal
 Haskell https://github.com/Lacaranian/surreal (only integers)
 Haskell https://github.com/serialhex/SurrealNumbers (not quite working)
 Haskell https://github.com/elfakyn/Haskellsurreals (only integers)
 Javascript http://www.iwriteiam.nl/SurrealNumbers.html
 Mathematica http://demonstrations.wolfram.com/GeneratingTheSurrealNumbers/ (only generation)
 Lua https://github.com/Ryan1729/luaSurrealnumbers?files=1
 Python https://github.com/codeinthehole/pythonsurreal (not much implemented)
 Python https://github.com/314eter/surrealnumbers
 Ruby http://raganwald.com/2009/03/07/eleganceandthesurreals.html
And some of these languages might be more appropriate in some ways for this task. But I wanted to learn Julia, and see how far I could take it here. Moreover, most of these are at least as incomplete as the code here. For instance, none (as far as I am aware) implement division (the Julia toolkit here has a very limited version of division included).
More information about Surreals
More Examples
Some examples of code are included here to make a little of this more real.
Creating surreals
Let's start by creating some surreals
x0 = SurrealShort("0", ∅, ∅)
x1 = SurrealShort("1", [x0], ∅)
x11 = SurrealShort("1", ∅, [x0])
x4 = SurrealShort( [x11], [x1,x0])
x5 = SurrealShort( [x11,x0,x4], [x1])
x21 = SurrealShort("2", [x1], ∅)
x22 = SurrealShort( [x0, x1], ∅)
x23 = SurrealShort( [x11,x1], ∅)
x24 = SurrealShort( [x11,x0,x1], ∅)
x25 = [x11,x0,x1] ≀ ∅
z = zero(x1)
z = zero(SurrealShort)
z = one(x1)
s1 = convert(SurrealShort, 1//2)
s2 = convert(SurrealShort, 3//4)
You can check the values using float
to convert them back to
standard floating point numbers. Note that ∅ is shorthand for an empty
array of surreal numbers, which is quite helpful in many places.
Try printing them out in various forms:
print("s1 = ", s1, " = " )
pf(s1)
println()
Operations on surreals and arrays of surreals
Test out some operations on the above variables with expected results
# comparisons
x0 <= x1
x0 <= x0
x1 >= x11
!( x11 ≅ x1 )
x21 ≅ x22 ≅ x23 ≅ x24
# addition and subtraction
x1 == x11
 x0 == x0
 x4 == x4
x1 + x11 ≅ x0
x1  x1 ≅ x0
x1 + x0 == x1
# multiplication
x1*x1 ≅ x1
x0*x1 ≅ x0
x4*x0 ≅ x0
convert(SurrealShort,2)*convert(SurrealShort,2) ≅ convert(SurrealShort,4)
# division
x11/one(x11) ≅ x1
x22 / x22 ≅ one(x22)
x22/x1 ≅ x22
x4/x1 ≅ x4
float(x1/x4) == float(x1)/float(x4)
convert(SurrealShort, 6)/ convert(SurrealShort, 3) ≅ convert(SurrealShort, 2)
Simple functions
Many of the simple numerical functions will work, but not the more
advanced ones such as log
. The two specific to surreals are
generation
which returns the generation or birth day of a surreal,
and canonicalise
which reduces it to the canonical form of the
surreal.
generation( zero(SurrealShort) ) == 0
canonicalise( convert(SurrealShort,2)*convert(SurrealShort,2) ) == convert(SurrealShort,4)
Converting back to standard real
Examples of how to convert a number back to rationals or floats.
convert(Rational, convert(SurrealShort, 1//8)) == 1//8
float( x4 ) == 0.5
convert( Rational, SurrealShort( [ 7//16 ], [ 15//16 ] ) ) == 1//2
convert(String, x4) == "1/2"
And of how conversion is automatically applied by promotion rules
x0 <= 1.0
x1 == 1//1
1.0 + x1 == 2.0
canonicalise(1//2 + x1) == 1.5
1//2 + x1 ≅ 1.5
Final notes
The total implementation here is a little more than 1,000 lines of code. No doubt an expert in Julia could make it a good deal tighter  I have concentrated on making the code easy (for me) to understand rather than super concise.
This little project would have been a lot easier if I knew more Julia, or more about the surreals. Trying to build something to learn about two moving parts at once wasn't a brilliant idea :)
But I was correct in thinking that (at least for me) this would have been almost impossible to build in Matlab. And at the least, it would have required a good deal more work without all the automagical pieces of Julia helping.