PackedArrays有快速的产品运营吗？

Question

在Mathematica中，包含所有机器大小的整数或浮点数的向量(或矩形数组)可以存储在压缩数组中。这些对象占用较少的内存，并且对它们的一些操作要快得多。

RandomReal会在可能的情况下生成压缩数组。可以使用Developer函数FromPackedArray对压缩数组进行解包

请考虑以下时间安排

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}

因此，在压缩数组的情况下，Total的速度比Plus @@快很多倍，但对于非压缩数组则大致相同。请注意，在压缩数组上，Plus @@实际上会稍微慢一点。

现在考虑一下

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}

最后，我的问题是: Mathematica中有没有类似于Total的压缩数组列表乘积的快速方法

我怀疑这可能是不可能的，因为数字误差与乘法复合的方式。此外，该函数需要能够返回有用的非机器浮点数。

Answer 1

首先，为了避免混淆，请看一个示例，它的结果都可以表示为硬件机器精度数字，这些数字必须都小于

In[1]:= $MaxMachineNumber

Out[1]= 1.79769*10^308

你的Total示例已经有了这个很好的(并且快速的)属性。以下是使用机器编号的Times示例的变体：

In[2]:= lst = RandomReal[{0.99, 1.01}, 5000000];
Times @@ lst // Timing

Out[3]= {1.435, 1.38851*10^-38}

现在，我们可以使用Compile创建一个编译函数来高效地执行此操作：

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-]

它要快得多：

In[5]:= listproduct[lst] // Timing

Out[5]= {0.141, 1.38851*10^-38}

假设你有一个C编译器和Mathematica 8，你也可以自动编译成C代码。创建一个临时DLL，并在运行时链接回Mathematica。

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-]

这使得性能与内置的Mathematica函数没有太大区别：

In[7]:= compiledlistproduct[lst] // Timing

Out[7]= {0.015, 1.38851*10^-38}

请注意，如果您的产品真的超越了$MaxMachineNumber (或$MinMachineNumber)，那么您最好坚持使用Apply[Times, list]。同样的评论也适用于Total，如果你的结果能达到那么大的话：

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 2

Simon的方法很快，但在负值时失败。将它与his answer to my other question相结合，这里是一个处理负片的快速解决方案。谢谢，西蒙。

函数

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

测试

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}