Bezier evaluation in O(n)

The really stupid way to evaluate bezier curves is with recursion, which is \$O(2^n)\$, which can be lowered to \$O(n^2)\$ with memoization. De Casteljau's algorithm improves on this, and is still \$O(n^2)\$ but faster.

I misinterpreted De Casteljau's and thought it was the Bernstein form:

$$\sum_{i=0}^{n-1}{\binom{n}{i}p_i (1-t)^{n-i-1} t^i}$$

Where there are \$n\$ control points including the start and end named \$p_0, p_1,p_2, \dots, p_{n-1}\$.

Calculating the combination is normally \$O(n)\$ which would make this algorithm as a whole \$O(n^2)\$ since there are \$n\$ terms to sum, however, there is a cheat:

$$\binom{n}{k}=\prod_{i=1}^{k}{\frac{n+1-i}{k}}$$

We can store the result after each iteration, resulting in the entire sequence calculated in \$O(n)\$ time, making the algorithm as a whole \$O(n)\$.

\$2\leq n\leq 5\$ is the general use case (2 and 4 probably actually), so it's hardcoded. For \$n\geq 6\$, which uses this algorithm, I've done a ridiculous number of optimizations.

The smallest double increment (one mantissa bit) is \$2.2E{-16}\$, and based on some quick benchmarks with random data it never was more inaccurate than \$1E{-14}\$, so the inaccuracy is acceptable.

For very large \$n\$ (hundreds, 1200 was where I first saw it) the combination function exceeded the max value for doubles, causing the algorithm to return NaN, whereas De Casteljau's, though much slower, returned what was more or less the correct value.

The final code:

```
public static double bezier2(double a,double b,double t){
//Total of 3 floating point operations
// 1 2 3
return a+t*(b-a);
}

public static double bezier(double t,double... ds){
return bezier(ds,t);
}

public static double bezier(double[] ds,double t){
int count = ds.length;
switch(count){
case 0:throw new IllegalArgumentException("Must have at least two items to interpolate between");
case 1:return ds[0];
case 2:return bezier2(ds[0],ds[1],t);
case 3:{

You need to look at your unit tests, because there's a bad bug in chooseLongRange. int product = 1 should be long product = 1.

if(count<29){//See note 1
        choose = new double[halfn1];
        int[] chooseInt = chooseIntRange(n1,halfn);
        for(int i=0;i<halfn1;i++){
            choose[i]=chooseInt[i];
        }
    }else if(count<60){//See note 2
        choose = new double[halfn1];
        long[] chooseLong = chooseLongRange(n1,halfn);
        for(int i=0;i<halfn1;i++){
            choose[i]=chooseLong[i];
        }
    }else{
        choose = chooseDoubleRange(n1,halfn);
    }

IMO this is excessively complicated. Unless you're reusing chooseIntRange and chooseLongRange in contexts which require int[] or long[], why not make them return double[]? And then you can push in the case split for the type to use, exposing a single method and leaving the responsibility in the method where it really belongs.

I'm slightly puzzled by "But as it turns out, the extra cost of that single if statement and a redirect call meant that it was only worth redirecting for 6≤n≤9". What "single if statement"? Do you mean the overhead of the switch table lookup?

I get the impression that you care more about performance than numerical analysis, but I think it's worth making this point anyway.

There's an important difference between

a + (b - a) * t

and

a * (1 - t) + b * t

If you substitute \$(a, b, t) = (b, a, 1-t)\$ then the second will give you the same result, but the first won't necessarily. In particular, when \$t \approx 1\$ and \$b \ll a\$ the result should be very close to \$b\$, but in the first expression the subexpression \$b - a\$ will lose the precision of \$b\$ and it can't be recovered in the final result.

StackExchange Code Review Q#155781, answer score: 2