Iterability of exp(x)-1

[andydude]

I am able to compute them so fast because of 4 secrets:

1. I have Mathematica.
   
2. I know that (IDM+Interpolation) is faster than (PDM+Powers)
   
3. I use the Double-Binomial instead of Lagrange Interpolation.
   
4. I use a custom IDM method.
   

With (PDM+Powers), you are effectively computing the interpolation of a triangular matrix, which is \( \frac{n(n+1)}{2} \) interpolations when we only need \( n \) interpolations. This indicated to me that the IDM method was a better way to go. After I realized that, it was only a matter of comparing the timings of Lagrange Interpolation vs. Double-Binomial expansion to see which was faster.

First, I will start with the straight-forward Mathematica code:

Code:

size = 10;

coeff = Join[{0}, Map[InterpolatingPolynomial[#, t] &,
  Transpose@Table[Series[Nest[(Exp[#] - 1) &, x, k],
  {x, 0, size}][[3]], {k, 1, size}]]];

TableForm@Expand@coeff

and now for the highly-optimized Mathematica code:

Code:

(* This is not built-in for some reason... *)
Unprotect[Binomial];
Binomial[t_, k_Integer?Positive] := Product[t - j, {j, 0, k - 1}]/k!;
Protect[Binomial];

(* This assumes x0 == 0 *)
ParabolicIDM[expr_, {x_, 0, size_}] := Module[{ret, ser},
  ret = {};
  ser = Series[x, {x, 0, size}];
  Do[
    ser = Series[expr /. x -> Normal[ser], {x, 0, size}];
    AppendTo[ret, Join[{0}, ser[[3]]]],
  {size + 1}];
  Join[{UnitVector[size + 1, 2]}, ret]
];

(* Then use the custom IDM as follows *)
size = 50;
matrix = ParabolicIDM[Exp[x] - 1, {x, 0, size}];
coeff = Table[Sum[(-1)^(n - k - 1) matrix[[k + 1, n + 1]]
  Binomial[t, k] Binomial[t - k - 1, n - k - 1],
  {k, 0, n - 1}], {n, 0, size}];

(* Then do what you like to 'coeff', for example: *)
TableForm@Expand@coeff[[1 ;; 20]]

Using this code, I find that I can get 100 coefficients in under 25 sec, and 150 coefficients in under 2 mins. As you can see, the custom IDM method is truncating the series at every step, so you're not saving the iteration of a 50-degree polynomial composed with a 50-degree polynomial which otherwise would be 2500-degree polynomial only after the second iterate. This "chopping" keeps the iterations within the size you are interested in, so you aren't calculating things that you don't care about.

Andrew Robbins