Exercise for the Reader

November 25, 2008

F# Arithmetic Specialization Rocks!

Filed under: Functional programming, Numerical methods — Tags: , , — Seth Porter @ 2:17 am

With a title like that, you know you’re in for a rocking good time.

Anyway, every good side project needs at least one new language, and after a long time wandering in the imperative wasteland (C# and, more recently, Java), I decided to give F# a look. It’s amazingly production-ready and very nicely integrated with Visual Studio (syntax highlighting and IntelliSense in a type-theoretic language?!). I know, strong type models and a good interactive mode should leave you happy to work in Notepad, but I’ve gotten awfully used to support from a rich IDE. But what a thrill it was to see, in a nicely formatted tooltip no less, my function’s inferred type:

val myFun : 'a * 'a -> ('a -> 'b) -> ('b * 'b)

Anyway, fun is fun, but my theoretical reason for all this was to try out some numerical recipes without tying my hands too early (as far as working in single / double resolution, or even something more exotic like interval arithmetic). Without getting side-tracked with numerical methods, consider good old linear interpolation:

f(p0, p1, u) = p0 * (1-u) + p1 * u

where ‘u’ varies from 0 to 1, and ‘p0′ and ‘p1′ are the values you’re interpolating between. One thing to note is that ‘p0′ and ‘p1′ can, in general, be anything you want β€” floating point values, vectors, or anything that can be scaled (by ‘u’) and added together.

So, my first, simplest approximation was very straightforward:

let linearInterp0(p0, p1, u) = p0 * (1-u) + u * p1

Looks good, but the inferred type comes out as ‘val linearInterp0 : int * int * int -> int’. Hmm. Maybe I need to explicitly tell it to be more generic. But ‘let linearInterp0(p0 : ‘a, p1 : ‘a, u) = p0 * u + (1-u) * p1′ (where ‘a is a generic type; loosely equivalent to ‘T myFun<T>(T p0, T p1…)’ in Java) fails to compile, with the helpful warning that ’1-u’ “causes code to be less generic than indicated by the type annotations”.

Fair enough, looks like that “1″ is forcing everything to be an integer. (I discovered later that int is also the default type for under-specified math operations, but let’s keep things simple.) Changing it to ’1.0′ gives me a type of ‘float * float * float -> float’, or ’1.0f’ for ‘float32 * float32 * float32 -> float32′ (apparently F# uses ‘float’ for what I think of as ‘double’, and ‘float32′ for single-precision floats). But I want to operate on any of the above (and don’t forget about vectors!). There was also no obvious way to manage one type for P0 and P1, and also retain flexibility on whether to use singles or doubles for ‘u’, but one step at a time.

So, some Googling lead to this thread, which almost broke my spirit right there. (Note that the last reply, from “Andyman”, would have helped a lot, but wasn’t posted at the time.) Either you create your own “numeric” interface, and implement it for all the types you care about, or you use magic global associations functions and create pseudo-generic methods which fail at runtime if called with an inappropriate type. Anyway, I tried a trivial “sum” method using the GetNumericAssociation approach (if you’re following along at home, you’ll need to include a reference to FSharp.PowerPack.dll, and open the ‘Microsoft.FSharp.Math’ namespace):

let sumNA(p0 : 'a, p1 : 'a) =
    let numeric = GlobalAssociations.GetNumericAssociation<'a>()
    numeric.Add(p0, p1)

It has a wonderfully generic type, ‘a * ‘a -> ‘a (again, in Java terms something like ‘T myFunc<T>(T p0, T p1)’). In fact, it’s too generic; it will fail at runtime if GetNumericAssociation fails for the provided type. However, it works for float32 and float, so it’s worth seeing what we get out of the compiler. Using ILDASM on the compiled assembly, the truth comes out:

  IL_0000:  call       class [FSharp.Core]Microsoft.FSharp.Core.Option`1<
         class [FSharp.PowerPack]Microsoft.FSharp.Math.INumeric`1<!!0>>
         [FSharp.PowerPack]Microsoft.FSharp.Math.GlobalAssociations::TryGetNumericAssociation<!!0>()
  IL_0005:  call       instance !0 class [FSharp.Core]Microsoft.FSharp.Core.Option`1<
         class [FSharp.PowerPack]Microsoft.FSharp.Math.INumeric`1<!!A>>::get_Value()
  IL_000a:  ldarg.0
  IL_000b:  ldarg.1
  IL_000c:  tail.
  IL_000e:  callvirt   instance !0 class [FSharp.PowerPack]
         Microsoft.FSharp.Math.INumeric`1<!!A>::Add(!0,!0)

Without worrying about the details of the .Net assembly generics syntax, the key is in that last line β€” ‘callvirt’. We’re calling a virtual method defined in the ‘INumeric’ interface, albeit with some syntactic sugar. I could do that in C# (or even Java), with less confusion. We’ve also lost the nice equivalence of the mathematical definition and the source code. Sigh. I want something more like C++ templating, but with actual types (instead of almost-textual expansion and all the other downsides of C++).

But before giving up, I googled some more, and ran into the ‘inline’ keyword. After some ILDASM-driven tweaking (in particular, if you want to curry one argument, you should curry them all, or else you’ll end up creating and dissecting tuples at runtime), I came up with this version:

let inline linearInterp one p0 p1 u = p0 * (one-u) + u * p1

This doesn’t look like much of a win; it’s much like the ‘numeric association’ version above, except now we have to pass an explicit value for “1″. At least we’ve gotten back to being able to use infix operators instead of function call syntax, which is something. However, look at the inferred type:

val inline linearInterp :
   ^a ->  ^b ->  ^c ->  ^d ->  ^e
    when ( ^a or  ^d) : (static member ( - ) :  ^a *  ^d ->  ^h) and
         ( ^b or  ^h) : (static member ( * ) :  ^b *  ^h ->  ^j) and
         ( ^j or  ^k) : (static member ( + ) :  ^j *  ^k ->  ^e) and
         ( ^d or  ^c) : (static member ( * ) :  ^d *  ^c ->  ^k)

Whew. (If you’re following along at home, I annotated the arguments with explicit type variables, so that they’re in alphabetical order.) In English, this means (roughly) that:

  • ‘linearInterp’ is a function of four values, each with a generic type, returning a value of a fifth generic type.
  • The first argument (‘one’) and the fourth argument (‘u’) must provide a subtraction operator returning some type ‘^h’
  • The second argument (‘p0′) and that ‘^h’ type must provide a multiplication operator, returning some type ‘^j’
  • (Skipping a line) the third argument (‘p1′) and the fourth argument (‘u) must provide a multiplication operator, returning some type ^k
  • The types ^j and ^k (returned from the scaling operations on p0 and p1) must provide an addition operator returning some type ^e, which is the return type of the whole method

That’s precisely as generic as this method can be. It doesn’t have the false genericism of the GetNumericAssociation variant; instead, the type spells out exactly what’s required for the function to make sense, and no more. By providing different values for ‘one’, such as 1.0 (double-precision), or 1.0f (single-precision), or even a hypothetical IntervalArithmetic.One, we can produce a family of specialized versions. For instance,

let linearFloat p0 p1 u = linearInterp 1.0 p0 p1 u
let linearSingle p0 p1 u = linearInterp 1.0f p0 p1 u
let bar = linearFloat 3.0 10.0 0.75

(and yes, bar = 8.25). (For the functional fans in the audience β€” yes, you could specify linearFloat just as well as ‘let linearFloat = linearInterp 1.0′, but the code generation is different in this case; it produces a closure rather than a static function.) But the proof is in the pudding; it’s quite possible that the compiler is (quite cleverly) inferring an unnamed interface much like the “INumeric” interface above, with Multiply and Add methods, and just providing syntactic sugar to make the calls look like infix operators. Hell, that’s what I’d do if I had to write a compiler for this sort of spontaneous hyper-abstraction. So, fire up ILDASM again and let’s see what the new methods look like:

.method public static float64  linearFloat(float64 p0,
                                           float64 p1,
                                           float64 u) cil managed
{
  // Code size       18 (0x12)
  .maxstack  5
  IL_0000:  ldarg.0
  IL_0001:  ldc.r8     1.
  IL_000a:  ldarg.2
  IL_000b:  sub
  IL_000c:  mul
  IL_000d:  ldarg.2
  IL_000e:  ldarg.1
  IL_000f:  mul
  IL_0010:  add
  IL_0011:  ret
}

Again, don’t worry too much about the details of the assembly (hint: it’s a stack-based computation model), but note that there isn’t a function call in sight. Those are real floating point (double-precision) operations, right there in the generated code. Also note that the “one” constant has been inlined as a real constant (loaded on line IL_0001), not a closure parameter or anything funky. In fact, I wrote an equivalent method in C#, explicitly on doubles, and the assembly is identical! (I’ve never seen that happen between two high level languages, though to be fair the common calling conventions etc. of the CLR make it easier.)

So we have all the benefits of writing functional, lexically nice code, with extreme levels of genericism, but at the end of the day we’re emitting the same bytecode that we’d get with hand-written C#. The very same method can be specialized for single-precision floats (perhaps for real-time rendering tasks), or used with any user-defined types which provide the appropriate operators.

To live up to the title, an interesting exercise will be left for the reader. Consider this minimal discriminated type, which might be the (much-abbreviated) basis of a symbolic math library:

type 'a expression =
  | Constant of 'a
  | Variable of string
  | Negate of 'a expression
  | Times of 'a expression * 'a expression
  | Plus of 'a expression * 'a expression
  static member (+) (e1 : 'a expression, e2 : 'a expression) = Plus(e1, e2)
  static member (*) (e1 : 'a expression, e2: 'a expression) = Times(e1, e2)
  static member (-) (e1: 'a expression, e2 : 'a expression) = Plus(e1, Negate(e2))

Given this type definition, which includes the *, + and – operators, we can define

let linearSymbolic p0 p1 u = linearInterp exprOne p0 p1 u

and exercise it a little (the > is the prompt from an interactive F# session, the double-semicolon triggers evaluation):

> linearSymbolic (Variable("p0")) (Variable("p1")) (Variable("u"));;
val it : float expression
= Plus
    (Times (Variable "p0",Plus (Constant 1.0,Negate (Variable "u"))),
     Times (Variable "u",Variable "p1"))

This may whet your appetite if you’ve ever had to do this sort of thing in C++; at this point I’d still be thinking about how I was going to write the parser.

So, the exercise for the reader is… well basically it’s to consider just how freaking cool it is that the exact same function definition can be used in full-on bulk computation code or in classic functional-language-style symbolic work. But while you’re at it, try to think of an equivalently cool trick to get back from the symbolic rep to the evaluatable version (for example, after performing some nice simplifications in symbolic form). I’m not sure it can be done without getting into compiler parse-tree structures (which is appears that F# makes available, so I may pursue this). It’s easy to write the classic evaluate method, but performance won’t even resemble the fully-inlined version we had above.

Ah, but it wouldn’t be an exercise for the reader if it was easy.

About these ads

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

The Silver is the New Black Theme. Blog at WordPress.com.

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: