Coefficient function is slow

Question

Please consider:

Clear[x]
expr = Sum[x^i, {i, 15}]^30;

CoefficientList[expr, x]; // Timing
Coefficient[Expand@expr, x, 234]; // Timing
Coefficient[expr, x, 234]; // Timing

{0.047, Null}

{0.047, Null}

{4.93, Null}

Help states:

Coefficient works whether or not expr is explicitly given in expanded form.

Is there a reason why Coefficient needs to be so slow in the last case?

Answer 1

Here is a hack which may enable your code to be fast, but I don't guarantee it to always work correctly:

ClearAll[withFastCoefficient];
SetAttributes[withFastCoefficient, HoldFirst];
withFastCoefficient[code_] :=
   Block[{Binomial},
     Binomial[x_, y_] := 10 /; ! FreeQ[Stack[_][[-6]], Coefficient];
     code]

Using it, we get:

In[58]:= withFastCoefficient[Coefficient[expr,x,234]]//Timing
Out[58]= {0.172,3116518719381876183528738595379210}

The idea is that, Coefficient is using Binomial internally to estimate the number of terms, and then expands (calls Expand) if the number of terms is less than 1000, which you can check by using Trace[..., TraceInternal->True]. And when it does not expand, it computes lots of sums of large coefficient lists dominated by zeros, and this is apparently slower than expanding, for a range of expressions. What I do is to fool Binomial into returning a small number (10), but I also tried to make it such that it will only affect the Binomial called internally by Coefficient:

In[67]:= withFastCoefficient[f[Binomial[7,4]]Coefficient[expr,x,234]]
Out[67]= 3116518719381876183528738595379210 f[35] 

I can not however guarantee that there are no examples where Binomial somewhere else in the code will be computed incorrectly.

EDIT

Of course, a safer alternative that always exists is to redefine Coefficient using the Villegas - Gayley trick, expanding an expression inside it and calling it again:

Unprotect[Coefficient];
Module[{inCoefficient},
   Coefficient[expr_, args__] :=
      Block[{inCoefficient = True},
         Coefficient[Expand[expr], args]] /; ! TrueQ[inCoefficient]
];
Protect[Coefficient];

EDIT 2

My first suggestion had an advantage that we defined a macro which modified the properties of functions locally, but disadvantage that it was unsafe. My second suggestion is safer but modifies Coefficient globally, so it will always expand until we remove that definition. We can have the best of both worlds with the help of Internal`InheritedBlock, which creates a local copy of a given function. Here is the code:

ClearAll[withExpandingCoefficient];
SetAttributes[withExpandingCoefficient, HoldFirst];
withExpandingCoefficient[code_] :=
   Module[{inCoefficient},
      Internal`InheritedBlock[{Coefficient},
         Unprotect[Coefficient];
         Coefficient[expr_, args__] :=
            Block[{inCoefficient = True},
               Coefficient[Expand[expr], args]] /; ! TrueQ[inCoefficient];
         Protect[Coefficient];
         code
      ]
   ];

The usage is similar to the first case:

In[92]:= withExpandingCoefficient[Coefficient[expr,x,234]]//Timing
Out[92]= {0.156,3116518719381876183528738595379210}

The main Coefficient function remains unaffected however:

In[93]:= DownValues[Coefficient]
Out[93]= {}

Answer 2

Coefficient will not expand unless it deems it absolutely necessary to do so. This does indeed avoid memory explosions. I believe it has been this way since version 3 (I think I was working on it around 1995 or so).

It can also be faster to avoid expansion. Here is a simple example.

In[28]:= expr = Sum[x^i + y^j + z^k, {i, 15}, {j, 10}, {k, 20}]^20;

In[29]:= Coefficient[expr, x, 234]; // Timing

Out[29]= {0.81, Null}

But this next appears to hang in version 8, and takes at least a half minute in the development Mathematica (where Expand was changed).

Coefficient[Expand[expr], x, 234]; // Timing

Possibly some heuristics should be added to look for univariates that will not explode. Does not seem like a high priority item though.

Daniel Lichtblau

Answer 3

expr = Sum[x^i, {i, 15}]^30; 

scoeff[ex_, var_, n_] /; PolynomialQ[ex, var] := 
   ex + O[var]^(n + 1) /. 
    Verbatim[SeriesData][_, 0, c_List, nmin_, nmax_, 1] :> 
     If[nmax - nmin != Length[c], 0, c[[-1]]]; 

Timing[scoeff[expr, x, 234]]

seems quite fast, too.

Answer 4

After some experimentation following Rolf Mertig's answer, this appears to be the fastest method on the type of expression such as Sum[x^i, {i, 15}]^30:

SeriesCoefficient[expr, {x, 0, 234}]