Friday, March 30, 2007

Evaluating a quadratic polynomial

Christopher Diggins writes:
I am about to release a version of Cat which supports optional named parameters. The turning point in making this decision to add support for named parameters to a version of Cat occured when I tried to implement a simple quadratic equation with names.

He shows the following program for evaluating a quadratic polynomial f(x)=ax+bx+c, which takes the values x, a, b, c as input:
define quadratic_point_free
[dupd [sqr] dip swap [*] dip] dipd [* +] dip +
This is not good stack code.

I'm not writing this to bash Christopher, but I think he has not had much experience actually programming code in a stack language. Implementing an interpreter for any language doesn't force one to practice the idioms and problem solving techniques employed when actually programming in the language; likewise, one doesn't need a deep understanding of the implementation of their language to use it. So I think Chris is committing the classical newbie mistake when learning a language -- trying to relate features of the language to what they already know, and are somewhat surprised to find that they don't map one to one. In Christopher's case, however, he has the freedom to change the language as he sees fit. Since I don't know Cat very well, I can't say if named parameters would be a good thing or not. However, I think a feature like this needs a more compelling use-case than encoding polynomial evaluations.

So why is this code in bad taste? First, there's no stack effect comment. It wasn't hard for me to see that the stack effect is ( x a b c -- f(x) ), but it should be documented. However, the real problem is that this code operates at a level of abstraction which is too low.

Since he is a compiler developer, I imagine Christopher wrote the above code by manually "compiling" it in his head. This is the wrong way to approach a stack language. Instead, I have found two alternative approaches to be far more fruitful:
  • Write your code to take advantage of stack language idioms, and combine existing words express your problem in a simple way. That is, implement a Domain Specific Language (DSL)!
  • Write code which writes tedious code and performs this compilation process for you. That is, implement a DSL with more hair!

I'll show how to apply both approaches to solving this problem in Factor. In both cases we use Factor features not found in Cat, but they would be trivial to add. Cat, Joy and Factor are all very closely related, and share the same evaluation model and representation of code as quotations. Factor is compiled and it's compiler transforms are unique among concatenative languages, but a hypothetical compiled Joy or Cat could implement the same functionality with very little effort. Almost everything I write below also applies to low-level Forth.

Without further ado, I'll start with the first approach.

Just because the formula at hand has three additions and two multiplications, doesn't mean the symbols in the program must map one-to-one. The set of polynomials with real (or complex) coefficients forms a vector space. Polynomials are vectors. Polynomial evaluation is the dot product:
ax^2+bx+c = (a,b,c).(x^2,x,1)
Integers are sequences in Factor, and generating a sequence of the form { 1 x x^2 ... } is easy by mapping over an integer with the ^ word:
: powers ( x n -- seq ) [ ^ ] map-with ;
As long as the coefficients of the quadratic polynomial are given as a single object on the stack, a sequence { c b a }, we are done:
: quadratic ( c,b,a x -- f(x) ) 3 powers v. ;
However, we may as well generalize to evaluating a polynomial of any degree;
: poly-eval ( x coeff -- f(x) ) [ length powers ] keep v. ;
For example, to compute f(x) = 1024x^3 + 115x^2 - 71x + 576, you write
{ 576 -71 115 1024 } poly-eval
Certainly a lot simpler than the equivalent using the + and * words directly. And of course it will work with bignums, rationals, floats, and complex numbers. Also the polynomial is now a first-class entity in the language; even though it is just a sequence of numbers, it still stands that we have identified the essence of a polynomial. We are not threading the polynomial through the expression as a series of scalars. This means we can now pass the polynomial around as a single entity, and manipulate it with sequence words; we can add two polynomials together, multiply a polynomial by a scalar, juggle the coefficients around.

So to sum things up, I would consider the following a better implementation of Christopher's function:
: powers ( x n -- seq ) [ ^ ] map-with ;
: poly-eval ( x coeff -- f(x) ) [ length powers ] keep v. ;

This is useful as a general library routine; and it is very simple considering what it does. Also note that the way we approached the problem was different. We solved a general problem and the special case fell out for free. (In fact, the libs/math library already implements polynomial evaluation. This whole discussion is moot.)

Now, let us consider the second approach to problem solving: metaprogramming.

Back to Chris's function; here is a literal translation into Factor:
: quadratic-eval ( x a b c -- f(x) )
>r >r dupd >r sq r> swap >r * r> r> r> >r * + r> + ;

This looks like the output of an infix DSL compiler. It would be easy to implement an infix math DSL in Factor (easier than embedding Factor features into applicative languages!) and it would probably generate something almost identical to the above output, given input like
INFIX" a*x*x+b*x+c"
However I won't show you how to build an infix DSL. Instead, I'll solve a more limited problem; polynomial evaluation using Horner's Rule. The result will be a poly-eval word called in the exact same manner:
{ 576 -71 115 1024 } poly-eval
But instead, a compile-time transform will generate efficient code which does not actually create or iterate over arrays, or call the ^ word at run time. Such an approach is overkill in most situations. However, if your program needs to evaluate a lot of polynomials, it can speed things up. It also shows that the Factor compiler is very easy to extend, because it is written in, and compiles, a stack language.

Since we will be defining a compile-time transform, we need to write a word which takes a sequence of coefficients, and generates a quotation -- the quotation receives the input on the stack and evaluates the polynomial at that input. This word will be recursive. Consider the base case first. We are evaluating a constant polynomial;
: poly-eval-base ( c -- quot )
[ nip ] curry ;

For example, the code to evaluate the constant polynomial { 3.5 } is 3.5 nip; keep in mind the top of the stack is the value of x (which is disregarded by constant functions).

Now if we have a way of computing some polynomial g(x), we write a word which computes x*g(x) + c:
: poly-eval-step ( quot c -- newquot )
[ >r dupd call * r> + ] curry curry ;

The recursion puts the two together:
: poly-eval-quot ( seq -- quot )
dup length 1 = [
first poly-eval-base
] [
unclip >r poly-eval-quot r> poly-eval-step
] if ;

The above three words primarily take sequences apart and create new ones. The recursion is really a left fold; a left fold combinator could be added to the Factor library, but I haven't needed it yet.

We can test this code first:
  3 { 1 2 2 } poly-eval-quot call .

In fact, it is true that 25 = 1 + 2*3 + 2*3*3. As it stands, this new code generating implementation is not any more efficient. Instead of creating an intermediate array of exponents, we create an intermediate quotation of code. The trick is that if the sequence of coefficients is known at compile time, we can generate the quotation at compile time too:
: poly-eval poly-eval-quot call ;

\ poly-eval 1 [ poly-eval-quot ] define-transform
Consider the following word:
: foo { 11.4 23.7 78.22 99.6 -174.67 } poly-eval ;

We can define a specializer for it. This tells the compiler that while the word can take any type of object as input, a certain type is anticipated to happen more than most others; in our case, we optimize for the case where the input is a floating point value:
\ foo { float } "specializer" set-word-prop
Now, take a look at the PowerPC disassembly of the case where the input is indeed a floating point value:
0x059ea500:     lwz     r11,0(r14)
0x059ea504: lfd f0,3(r11)
0x059ea508: lis r11,1438
0x059ea50c: ori r11,r11,42632
0x059ea510: lwz r11,0(r11)
0x059ea514: lfd f1,3(r11)
0x059ea518: fmul f0,f0,f1
0x059ea51c: lis r11,1438
0x059ea520: ori r11,r11,42648
0x059ea524: lwz r11,0(r11)
0x059ea528: lfd f1,3(r11)
0x059ea52c: fadd f0,f0,f1
0x059ea530: lwz r11,0(r14)
0x059ea534: lfd f1,3(r11)
0x059ea538: fmul f1,f1,f0
0x059ea53c: lis r11,1438
0x059ea540: ori r11,r11,42644
0x059ea544: lwz r11,0(r11)
0x059ea548: lfd f0,3(r11)
0x059ea54c: fadd f1,f1,f0
0x059ea550: lwz r11,0(r14)
0x059ea554: lfd f0,3(r11)
0x059ea558: fmul f0,f0,f1
0x059ea55c: lis r11,1438
0x059ea560: ori r11,r11,42640
0x059ea564: lwz r11,0(r11)
0x059ea568: lfd f1,3(r11)
0x059ea56c: fadd f0,f0,f1
0x059ea570: lwz r11,0(r14)
0x059ea574: lfd f1,3(r11)
0x059ea578: fmul f1,f1,f0
0x059ea57c: lis r11,1438
0x059ea580: ori r11,r11,42636
0x059ea584: lwz r11,0(r11)
0x059ea588: lfd f0,3(r11)
0x059ea58c: fadd f1,f1,f0
0x059ea590: lis r12,1
0x059ea594: ori r12,r12,39104
0x059ea598: lwz r12,0(r12)
0x059ea59c: lwz r11,4(r12)
0x059ea5a0: addi r11,r11,16
0x059ea5a4: stw r11,4(r12)
0x059ea5a8: addi r11,r11,-16
0x059ea5ac: li r12,43
0x059ea5b0: stw r12,0(r11)
0x059ea5b4: stfd f1,8(r11)
0x059ea5b8: ori r12,r11,5
0x059ea5bc: stw r12,0(r14)
0x059ea5c0: lis r12,1
0x059ea5c4: ori r12,r12,4492
0x059ea5c8: mtlr r12
0x059ea5cc: blrl
The generated assembly will look like gibberish to people who are not familiar with the internals of the Factor compiler (anyone except for me?), so you'll have to take my word for this, but in fact we have the following:
  • The compiler doesn't do a good job of eliminating redundant loads from memory, and the above code could be a lot better.
  • But...
  • There is no recursion, sequence iteration, or sequence allocation in the above assembly.
  • The only place where a floating point value is boxed on the heap is at the end, after the final value has been computed. The allocation is done with customized inline assembly.
  • There are no subroutine calls until the end, where the GC check is called.
  • The input value is unboxed (too many times) and the constant coefficient values are unnecessarily unboxed, because the compiled code refers to boxed float objects not unboxed float values. This could be improved. However unboxing a float is inlined, and only takes two instructions.
  • All this is despite the fact that the code was generated from the following high-level description of the problem:
    { 11.4 23.7 78.22 99.6 -174.67 } poly-eval
We can compare the performance of the two implementations using the following benchmark:
:my-poly { 11.4 23.7 78.22 99.6 -174.67 } poly-eval ;
\ my-poly { float } "specializer" set-word-prop
: benchmark 1000000 [ >float my-poly drop ] each ;

It evaluates a fixed polynomial on every integer between 0 and 1000000, discarding the results. All arithmetic is done in floating point. Here are the results:
ImplementationRun time (ms)
Compile-time transform102
This is a 27x speed improvement!

Here is the grand finale. Suppose the polynomial coefficients are not fixed at compile time, but compile from some input file which is loaded by the user while the program runs. The user might not be aware there is such a thing as a compiler. Well, in Factor, compile-time is any-time. At the time the file is loaded, you can generate the evaluation quotations and compile them just by calling compile-quot. This technique is widespread among high-level languages. There are documented cases where Common Lisp code can beat optimized C by using information only available at run-time to optimize code. In stack languages, higher-order functions and code generation come even more naturally than in Lisp.

Nothing we did called for named parameters. Instead, one has to learn to "context switch" into the stack language mindset. It is a very different way of looking at the world.


d166e8 said...

What you describe as a stack-based idiom, is in fact an array processing idiom. Not that there is anything wrong with that, but it is unfair to accuse me of not understanding the idioms of stack-based language usage by comparing apples and oranges. Of course I am aware of array-processing techniques, and while they can be effective I chose not to use it in this case. Particularly because I wanted the code to be representative of stack-based code.

Leaving out the comment was also intentional, because I was intending to show that stack shuffling operations obscure the intention of code.

Even though I was vaguely irritated by the criticisms, which I view as unfounded, I am still very appreciative of the attention you have given Cat, and I am impressed with your work on Factor. I wish you lots of continued success.

Slava Pestov said...

Christopher, I certainly didn't mean to irritate anybody with my post.

When I first started working on Factor, I was not very good at dealing with the stack either. I would drop down into Java (the implementation language at the time) to implement various things, even though they could be done in Factor, because I was lazy. It wasn't until I started working on the native Factor implementation that I really learned Factor, because then I actually had no choice but to write a lot of Factor code!

I guess my whole point was that your function is not representative of stack code, at least not stack code written by humans (as opposed to compilers, rewrite systems, etc).

Array processing and similar tools greatly minimize the need for local names; the fact that they can be expressed in a stack language without a lot of effort is a plus, in my book. As Cat matures, you might find that having a rich library offsets the need for complex language features. Good luck.

Unknown said...

Hi Slava!
I've been interested in Stack-based languages for a while now, and Factor sticks out as a extremely good implementation and language.

But the lack of variables still bother me. I've heard advice similar to yours on this quadratic evaluator before. (namely, to re-interpret the mathematical function). But that has always seemed to me to be avoiding the issue.

Consider the case where we have a derived formula already. It might involve many different parameters and functions... but these fell naturally out of the mathematical derivation.

Given that the function is already well defined, it seems to me that the process of coding it should be almost a systematic process.

The quadratic function was easy to re-factor. But I'm sure it's very possible to devise a formula that is very difficult to re-factor.

In these cases, I think named parameters are a powerful abstraction that shouldn't be overlooked.

Anyway, good job on Factor! I am continually amazed how clean and well designed everything is. From the language down to the libraries and even the editor!


PS: I'm interested right now in a language that avoids heap allocation but still has a good abstraction for dealing with lists (or generators). Do you know of any?

Slava Pestov said...

CuppoJava: Factor supports named variables also.