by Tom Switzer on Jul 07, 2013

technical

In this post I’d like to introduce you to what I have been calling *generic
numeric programming*.

What do we mean by generic numeric programming? Let’s take a simple example; we
want to add 2 numbers together. However, we don’t want to restrict ourselves to
a particular type, like `Int`

or `Double`

, instead we just want to work with
some *generic* type `A`

that can be added. For instance:

```
def add[A](x: A, y: A): A = x + y
```

Of course, this won’t compile since `A`

has no method `+`

. What we are really
saying is that we want `A`

to be some type that *behaves like a number*. The
usual OO way to achieve this is by creating an interface that defines our
desired behaviour. This is less than ideal, but if we were to go this route,
our `add`

function might look like this:

```
trait Addable[A] { self: A =>
def +(that: A): A
}
def add[A <: Addable[A]](x: A, y: A): A = x + y
```

We’ve created an interface that defines our `+`

method, and then bound our type
parameter `A`

to subsume this interface. The main problem with this is that we
can’t directly use types out of our control, like those that come in the
standard library (ie. `Int`

, `Long`

, `Double`

, `BigInt`

, etc). The only option
would be to wrap these types, which means extra allocations and either explicit
or implicit conversions, neither of which are good options.

A better approach is to use type classes. A discussion on type classes is out
of the scope of this post, but they let us express that the type `A`

must have
some desired behaviour, without inheritence. Using the type class pattern, we
could write something like this:

```
trait Addable[A] {
// Both arguments must be provided. Addable works with the type A, but
// does not extend it.
def plus(x: A, y: A): A
}
// This class adds the + operator to any type A that is Addable,
// by delegating to that Addable's `plus` method.
implicit class AddableOps[A](lhs: A)(implicit ev: Addable[A]) {
def +(rhs: A): A = ev.plus(lhs, rhs)
}
// We use a context bound to require that A has an Addable instance.
def add[A: Addable](x: A, y: A): A = x + y
```

We can then easily add implementations for any numeric type, regardless if we control it or not, or even if it is a primitive type:

```
implicit object IntIsAddable extends Addable[Int] {
def plus(x: Int, y: Int): Int = x + y
}
add(5, 4)
```

This is, more or less, the approach Spire takes.

Why be generic? The flippant answer I could give is: why not? I do hope that after reading this, that is an acceptable answer to you, but I know that’s not what you came here for.

The first reason is the obvious one; sometimes you want to run the same
algorithm, but with different number types. Euclid’s GCD algorithm is the same
whether you are using `Byte`

, `Short`

, `Int`

, `Long`

, or `BigInt`

. Why implement
it only for 1, when you could do it for all 5? Worse; why implement it 5
times, when you need only implement it once?

Another reason is that you want to push certain trade-offs, such as speed vs
precision to the user of your library, rather than making the decision for
them. `Double`

is fast, but has a fixed precision. `BigDecimal`

is slow, but
can have much higher precision. Which one do you use? When in doubt, let
someone else figure it out!

A last great reason is that it let’s you write less tests and can make testing much less hairy.

So, what does a generic version of Euclid’s GCD algorithm look like? Spire strives to make generic numeric code look more or less like what you’d write for a direct implementation. So, let’s let you compare; first up, the direct implementation:

```
def euclidGcd(x: Int, y: Int): Int =
if (y == 0) x
else euclidGcd(y, x % y)
```

With Spire, we can use the `spire.math.Integral`

type class to rewrite this as:

```
import spire.math.Integral
import spire.implicits._
def euclidGcd[A: Integral](x: A, y: A): A =
if (y == 0) x
else euclidGcd(y, x % y)
```

The 2 methods are almost identical, save the `Integral`

context bound.
`Integral`

gives us many methods we expect integers to have, like addition,
multiplication, and euclidean division (quotient + remainder).

Because Spire provides default implicit instances of `Integral`

for all of the
integral types that come in the Scala standard library, we can immediately use
`euclidGcd`

to find the GCD of many integer types:

```
euclidGcd(42, 96)
euclidGcd(42L, 96L)
euclidGcd(BigInt(42), BigInt(96))
```

This is much better than writing 5 different versions of the same algorithm!
With Spire, you can actually do away with `euclidGcd`

altogether, as `gcd`

comes with `Integral`

anyways:

```
spire.math.gcd(BigInt(1), BigInt(2))
```

Another benefit of generic numeric programming, is that you can push the choice
of numeric type off to someone else. Rather than hardcode a method or data
structure using `Double`

, you can simple require some `Fractional`

type.

I actually first found a need for generic numeric programming after I had
implemented a swath of algorithms with double precision floating point
arithmetic, only to find out that the minor precision errors were causing
serious correctness issues. The obvious fix was to just to use an exact type,
like `spire.math.Rational`

, which would’ve worked for many of my purposes.
However, many of the algorithms actually worked fine with doubles or even
integers, where as others required exact n-roots (provided by a number type
like `spire.math.Real`

). Being more precise meant everything got slower, even
when it didn’t need to be. Being less precise meant some algorithms would
occasionally return wrong answers. Abstracting out the actual number type
used meant I didn’t have to worry about these issues. I could make the choice
later, when I knew a bit more about my data, performance and precision
requirements.

We can illustrate this using a simple algorithm to compute the mean of some numbers.

```
import spire.math._
import spire.implicits._
// Note: It is generally better to use an incremental mean.
def mean[A: Fractional](xs: A*): A = xs.reduceLeft(_ + _) / xs.size
```

Here, we don’t care what type `A`

is, as long as it can be summed and divided.
If we’re working with approximate measurements, perhaps finding the mean of a
list of `Double`

s is good enough:

```
mean(0.5, 1.5, 0.0)
// = 0.6666666666666666
```

Or perhaps we’d like an exact answer back instead:

```
import spire.math.Rational
mean(Rational(1, 2), Rational(3, 2), Rational(0))
// = Rational(2, 3)
```

The main thing here is that as a user of the `mean`

function, I get to choose
whether I’d prefer the speed of `Double`

or the precision or `Rational`

. The
algorithm itself looks no different, so why not give the user the choice?

One of the best things is that if you write test code that abstracts over the number type, then you can re-use your tests for many different types. Spire makes great use of this, to ensure instances of our type classes obey the rules of algebra and that the number types in Spire (Rational, Complex, UInt, etc) are fundamentally correct.

There is another benefit though – you can ignore the subtleties of floating
point arithmetic in your tests if you want! If your code works with any number
type, then you can test with an exact type such as `spire.math.Rational`

or
`spire.math.Real`

. No more epsilons and NaNs. You shouldn’t let this excuse you
from writing numerically stable code, but it may save you many false negatives
in your build system, while also making you more confident that the fundamentals
are correct.

This is a big topic, deserving of its own blog post (you know who you are), so I’ll leave this here.

We’ve already seen `Integral`

, which can be used wherever you need something
that acts like an integer. We also saw the modulus operator, `x % y`

, but not
integer division. Spire differentiates between *integer division* and *exact
division*. You perform integer division with `x /~ y`

. To see it in action,
let’s use an overly complicated function to negate an integer:

```
import spire.math._
import spire.implicits._
def negate[A: Integral](x: A) = -(42 * (x /~ 42) + x % 42)
```

Instances of `Integral`

exist for `Byte`

, `Short`

, `Int`

, `Long`

and `BigInt`

.

Another type class Spire provides is `Fractional[A]`

, which is used for things
that have “exact” division. “Exact” is in quotes, since `Double`

or
`BigDecimal`

division isn’t really exact, but it’s close enough that we give
them a pass. `Fractional`

also provides `x.sqrt`

and `x nroot k`

for taking the
roots of a number.

```
def distance[A: Fractional](x: A, y: A): A = (x ** 2 + y ** 2).sqrt
```

Note that `Fractional[A] <: Integral[A]`

, so anything you can do with
`Integral`

, you can do with `Fractional[A]`

too. Here, we can use `distance`

to calculate the length of the hypotenuse with `Double`

s, `Float`

s,
`BigDecimal`

s, or some of Spire’s number types like `Real`

or `Rational`

.

Lastly, you often have cases where you just don’t care if `/`

means exact or
integer division, or whether you are taking the square root of an `Int`

or a
`Double`

. For this kind of catch-all work Spire provides `Numeric[A]`

.

If you’ve already hit the types of problems solved by generic numeric
programming, then you may have seen that `scala.math`

also provides `Numeric`

,
`Integral`

, and `Fractional`

, so why use Spire? Well, we originally created
Spire largely due to the problems with the type classes as they exist in Scala.

To start, Scala’s versions aren’t specialized, so they only worked with boxed versions of primitive types. The operators in Scala also required boxing, which means you need to trade-off performance for readability. They also aren’t very useful for a lot of numeric programming; what about nroots, trig functions, unsigned types, etc?

Spire also provides many more useful (and specialized) type classes. Some are
ones you’d expect, like `Eq`

and `Order`

, while others define more basic
algebras than `Numeric`

and friends, like `Ring`

, `Semigroup`

, `VectorSpace`

,
etc.

There are many useful number types that are missing from Scala in Spire, such
as `Rational`

, `Complex`

, `UInt`

, etc.

Spire was written by people who actually use it. I somewhat feel like Scala’s
Numeric and friends weren’t really used much after they were created, other
than for Scala’s NumericRange support (ie. `1.2 to 2.4 by 0.1`

). They miss
many little creature comforts whose need becomes apparent after using Scala’s
type classes for a bit.

One of Spire’s goals is that the performance of generic code shouldn’t suffer. Ideally, the generic code should be as fast as the direct implementation. Using the GCD implementation above as an example, we can compare Spire vs. Scala vs. a direct implementation. I’ve put the benchmarking code up as a Gist.

```
gcdDirect: 29.981 1.00 x gcdDirect
gcdSpire: 30.094 1.00 x gcdDirect
gcdSpireNonSpec: 36.903 1.23 x gcdDirect
gcdScala: 38.989 1.30 x gcdDirect
```

For another example, we can look at the code to find the mean of an array.

```
meanDirect: 10.592 **1.00 x gcdDirect**
meanSpire: 10.638 **1.00 x gcdDirect**
meanSpireNonSpec: 13.434 **1.26 x gcdDirect**
meanScala: 19.388 **1.83 x gcdDirect**
```

Spire achieves these goals fairly simply. All our type classes are
`@specialized`

, so when using primitives types you can avoid boxing. We then use macros to remove
the boxing normally required for the operators by the implicit conversions.

Using `@specialized`

, both `gcdSpire`

and `meanSpire`

aren’t noticably slower
than the direct implementation. We can see the slow down caused by dropping
`@specialized`

in `gcdSpireNonSpec`

and `meanSpireNonSpec`

. The difference
between `gcdSpireNonSpec`

and `gcdScala`

is because Spire doesn’t allocate an
object for the `%`

operator (using macros to remove the allocation). The
difference is even more pronounced between `meanSpireNonSpec`

and `meanScala`

.

`Numeric`

, `Integral`

, and `Fractional`

The 3 type classes highlighted in this post are just the tip of the iceberg.
Spire provides a whole slew of type classes in `spire.algebra`

. This package
contains type classes representing a wide variety of algebraic structures,
such as `Monoid`

, `Group`

, `Ring`

, `Field`

, `VectorSpace`

, and more. The 3 type
classes discussed above provide a good starting point, but if you use Spire in
your project, you will probably find yourself using `spire.algebra`

more and
more often. If you’d like to learn more, you can watch my talk on abstract
algebra in Scala.

As an example of using the algebra package, `spire.math.Integral`

is simply
defined as:

```
import spire.algebra.{ EuclideanRing, IsReal }
trait Integral[A] extends EuclideanRing[A]
with IsReal[A] // Includes Order[A] with Signed[A].
with ConvertableFrom[A]
with ConvertableTo[A]
```

Whereas `spire.math.Fractional`

is just:

```
import spire.algebra.{ Field, NRoot }
trait Fractional[A] extends Integral[A] with Field[A] with NRoot[A]
```

Spire also adds many new useful number types. Here’s an incomplete list:

**spire.math.Rational**is a fast, exact number type for working with rational numbers,**spire.math.Complex[A]**is a parametric number type for complex numbers,**spire.math.Number**is a boxed number type that strives for flexibility of use,**spire.math.Interval**is a number type for interval arithmetic,**spire.math.Real**is a number type for exact geometric computation that provides exact n-roots, as well as exact division,**spire.math.{UByte,UShort,UInt,ULong}**unsigned integer types, and**spire.math.Natural**an arbitrary precision unsigned integer type.

Spire also provides better operator integration with `Int`

and `Double`

. For
instance, `2 * x`

or `x * 2`

will just work for any `x`

whose type has an
`Ring`

. On the other hand, Scala requires something like
`implicitly[Numeric[A]].fromInt(2) * x`

which is much less readable. This
also goes for working with fractions; `x * 0.5`

will just work, if `x`

has a
`Field`

.

Spire has a basic algebra that let’s us work generically with numeric types. It does this without sacrificing readability or performance. It also provides many more useful abstractions and concrete number types. This means you can write less code, write less tests, and worry less about concerns like performance vs. precision. If this appeals to you, then you should try it out!

There is some basic information on getting up-and-running with Spire in SBT on
Spire’s project page. If you have any further
questions, comments, suggestions, criticism or witticisms you can say what you
want to say on the Spire mailing list
or on IRC on Freenode in `#spire-math`

.