Evaluation of a power series

I feel bad for people who have to use programming languages without higher-order functions. I think this is most programmers.

For a statistics project I'm working on right now, I implemented a few statistical functions in Factor. Most non-trivial functions are evaluated in a power series, I found from some quick Wikipedia research. That's an infinite summation that can't be reduced to anything finite, but it (hopefully) converges. So the evaluation strategy is to keep calculating the sum until the change in the sum goes below a particular value--say 1e-8. And, just to make sure, we'll stop it after a number of steps in case it diverges--say 2000.

For me, the immediate instinct was to write a higher-order function, which I called sigma, that would take a quotation and feed it values 0 to 2000 (or whatever the maximum step was set to) and sum up the results until they got too small. So I came up with this:

SYMBOL: max-steps
SYMBOL: delta
2000 max-steps set-global
1.0e-8 delta set-global

: sigma ( quot -- sum )
    0 swap max-steps get [ swap >r ] rot
    [ r> swap [ + ] keep abs delta get > ] 3append all?
    [ "Series doesn't converge quickly enough" throw ] when ; inline

: sigma-with ( obj quot -- sum )
    >r [ swap ] curry r> append sigma ; inline

[I first wrote this combinator in the old, pass-everything-around style, but decided to go for a different, less-efficient-in-the-current-implementation style where I mutate quotations. Slava wrote an email about this, which I never got around to responding to. Eventually, this will be fast in Factor.]

Building on this, an (extremely naive and slow) implementation of sin on all complex numbers could be written this way:

: sin ( x -- sin(x) )
    [ 
        [ [ 2 * 1+ ^ ] keep -1 swap ^ * ] keep
        2 * 1+ factorial /
    ] sigma-with ;

As another example, here's a calculation of e. Here, factorial is calculated incrementally, so it's not ridiculously stupid like the calculation of sin above.

: e ( -- e )
    1 1.0 [ dup zero? [ drop ] [ * dup recip ] if ] sigma nip ;

Depending on how the delta is set, this approximation can be as accurate as floating point will allow. To test this, just do [ epsilon delta set e ] with-scope. Have you had enough Greek letters yet?

Anyway, the point is, in Java, this would be practically impossible. No, it wouldn't be impossible to calculate any of these things. But the abstraction of sigma could not exist. It would have to be repeated each time the power series behavior is needed. It surprises me to realize that, after so many years of programming language research, most people use programming languages that are incapable of expressing sigma in a way convenient enough to use.