【问题标题】:Generate a list in Mathematica with a conditional tested for each element在 Mathematica 中生成一个列表,对每个元素进行条件测试
【发布时间】:2011-06-16 06:28:15
【问题描述】:

假设我们要生成一个素数列表p,其中p + 2 也是素数。

一个快速的解决方案是生成前 n 个素数的完整列表,并使用 Select 函数返回满足条件的元素。

Select[Table[Prime[k], {k, n}], PrimeQ[# + 2] &]

但是,这是低效的,因为它会在返回过滤后的列表之前将一个大列表加载到内存中。带有 Sow/Reap(或 l = {}; AppendTo[l, k])的 For 循环解决了内存问题,但它远非优雅,而且在 Mathematica 脚本中多次实现很麻烦。

Reap[
  For[k = 1, k <= n, k++,
   p = Prime[k];
   If[PrimeQ[p + 2], Sow[p]]
  ]
 ][[-1, 1]]

理想的解决方案是内置函数,它允许类似的选项。

Table[Prime[k], {k, n}, AddIf -> PrimeQ[# + 2] &]

【问题讨论】:

  • 所以本质上你想要一个更复杂的list comprehension...
  • 顺便说一句:NextPrime 文档中有一个使用 Select 的简单孪生素数生成器。
  • 很好的发现,我没听懂。当然,找到孪生素数只是我能想到的一个简单的例子。

标签: wolfram-mathematica


【解决方案1】:

我将把它更多地解释为一个关于自动化和软件工程的问题,而不是关于手头的具体问题,并且已经发布了大量的解决方案。 ReapSow 是收集中间结果的好方法(可能是符号设置中的最佳方法)。让我们把它概括一下,以避免代码重复。

我们需要的是写一个高阶函数。我不会做任何全新的事情,只会简单地打包您的解决方案以使其更普遍适用:

Clear[tableGen];
tableGen[f_, iter : {i_Symbol, __}, addif : Except[_List] : (True &)] :=
 Module[{sowTag},   
  If[# === {}, #, First@#] &@
       Last@Reap[Do[If[addif[#], Sow[#,sowTag]] &[f[i]], iter],sowTag]];

使用Do 优于For 的优点是循环变量是动态本地化的(因此,在Do 范围之外不会对其进行全局修改),而且Do 的迭代器语法更接近与Table 相比(Do 也稍快一些)。

现在,这里是用法

In[56]:= tableGen[Prime, {i, 10}, PrimeQ[# + 2] &]

Out[56]= {3, 5, 11, 17, 29}

In[57]:= tableGen[Prime, {i, 3, 10}, PrimeQ[# + 1] &]

Out[57]= {}

In[58]:= tableGen[Prime, {i, 10}]

Out[58]= {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

编辑

这个版本更接近你提到的语法(它需要一个表达式而不是一个函数):

ClearAll[tableGenAlt];
SetAttributes[tableGenAlt, HoldAll];
tableGenAlt[expr_, iter_List, addif : Except[_List] : (True &)] :=
 Module[{sowTag}, 
  If[# === {}, #, First@#] &@
    Last@Reap[Do[If[addif[#], Sow[#,sowTag]] &[expr], iter],sowTag]];

它还有一个额外的优势,即您甚至可以全局定义迭代器符号,因为它们在传递时是未经计算的并且是动态本地化的。使用示例:

In[65]:= tableGenAlt[Prime[i], {i, 10}, PrimeQ[# + 2] &]

Out[65]= {3, 5, 11, 17, 29}

In[68]:= tableGenAlt[Prime[i], {i, 10}]

Out[68]= {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

请注意,由于现在语法不同,我们必须使用 Hold-attribute 来防止传递的表达式 expr 过早评估。

编辑 2

根据@Simon 的要求,这里是多维度的概括:

ClearAll[tableGenAltMD];
SetAttributes[tableGenAltMD, HoldAll];
tableGenAltMD[expr_, iter__List, addif : Except[_List] : (True &)] :=
Module[{indices, indexedRes, sowTag},
  SetDelayed @@  Prepend[Thread[Map[Take[#, 1] &, List @@ Hold @@@ Hold[iter]], 
      Hold], indices];
  indexedRes = 
    If[# === {}, #, First@#] &@
      Last@Reap[Do[If[addif[#], Sow[{#, indices},sowTag]] &[expr], iter],sowTag];
  Map[
    First, 
    SplitBy[indexedRes , 
      Table[With[{i = i}, Function[Slot[1][[2, i]]]], {i,Length[Hold[iter]] - 1}]], 
    {-3}]];

这要简单得多,因为我必须Sow 索引和添加的值,然后根据索引拆分生成的平面列表。下面是一个使用示例:

{i, j, k} = {1, 2, 3};
tableGenAltMD[i + j + k, {i, 1, 5}, {j, 1, 3}, {k, 1, 2}, # < 7 &]

{{{3, 4}, {4, 5}, {5, 6}}, {{4, 5}, {5, 6}, {6}}, {{5, 6}, {6}}, {{6}}}

我将值分配给i,j,k 迭代器变量以说明此函数确实本地化了迭代器变量并且对它们可能的全局值不敏感。为了检查结果,我们可以使用Table,然后删除不满足条件的元素:

In[126]:= 
DeleteCases[Table[i + j + k, {i, 1, 5}, {j, 1, 3}, {k, 1, 2}], 
    x_Integer /; x >= 7, Infinity] //. {} :> Sequence[]

Out[126]= {{{3, 4}, {4, 5}, {5, 6}}, {{4, 5}, {5, 6}, {6}}, {{5, 6}, {6}}, {{6}}}

请注意,我没有进行广泛的检查,因此当前版本可能包含错误并需要更多测试。

编辑 3 - 错误修复

注意重要的错误修复:在所有函数中,我现在使用带有自定义唯一标签的Sow,以及Reap。如果没有此更改,当它们评估的表达式也使用Sow 时,这些函数将无法正常工作。这是Reap-Sow 的一般情况,类似于异常情况 (Throw-Catch)。

编辑 4 - SyntaxInformation

由于这是一个潜在有用的函数,让它更像一个内置函数是很好的。首先我们通过

添加语法高亮和基本参数检查
SyntaxInformation[tableGenAltMD] = {"ArgumentsPattern" -> {_, {_, _, _., _.}.., _.},
                                    "LocalVariables" -> {"Table", {2, -2}}};

然后,添加使用消息允许菜单项“制作模板”(Shift+Ctrl+k) 工作:

tableGenAltMD::usage = "tableGenAltMD[expr,{i,imax},addif] will generate \
a list of values expr when i runs from 1 to imax, \
only including elements if addif[expr] returns true.
The default of addiff is True&."

可以在this gist 中找到更完整和格式化的使用消息。

【讨论】:

  • +1。这是一个非常方便的功能。是否有可能使它像Table 一样用于多个迭代器并输出多维列表?目前,用iter__ 简单地替换iter_ 是可行的,但由于显而易见的原因,会产生一个扁平列表。
  • @Simon 我添加了一个多维版本,请参阅我的编辑。然而,它要复杂得多/晦涩难懂,也许也不那么优雅。
  • @Leonid:哎哟!少了很多优雅。它看起来像是您在编写它时才理解的代码类型。我了解它的工作原理(感谢您的描述),但仍然......虽然它比您的其他解决方案复杂得多,但速度也不算太差。
  • @Leonid:顺便说一句,tableGen[i + j, {i, 2}, {j, 2}] 不起作用,因为它认为{j,2}addif 函数。要使其工作,您需要限制 addif 的头部或执行一些其他模式检查。
  • @Vortico:我怀疑你将它添加到你的 Mma 内核中(除非你在 WRI 工作?) - 但我同意,我可以看到自己经常使用它。再次感谢列昂尼德!
【解决方案2】:

我认为 Reap/Sow 方法在内存使用方面可能是最有效的。一些替代方案可能是:

DeleteCases[(With[{p=Prime[#]},If[PrimeQ[p+2],p,{}] ] ) & /@ Range[K]),_List]

或者(这可能需要某种 DeleteCases 来消除 Null 结果):

FoldList[[(With[{p=Prime[#2]},If[PrimeQ[p+2],p] ] )& ,1.,Range[2,K] ]

两者都在内存中保存了一个从 1 到 K 的大整数列表,但 Prime 的范围在 With[] 构造内。

【讨论】:

    【解决方案3】:

    是的,这是另一个答案。包含 Reap/Sow 方法和 FoldList 方法的另一种替代方法是使用 Scan。

    result = {1};
    Scan[With[{p=Prime[#]},If[PrimeQ[p+2],result={result,p}]]&,Range[2,K] ];
    Flatten[result]
    

    同样,这涉及一长串整数,但中间的 Prime 结果没有被存储,因为它们在 With 的本地范围内。因为 p 是 With 函数范围内的常数,所以可以使用 With 而不是 Module,并获得一点速度。

    【讨论】:

    • 您可以通过使用纯函数来避免使用With 以避免重复评估:If[PrimeQ[# + 2], result = {result, #}] &amp;[Prime[#]] &amp;。这并不总是可能的,但在这里是可能的,因为您的 If 的主体不包含槽变量(#1 等)。
    • 感谢@Leonid:这些是我无法在工作中尝试的东西,因为(a)我没有安装 Mma 并且(b)我打算做其他事情.不过,是否有特别的理由要避免使用 With?
    • 我只是认为基于纯函数的方法更简单一些。此外,由于With 是一个作用域构造,并且必须解决可能的名称冲突,匿名纯函数可能会稍微快一些。
    • @Leonid:很公平,我只是不相信自己能在盲目的情况下编写这样的嵌套纯函数。
    • 这对我来说可能也是一件正确的事情。从长远来看,我完全不确定我的解决方案中失去可读性是否合理。
    【解决方案4】:

    你或许可以试试这样的:

    Clear[f, primesList]
    f = With[{p = Prime[#]},Piecewise[{{p, PrimeQ[p + 2]}}, {}] ] &;
    primesList[k_] := Union@Flatten@(f /@ Range[k]);
    

    如果你想要素数p和素数p+2,那么解决方案是

    Clear[f, primesList]
    f = With[{p = Prime[#]},Piecewise[{{p, PrimeQ[p + 2]}}, {}] ] &;
    primesList[k_] := 
      Module[{primes = f /@ Range[k]}, 
       Union@Flatten@{primes, primes + 2}];
    

    【讨论】:

    • 您为每个案例评估 Prime[p] 两次。我建议将函数 f 更改为使用 With[{p=Prime[#]},Piecewise... etc
    • @Verbeia:感谢您指出这一点,我忘了这样做。
    【解决方案5】:

    好吧,有人必须在某个地方为整个表的大小分配内存,因为事先不知道最终的大小是多少。

    在函数式编程之前的美好时光:),这种事情是通过分配最大数组大小来解决的,然后使用单独的索引插入它,这样就不会产生任何漏洞。像这样

    x=Table[0,{100}];  (*allocate maximum possible*)
    j=0;
    Table[ If[PrimeQ[k+2], x[[++j]]=k],{k,100}];
    
    x[[1;;j]]  (*the result is here *)
    
    {1,3,5,9,11,15,17,21,27,29,35,39,41,45,51,57,59,65,69,71,77,81,87,95,99}
    

    【讨论】:

    • 我不同意你的第一句话。为一个完整的表分配内存是最后的手段,也是一种资源浪费。在 Mathematica 中这样做的唯一原因是当我们知道大小的上限并准备好用内存换取速度时(考虑到 mma 性能调整的性质)。一般来说,当事先不知道大小时,解决办法是分配巨大的内存,而是动态地构造一个列表。在像 C 这样的低级语言中,可能会为此使用链表。链表也可以在 Mathematica 中使用,但 Reap-Sow 方法更胜一筹。
    • 顺便说一句,函数式编程几乎和命令式一样古老。 LISP 只比 Fortran 小一岁,比目前使用的大多数其他语言都要老。
    • @Leonid,当我说要分配最大大小时,我显然是在谈论这个问题。我们知道这里的最大大小是 N,在这种情况下 N=100。对于这样的问题,我看不出还有什么比这更有效的。对数组的访问时间为 O(1)。没有比这更快的了?
    • @me,经过深思熟虑,我意识到这是一个素数的事情。而且由于素数非常稀疏,因此分配最大数是不明智的。除非一个人确切地知道具有这种性质的素数总数是多少(即 p+2 也是素数),否则我不知道素数理论,所以不知道这是否是已知的。我的意思是,如果事先知道需要多少存储空间,那么预先分配一个数组来保存结果是非常有效的。我想这是时间和空间之间的典型权衡。
    • @Nasser 目前假设有无数个孪生素数。见en.wikipedia.org/wiki/Twin_prime
    【解决方案6】:

    这里有另外几个使用NextPrime的替代方案:

    pairs1[pmax_] := Select[Range[pmax], PrimeQ[#] && NextPrime[#] == 2 + # &]
    
    pairs2[pnum_] := Module[{p}, NestList[(p = NextPrime[#];
                          While[p + 2 != (p = NextPrime[p])]; 
                          p - 2) &, 3, pnum]] 
    

    以及对 Reap/Sow 解决方案的修改,可让您指定最大素数:

    pairs3[pmax_] := Module[{k,p},
                       Reap[For[k = 1, (p = Prime[k]) <= pmax, k++,
                            If[PrimeQ[p + 2], Sow[p]]]][[-1, 1]]]
    

    以上是按速度递增的顺序。

    In[4]:= pairs2[10000]//Last//Timing
    Out[4]= {3.48,1261079}
    In[5]:= pairs1[1261079]//Last//Timing
    Out[5]= {6.84,1261079}
    In[6]:= pairs3[1261079]//Last//Timing
    Out[7]= {0.58,1261079}
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2012-01-20
      • 1970-01-01
      • 2019-01-19
      • 1970-01-01
      • 2021-08-04
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多