Is there an operating system for PackedArrays?

Question

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

RandomReal creates a packed array RandomReal possible. A packed array can be unzipped using the Developer FromPackedArray function

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 an unpackaged array. Note that Plus @@ is actually a bit slower in a 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: is there a quick method in Mathematica for a packed array list product similar to Total ?

I suspect this may not be possible due to the fact that numerical errors are related to multiplication. In addition, the function will need to be able to return non-automatic floats to be useful.

Answer 1

I also wondered if there was a multiplicative equivalent of Total .

A really not that bad decision

 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} 

While the numbers are positive and not too large or small, rounding errors are not so bad. The hunch about what might happen in the evaluation is as follows: (1) If the numbers are positive floats, the Log operation can be quickly applied to a packed array. (2) Then the numbers can be quickly added using the Total method of packed arrays. (3) Then this is only the last step at which a non-plate float will appear.

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

Let me quickly check that this solution works with floats that give a non-standard answer. Compare with Andrew (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 select larger ones ( >1 ), then compiledListProduct will give a warning CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. and it will take some time to give a result ...

---

One curio is that both Sum and Product can accept arbitrary lists. Sum works great

 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 PackedArray s, such as the test examples here, Product fails because the automatically compiled code (in version 8.0) will not catch the overflow / overflow 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 provided by useful WRI technical support is to disable product compilation using SetSystemOptions["CompileOptions" -> {"ProductCompileLength" -> Infinity}] . Another option is to use lst=Developer`FromPackedArray[lst] .

Answer 2

First, to avoid confusion, take a look at an example whose results all appear as precision numbers of hardware machines that should be less

 In[1]:= $MaxMachineNumber Out[1]= 1.79769*10^308 

In your general example, this was already a nice (and quick) property. Here's an example from your Times example using machine numbers:

 In[2]:= lst = RandomReal[{0.99, 1.01}, 5000000]; Times @@ lst // Timing Out[3]= {1.435, 1.38851*10^-38} 

Now we can use Compile to make a compiled function for this operation to work effectively:

 In[4]:= listproduct = Compile[{{l, _Real, 1}}, Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot]] Out[4]= CompiledFunction[{l},Module[{tot=1.},Do[tot*=x,{x,l}];tot],-CompiledCode-] 

This is much faster:

 In[5]:= listproduct[lst] // Timing Out[5]= {0.141, 1.38851*10^-38} 

Assuming you have a C compiler and Mathematica 8, you can also automatically compile the entire path to the C code. A temporary DLL is created and bound to Mathematica at runtime.

 In[6]:= compiledlistproduct = Compile[{{l, _Real, 1}}, Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot], CompilationTarget -> "C"] Out[6]= CompiledFunction[{l},Module[{tot=1.},Do[tot*=x,{x,l}];tot],-CompiledCode-] 

This gives a performance that is not much different from what the Mathematica built-in function would have:

 In[7]:= compiledlistproduct[lst] // Timing Out[7]= {0.015, 1.38851*10^-38} 

Please note that if your product really goes beyond $ MaxMachineNumber (or $ MinMachineNumber ), then you better stick to Apply[Times, list] . The same comment applies to Total if your results can be so large:

 In[11]:= lst = RandomReal[10^305, 5000000]; Plus @@ lst // Timing Out[12]= {1.435, 2.499873364498981*10^311} In[13]:= lst = RandomReal[10^305, 5000000]; Total[lst] // Timing Out[14]= {1.576, 2.500061580905602*10^311} 

Answer 3

Simon's method is fast, but it does not work with negative values. Combining it with your answer to my other question , here is a quick solution that handles negatives. Thanks, Simon.

Functions

 f = (-1)^(-1 /. Rule @@@ Tally@Sign@# /. -1 -> 0) * Exp@Total@Log@Abs@# &; 

---

Testing

 lst = RandomReal[{-50, 50}, 5000000]; Times @@ lst // Timing f@lst // Timing lst = Developer`FromPackedArray[lst]; Times @@ lst // Timing f@lst // Timing 

---

 {0.844, -4.42943661963*10^6323240} {0.062, -4.4294366*10^6323240} {0.64, -4.42943661963*10^6323240} {0.203, -4.4294366*10^6323240}