F# Math (IV.) - Writing generic numeric code

Generic numeric code is some calculation that can be used for working with multiple different numeric types including types such as int, decimal and float or even our own numeric types (such as the type for clock arithmetic from the previous article of the series). Generic numeric code differs from ordinary generic F# code such as the 'a list type or List.map function, because numeric code uses numeric operators such as + or >= that are defined differently for each numeric type.

When writing simple generic code that has some type parameter 'T, we don’t know anything about the type parameter and there is no way to restrict it to a numeric type that provides all the operators that we may need to use in our code. This is a limitation of the .NET runtime and F# provides two ways for overcoming it.

• Static member constraints can be used to write generic code where the actual numeric operations are resolved at compile-time (and a generic function is specialized for all required numeric types). This approach makes resulting code very efficient and is generally very easy to use when writing a function such as List.sum.
• Global numeric associations (available in F# PowerPack) give us a way to obtain an interface implementing required numeric operations dynamically at runtime. This approach has some runtime overhead, but can be used for complex numeric types (such as Matrix<'T>).
• Combination of both techniques can be used to implement complex numeric type that is generic over the contained numeric values and has only a minimal runtime overhead.

Static member constraints are a unique feature of F# that is not available in other .NET languages, so if you're interested in writing numeric code for .NET, this may be a good reason for choosing F#. In C# or Visual Basic, you would be limited to the second option (which can be implemented in C#). In dynamic languages (like IronPython), everything is dynamic, so numeric computations can work with any numeric type, but will be significantly less efficient. In the rest of the article, we look at the three options summarized above.

This article is a part of a series that covers some F# and F# PowerPack features for numerical computing. Other articles in this series discuss matrices, defining custom numeric types and writing generic code. For links to other parts, see F# Math - Overview of F# PowerPack.

Using static member constraints

In this section, we’ll write a function with statically resolved type parameter. This type of parameter differs from usual generic type parameters in that it is resolved (and replaced with the actual type parameter) at compile-time.

When using statically resolved type parameter, we can use member constraints to specify that the actual type parameter must provide methods with certain types. These methods can then be called from the body of our generic function. Statically resolved type parameters can be used only when writing generic functions that are marked as inline. This means that a call to such function will be replaced with the body of the function during compilation, so this technique is useable only for relatively simple functions.

You can find more information about statically resolved type parameters [1] and inline functions [2] in the official MSDN documentation for F#, although I'll explain all the important bits in this article too.

Using built-in numeric operations

Even though it is possible to write a static member constraint manually, we can start more easily just by relying on the F# type inference. When writing an inline function that calls another inline function with member constraint, the constraint is automatically propagated. This means that all we need to do is to mark our function as inline and use only operations that are themselves written using static member constraints. This includes functions in the LanguagePrimitives module as well as standard F# operators such as * and +:

1: let inline halfSquare num =
2:   let res = LanguagePrimitives.DivideByInt num 2
3:   res * res
4: val inline halfSquare : ^a ->  ^b
5:   when ^a: (static member DivideByInt : ^a * int -> ^a)
6:    and ^a: (static member ( * ) : ^a * ^a -> ^b)

The function takes a numeric argument, divides it by two and then calculates the square of the result. The division cannot be written simply as num / 2, because this would use standard integer division and the compiler would infer num to be of type int. Instead, we use the DivideByInt function. This function is defined in the core libraries for all standard numeric types. For non-standard types, it uses a static member constraint and requires the type to have a DivideByInt member.

We can see the inferred static member constraints in the printed type signature of the function. The type parameter name is preceded with the caret (^) symbol denoting that it is a statically resolved type parameter. The when clause specifies the member constraints. As we can see, the type needs to provide DivideByInt method (which takes integer as the second parameter) and a multiplication operator.

The following two examples show that the function can be used with multiple different numeric types:

1: halfSquare 20.0f
2: val it : float32 = 100.0f
3: halfSquare 11M
4: val it : decimal = 30.25M

In the first example, we use the float32 type (which corresponds to System.Single or float in C#), while the second example uses decimal. It is worth noting that the DivideByInt member is not provided for integers, because it represents a precise floating point division. An alternative way to write the function above that would work for integers as well would use standard operators for division and addition (/ and +) and the GenericOne value from the LanguagePrimitives module.

Writing custom constraints

In the previous example, we relied on the F# type inference and let the compiler figure out the static constraints automatically, based on the operations that were used in our inline function. This worked well, because operators like / or numeric functions from the LanguagePrimitives module generate constraints.

However, the LanguagePrimitives module contains only a limited set of functions. If you want to write a function that works with any type with a specific member, you may need to write static member constraints explicitly. For example, both float32 (System.Single) and float (System.Double) have several members for dealing with IEEE floating point arithmetic.

The following function uses a static method IsInfinity to return an option value that is None when the floating-point number is infinity and Some otherwise:

1: let inline check< ^T when ^T : 
2:     (static member IsInfinity : ^T -> bool)> (num:^T) : option< ^T > =
3:   if (^T : (static member IsInfinity : ^T -> bool) (num)) then None
4:   else Some num

The declaration is quite verbose, but writing constraints explicitly is not usually needed very often - once you write a function that does a basic operation, it can be easily used thanks to type inference. The constraint is written as part of the generic parameter declaration. The syntax < ^T, when ... > specifies that the function has a (statically resolved) parameter ^T which must have a specified static member (IsInfinity) of a specified type.

To call the method inside the body of the function, we need to repeat the signature of the member constraint and we give it num as an argument. We can now call the function with different arguments:

1: check 42.0
2: val it : float option = Some(42)
3: check (1.0f / 0.0f)
4: val it : float32 option = null
5: check (1 / 2)
6: error FS0001: The type 'int' does not support 
7:               any operators named 'IsInfinity'

The first two calls are correct, because both float and float32 provide the required (static) member. In the third case, we get a compile-time error, because the F# compiler detects that the type int does not implement the required IsInfinity method. It is worth noting that the check function is not limited to the two types - if you define a custom numeric type that has the required static method, the check function will work as well. However, it is not possible to add IsInfinity method to existing types. Even if you implement it as an extension method, the F# compiler will ignore it (at least, in version 2.0).

As discussed in the introduction, static member constraints are useable only when writing an inline function. In the next section, we’ll explore the second technique of writing generic numerical code in F#.

Using global numeric associations

Using global numeric associations, we can obtain an interface providing implementations of basic operations for a specified numeric type at runtime. This is done using a global table of associations maintained by the F# PowerPack library [3]. The table contains numeric associations for all standard F# numeric types automatically and when you define a new type, you can register it with the table (see the previous article for an example). As a result any numeric code implemented using numeric associations will work with the newly defined type as well.

In this section, we’ll implement a type Quadruple<'T>. The type stores four (possibly) different values of type 'T and provides operations for point-wise addition and multiplication. You can think of it as a very simple vector type. There are two ways to implement the operations of the type:

• Operators of a type that are defined as static methods can be marked as inline and so they can use static member constraints. However, this approach doesn't work for instance members.
• When implementing a more complex type, we can pass around an interface INumeric<'T> that provides operations for the underlying numeric type.

The quadruple type implemented in this article is quite simple and we could use the first approach. However, I'd like to demonstrate the second alternative (which would be needed for more complex types like matrices). We'll get back to using static member constraints with quadruple in the last section.

As already mentioned, global associations are provided by the F# PowerPack library, so we first need to add a reference to FSharp.PowerPack.dll. The core part of the Quadruple<'T> type is straightforward. The constructor takes four parameters of type 'T and we expose them as members. More interestingly, the constructor also takes a value of type INumeric<'T> which provides numeric operations for the type 'T:

 1: type Quadruple<'T>(a:'T, b:'T, c:'T, d:'T, ops:INumeric<'T>) =
 2:   member x.Item1 = a
 3:   member x.Item2 = b
 4:   member x.Item3 = c
 5:   member x.Item4 = d 
 6: 
 7:   /// Expose implementation of operations for 'T
 8:   member x.Operations = ops
 9: 
10:   /// Create quadruple and retrieve numeric operations dynamically
11:   static member CreateRuntime(a, b, c, d) = 
12:     let ops = GlobalAssociations.GetNumericAssociation<'T>()
13:     new Quadruple<_>(a, b, c, d, ops)  
14: 
15:   /// Format quadruple as a string
16:   override x.ToString() = 
17:     sprintf "%A" (a, b, c, d)

We don’t expect the users of our type to obtain the numeric association themselves, so the type will be typically created using a method CreateRuntime. The method obtains numeric association dynamically using GetNumericAssociation. The function takes a single type parameter, which should be a numeric type and returns an implementation of the INumeric<'T> interface. In the above example, we don’t use the interface for anything yet. However, we expose it via a public property, so that we can use it for implementing overloaded operators.

The following listing shows code for * and + operators that we can add to our Quadruple<'T> type. Both of them are implemented as point-wise operations meaning that they multiply first element of the first quadruple with the first element of the second quadruple and so on:

 1: static member (+) (q1:Quadruple<_>, q2:Quadruple<_>) =
 2:   let inline ( +? ) a b = q1.Operations.Add(a, b)
 3:   Quadruple<_>(q1.Item1 +? q2.Item1, q1.Item2 +? q2.Item2, 
 4:                q1.Item3 +? q2.Item3, q1.Item4 +? q2.Item4,
 5:                q1.Operations)
 6: 
 7: static member (*) (q1:Quadruple<_>, q2:Quadruple<_>) =
 8:   let inline ( *? ) a b = q1.Operations.Multiply(a, b)
 9:   Quadruple<_>(q1.Item1 *? q2.Item1, q1.Item2 *? q2.Item2, 
10:                q1.Item3 *? q2.Item3, q1.Item4 *? q2.Item4,
11:                q1.Operations)

The operators have exactly the same structure. Both of them take two quadruples as parameters. Next, we define a simple local operator (just to make the syntax more convenient) that adds or multiplies two numbers (of type 'T) using the Add or Multiply method of the numeric association stored in the first quadruple. Since both of the quadruples have the same type, we could as well use the numeric association from the second one, because they are the same. Finally, we use the local helper operator to implement the point-wise operation and construct a new quadruple as the result. Note that we pass the numeric association as an argument to the newly constructed quadruple, so that it doesn’t need to access the associations table again.

The following listing demonstrates how quadruple works with integer and floating-point values:

 1: let quad a b c d = Quadruple<_>.CreateRuntime(a, b, c, d)
 2: val quad : 'a -> 'a -> 'a -> 'a -> Quadruple<'a>
 3: 
 4: (quad 1 2 3 4) + (quad 4 3 2 1)
 5: val it : Quadruple<int> = (5, 5, 5, 5)
 6: 
 7: let q1 = quad 1.0 2.0 3.0 4.0
 8: let q2 = quad 0.4 0.3 0.2 0.1
 9: (q1 + q2) * q1
10: val it : Quadruple<float> = (1.4, 4.6, 9.6, 16.4)

The snippet first defines a helper function quad that creates a quadruple. The type isn’t restricted by any member constraints. However, the implementation works only with numeric types, which means that when we call quad with non-numeric parameters, it will throw an exception.

It is also worth noting that the GetNumericAssociation function is called only when creating quadruples. Both of the operators extract the numeric association from one of the parameters, so once we construct quadruples, the rest of the computation can be efficient. This explains why this technique can be used for example by F# Matrix<'T> type, where we need to avoid unnecessary overheads.

Combining member constraints and INumeric

Using the GetNumericAssociations method to obtain an implementation of INumeric<'T> interface from a table maintained by F# PowerPack has some runtime overhead and it requires the numeric type to be registered in the associations table. However, it is possible to combine the approach based on static member constraints and the idea to capture the operations using the INumeric<'T> interface.

In this section, we'll extend the Quadruple<'T> type with a static member Create that constructs a value of the type more efficiently than CreateRuntime that we implemented previously. The CreateRuntime method may still be useful if you wanted to create quadruples dynamically using reflection. Before looking at a better way to implement the Create method, let's measure the performance of CreateRuntime using a simple F# Interactive command:

1: #time
2: 
3: // Add quadruples in a simple imperative loop
4: let mutable sum = quad 0 0 0 0
5: for i in 0 .. 10000000 do
6:   sum <- sum + (quad i (i/2) (i/3) (i/4))
7: Real: 00:00:02.536, CPU: 00:00:02.511, (...)
8: val mutable sum : Quadruple<int> = 
9:   (-2004260032, -1004630016, -2103075776, 1642668640)

The snippet creates 10 million quadruples and adds them in an imperative for loop. When executed as an unoptimized F# code in F# Interactive, the operation takes 2 and half seconds.

To add the Create method, we can just modify the existing type. The method is marked as inline and it has a statically resolved type parameter (hat-type) named ^TStatic. In the body, we use an object expression to implement the interface INumeric<^TStatic> and then we use it to construct a quadruple:

 1: /// Initialize a quadruple and capture operations for working with the
 2: /// type ^TStatic using static member constraints
 3: static member inline Create
 4:     (a : ^TStatic, b: ^TStatic, c : ^TStatic, d : ^TStatic) =
 5:   let ops = 
 6:     // Implement INumeric interface using generic operators
 7:     { new INumeric< ^TStatic > with
 8:         member x.Add(a, b) = a + b 
 9:         member x.Abs(a) = abs a
10:         member x.Compare(a, b) = compare a b
11:         member x.Equals(a, b) = a = b
12:         member x.Multiply(a, b) = a * b
13:         member x.Negate(a) = -a
14:         member x.Subtract(a, b) = a - b 
15:         member x.Sign(a) = sign a
16:         member x.Zero = LanguagePrimitives.GenericZero 
17:         member x.One = LanguagePrimitives.GenericOne
18:         member x.ToString(_, _, _) = failwith "not implemented"
19:         member x.Parse(_, _, _) = failwith "not implemented" }
20: 
21:   // Create a quadruple and pass it an INumeric implementation 
22:   new Quadruple<_>(a, b, c, d, ops)

When creating an inline static method, the inline keyword needs to be placed after the member keyword. We don't need to explicitly define the statically resolved type parameter - if we just use it in the type signature, the compiler will make the method generic for us.

The key idea of the trick that we used here is that the object expression that implements INumeric<^TStatic> is also going to be inlined when someone calls Create. This means that the implementation will be specialized for the specific numeric type and all operations that are used to implement it (+, *, functions like compare etc.) will be replaced with the optimized code for a specific numeric type. The type signature of the Create method (hover over its name to see a tool tip) shows all the requirements on the ^TStatic type parameter.

The only overhead of the Create method is that it creates an instance of a simple object that implements the INumeric<'T>interface. On the other hand, the CreateRuntime method performed lookup in a dictionary, based on some token representing a System.Type (that had to be retrieved based on a generic type parameter). Let's measure the performance of our previous example using a new quad function:

1: let inline quad a b c d = Quadruple<_>.Create(a, b, c, d)
2: 
3: // Imperative loop using statically resolved 'quad' function
4: let mutable sum = quad 0 0 0 0
5: for i in 0 .. 10000000 do
6:   sum <- sum + (quad i (i/2) (i/3) (i/4))
7: Real: 00:00:00.758, CPU: 00:00:00.748, (...)
8: val mutable sum : Quadruple<int> = 
9:   (-2004260032, -1004630016, -2103075776, 1642668640)

The snippet starts by declaring a new version of quad helper function. This time, it is implemented as inline function that calls the new inline method Create. When we run the for loop using our new implementation, the running time drops from 2.5 seconds to about 0,7 seconds, which means that the new version is between three- and four-times faster! Moreover, the inlined code only creates an instance of a specialized class, so the resulting binary isn't larger than the original version .

Summary

In .NET, there is no way to say that a generic type parameter should be a numeric type that supports basic numeric operators and functions. This makes it difficult to write numeric code that performs some calculation, but works with multiple different types. Ideally, we'd like to be able to write a calculation (or a generic type) once and use it with float32 (System.Single) and float (System.Double) and possibly also decimal and int if the calculation does not require floating-point arithmetic. Attempts to write generic numeric code in C# [5, 6, 7] have either runtime overhead or make simple code look very complicated.

In this article, we've seen how to solve the problem in F#. Thanks to static member constraints and the inline keyword, it is easy to write simple numeric functions that are generic over the numeric type used. The F# library itself uses this technique to implement functions like List.sum (to sum elements of a list). Such functions work with all standard types, but also with custom numeric types (as discussed in previous article).

For more complex generic types, we need to store the numeric operations in an object. In F#, this can be easily done using INumeric<'T> type that is available in F# PowerPack. If we need to obtain the implementation of the interface only rarely, we can use GlobalAssociations module that provides a runtime lookup. However, we've also seen that it is possible to make this more efficient by combining the INumeric<'T> interface and static member constraints. As far as I'm aware, this is probably the best solution available for .NET programmers in terms of succinctness and efficiency.

References & Links

• [1] Statically Resolved Type Parameters (F#) - MSDN Documentation
• [2] Inline Functions (F#) - MSDN Documentation
• [3] F# PowerPack (source code and binary packages) - CodePlex
• [4] Numeric Computing - Real World Functional Programming on MSDN
• [5] Using generics for calculations - Rüdiger Klaehn at CodeProject
• [6] Generic operators - Marc Gravell, Jon Skeet
• [7] Simulating INumeric with dynamic in C# 4.0 - Luca Bolognese