Is there a fast product operation for PackedArrays?

Question

In Mathematica a vector (or rectangular array) containing all machine size integers or floats may be stored in a packed array. These objects take less memory, and some operations are much faster on them.

RandomReal produces a packed array when possible. A packed array can be unpacked with the Developer function FromPackedArray

Consider these timings

lst = RandomReal[1, 5000000];

Total[lst] // Timing
Plus @@ lst // Timing

lst = Developer`FromPackedArray[lst];

Total[lst] // Timing
Plus @@ lst // Timing

Out[1]= {0.016, 2.50056*10^6}

Out[2]= {0.859, 2.50056*10^6}

Out[3]= {0.625, 2.50056*10^6}

Out[4]= {0.64, 2.50056*10^6}

Therefore, in the case of a packed array, Total is many times faster than Plus @@ but about the same for a non-packed array. Note that Plus @@ is actually a little slower on the packed array.

Now consider

lst = RandomReal[100, 5000000];
Times @@ lst // Timing

lst = Developer`FromPackedArray[lst];
Times @@ lst // Timing

Out[1]= {0.875, 5.8324791357*10^7828854}

Out[1]= {0.625, 5.8324791357*10^7828854}

Finally, my question: is there a fast method in Mathematica for the list product of a packed array, analogous to Total?

I suspect that this may not be possible because of the way that numerical errors compound with multiplication. Also, the function will need to be able to return non-machine floats to be useful.

Answer 1

I've also wondered if there was a multiplicative equivalent to Total.

A really not so bad solution is

In[1]:= lst=RandomReal[2,5000000];
        Times@@lst//Timing
        Exp[Total[Log[lst]]]//Timing
Out[2]= {2.54,4.370467929041*10^-666614}
Out[3]= {0.47,4.370467940*10^-666614}

As long as the numbers are positive and aren't too big or small, then the rounding errors aren't too bad. A guess as to what might be happening during evaluation is that: (1) Provided the numbers are positive floats, the Log operation can be quickly applied to the packed array. (2) The numbers can then be quickly added using Total's packed array method. (3) Then it's only the final step where a non-machine sized float need arise.

See this SO answer for a solution that works for both positive and negative floats.

Let's quickly check that this solution works with floats that yield a non-machine sized answer. Compare with Andrew's (much faster) compiledListProduct:

In[10]:= compiledListProduct = 
          Compile[{{l, _Real, 1}}, 
           Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot], 
           CompilationTarget -> "C"]

In[11]:= lst=RandomReal[{0.05,.10},15000000];
         Times@@lst//Timing
         Exp[Total[Log[lst]]]//Timing
         compiledListProduct[lst]//Timing
Out[12]= {7.49,2.49105025389*10^-16998863}
Out[13]= {0.5,2.4910349*10^-16998863}
Out[14]= {0.07,0.}

If you choose larger (>1) reals, then compiledListProduct will yield the warning CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. and will take some time to give a result...

---

One curio is that both Sum and Product can take arbitrary lists. Sum works fine

In[4]:= lst=RandomReal[2,5000000];
         Sum[i,{i,lst}]//Timing
         Total[lst]//Timing
Out[5]= {0.58,5.00039*10^6}
Out[6]= {0.02,5.00039*10^6}

but for long PackedArrays, such as the test examples here, Product fails since the automatically compiled code (in version 8.0) does not catch underflows/overflows properly:

In[7]:= lst=RandomReal[2,5000000];
         Product[i,{i,lst}]//Timing
         Times@@lst//Timing
Out[8]= {0.,Compile`AutoVar12!}
Out[9]= {2.52,1.781498881673*10^-666005}

The work around supplied by the helpful WRI tech support is to turn off the product compilation using SetSystemOptions["CompileOptions" -> {"ProductCompileLength" -> Infinity}]. Another option is to use lst=Developer`FromPackedArray[lst].