在 Mathematica 中更快地计算大整数的算术函数? [英] Compute the arithmetic functions for large integer in Mathematica faster?
问题描述
设 f 是一个算术函数，A={k1,k2,...,kn} 是递增顺序的整数.
Let f be an arithmetic function and A={k1,k2,...,kn} are integers in increasing order.
现在我想从 k1 开始并将 f(ki) 与 f(k1) 进行比较.如果f(ki)>f(k1)，把ki作为k1.
Now I want to start with k1 and compare f(ki) with f(k1). If f(ki)>f(k1), put ki as k1.
现在从ki开始，比较f(kj)和f(ki)，对于j>i.如果f(kj)>f(ki)，把kj作为ki，重复这个过程.
Now start with ki, and compare f(kj) with f(ki), for j>i. If f(kj)>f(ki), put kj as ki, and repeat this procedure.
最后我们将有一个A的子序列B={L1,...,Lm}，通过这个属性:
At the end we will have a sub sequence B={L1,...,Lm} of A by this property:
L1=k1
L2=ki
L3=kj
...
哪里
f(L(i+1))>f(L(i))，对于任意1<=i<=m-1
例如，设 f 是整数的除数函数.
For example, let f be the divisor function of integers.
我认为应该有一些方法可以比我做的更高效、更快.
I think there should be some way to do more efficient and faster than I did.
你知道如何在 Mathematica 或 Matlab 中为我的目的编写代码.
Do you know how to write a code for my purpose in Mathematica or Matlab.
最好使用 Mathematica.
Mathematica is preferable.
«««««««««««««««««««««««««««««««
«««««««««««««««««««««««««««««««
我用 Mathematica 为这个程序编写了代码，计算 ki 的 f 或大数的集合 B 需要几个小时.
I have written a code for this program with Mathematica, and it take some hours to compute f of ki's or the set B for large numbers.
这里我放了我的代码的一部分，这只是一个示例，我程序中的问题可能比这些更大:
Here I put some part of my code and this is just a sample and the question in my program could be more larger than these:
g 之间的空格是乘积.例如:
the space between g's are product. for example:
g[67757] g[353]=g[67757]*g[353]
g[67757] g[353]=g[67757]*g[353]
«««««««««««««««««««««««««««««««««««
««««««««««««««««««««««««««««««««««««
f[n_] := DivisorSigma[0, n];
g[n_] := Product[Prime[i], {i, 1, PrimePi[n]}];
k1 = g[67757] g[353] g[59] g[19] g[11] g[7] g[5]^2 6^3 2^7;
k2 = g[67757] g[353] g[59] g[19] g[11] g[7] g[5] 6^5 2^7;
k3 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5] 6^4 2^7;
k4 = g[67759] g[349] g[53] g[19] g[11] g[7] g[5] 6^5 2^6;
k5 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5] 6^4 2^8;
k6 = g[67759] g[349] g[53] g[19] g[11] g[7] g[5]^2 6^3 2^7;
k7 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5] 6^5 2^6;
k8 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5] 6^4 2^9;
k9 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5]^2 6^3 2^7;
k10 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5] 6^5 2^7;
k11 = g[67759] g[349] g[53] g[19] g[11] g[7] g[5]^2 6^4 2^6;
k12 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5]^2 6^3 2^8;
k13 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5]^2 6^4 2^6;
k14 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5]^2 6^3 2^9;
k15 = g[67757] g[359] g[53] g[19] g[11] g[7] g[5]^2 6^4 2^7;
k16 = g[67757] g[359] g[53] g[23] g[11] g[7] g[5] 6^4 2^8;
k17 = g[67757] g[359] g[59] g[19] g[11] g[7] g[5] 6^4 2^7;
k18 = g[67757] g[359] g[53] g[23] g[11] g[7] g[5] 6^4 2^9;
k19 = g[67759] g[353] g[53] g[19] g[11] g[7] g[5] 6^4 2^6;
k20 = g[67763] g[347] g[53] g[19] g[11] g[7] g[5] 6^4 2^7;
k = Table[k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16, k17, k18, k19, k20];
i = 1;
count = 0;
For[j = i, j <= 20, j++,
If[f[k[[j]]] - f[k[[i]]] > 0, i = j; Print["k",i];
count = count + 1]];
Print["count= ", count]
«««««««««««««««««««««««««««««««««««
««««««««««««««««««««««««««««««««««««
推荐答案
DivisorSigma 必须分解数字(它不知道它们是如何构造的).您可以通过删除列表的 gcd 来大大加快速度.详细说明:
DivisorSigma has to factor the numbers (it has no idea how they were constructed). You can speed this substantially by removing the gcd of the list. In detail:
将新列表计算为旧列表/gcd.
Compute new list as old list / gcd.
考虑 gcd.
使用一个函数，给定一对因数形式的整数，合并因式分解(因此您的乘积为因数形式).
Use a function that, given a pair of integers in factored form, merges the factorization (so you have their product in factored form).
然后对于缩减列表中的任意两个元素，通过将它们的因式分解与 gcd 的因式分解合并进行比较，并在给定因式分解形式时调用函数来计算除数的数量.最后一个只是指数的乘积，每个指数都增加一.
Then for any two elements in the reduced list, you compare by merging their factorizations each with that of the gcd, and invoking a function to compute the number of divisors when given the factored form. That last is simply the product of the exponents each increased by one.
在代码中:
kgcd = GCD @@ k;
newk = k/kgcd;
gcdfacs = FactorInteger[kgcd];
sumDivisors[faclist_] := Times @@ (1 + faclist[[All, 2]])
mergeFactorLists[fl1_, fl2_] :=
Flatten[GatherBy[Join[fl1, fl2], First] /.
{{p1_Integer,e1_Integer}, {p1_,e2_Integer}} -> {{p1,e1+e2}}, 1]
f2[v1_] := sumDivisors[mergeFactorLists[FactorInteger[v1], gcdfacs]]
这是您的示例，将 f2 应用于 newk 的元素.
Here is your example, with f2 applied to elements of newk.
Timing[i = 1;
count = 0;
For[j = i, j <= 20, j++,
If[f2[newk[[j]]] - f2[newk[[i]]] > 0, i = j; Print["k", i];
count = count + 1]];
Print["count= ", count]]
在评估 In[140]:= k2 期间
During evaluation of In[140]:= k2
在评估 In[140]:= k5 期间
During evaluation of In[140]:= k5
在评估 In[140]:= k7 期间
During evaluation of In[140]:= k7
在评估 In[140]:= k8 期间
During evaluation of In[140]:= k8
在评估 In[140]:= k9 期间
During evaluation of In[140]:= k9
在评估 In[140]:= k10 期间
During evaluation of In[140]:= k10
在评估 In[140]:= k12 期间
During evaluation of In[140]:= k12
在评估 In[140]:= k13 期间
During evaluation of In[140]:= k13
在评估 In[140]:= k14 期间
During evaluation of In[140]:= k14
在评估 In[140]:= k15 期间
During evaluation of In[140]:= k15
在评估 In[140]:= k16 期间
During evaluation of In[140]:= k16
在评估 In[140]:= k17 期间
During evaluation of In[140]:= k17
在评估 In[140]:= k18 期间
During evaluation of In[140]:= k18
在评估 In[140]:= count= 13 期间
During evaluation of In[140]:= count= 13
出[140]= {0.539918，空}
Out[140]= {0.539918, Null}
正如其他人评论的那样，您可能想要使用 SortBy 或者
As others commented, you might instead want to do SortBy or perhaps
sortedk = k[[Ordering[newk, All, f2[#1] < f2[#2] &]]];
--更新 2011-02-01--
--update 2011-02-01--
以下是各种请求的函数，用于对表示为素因数和相应幂的列表的整数进行运算.我们使用效用函数将两个或多个这样的表示相乘"，以便根据上述 g[] 的定义轻松构造它们.
Here are the various requested function, made to operate on integers represented as lists of their prime factors and corresponding powers. We use utility functions to "multiply" two or more such representations, so that they are easily constructed from the definition for g[] above.
对数[fl_] := fl[[All,2]] .日志[fl[[All,1]]]
logarithm[fl_] := fl[[All,2]] . Log[fl[[All,1]]]
divSigma[k_, fax_] := Times @@
((fax[[All, 1]]^(k*(fax[[All, 2]] + 1)) - 1)/(fax[[All, 1]]^k - 1))
mergeFactorLists[f1_,f2_,f3__] :=
mergeFactorLists[mergeFactorLists[f1,f2],f3]
mergeFactorLists[fl1_, fl2_] :=
Flatten[GatherBy[Join[fl1, fl2], First] /.
{{p1_Integer,e1_Integer}, {p1_,e2_Integer}} -> {{p1,e1+e2}}, 1]
eulerPhi[fl_] :=
Times @@ ((fl[[All, 1]] - 1)*fl[[All, 1]]^(fl[[All, 2]] - 1))
我使用 factorlist 的方式类似于上面使用 g[] 的方式，但是是为了获得分解后的列表而不是整数本身.为了方便转换代码，你可以做如下操作.
I use factorlist in a manner similar to use of g[] above, but to obtain the factored lists rather than the integer itself. For ease of converting code, you might do as below.
g[n__] := factorList[n]
那么你将构建 k1 等:
Then you would construct k1 et al as:
k1 = mergeFactorLists[g[67757], g[353], g[59], g[19], g[11], g[7],
g[5, 2], g[4, 3], g[2, 7]];
我注意到使用索引可能会更好，例如k[1]、k[2] 等.这样您就可以存储索引而不是数字(无论是表示为分解列表还是完全展开).这是您的评论或私人电子邮件中的一个问题，我不确定.
I remark that it might be better to use indexing e.g. k[1], k[2], etc. This way you can store the index instead of the number (whether represented as a factored list or fully expanded). This was a concern in either your comments or private email, I'm not sure.
这是一个简短的例子，表明这些功能可能会像宣传的那样工作.
Here is a short example to indicate that the functions might be working as advertised.
在[77]:= 示例 =mergeFactorLists[g[59], g[19], g[11], g[7], g[5, 2], g[4, 3], g[2, 7]]出[77]= {{2, 16}, {3, 9}, {5, 6}, {7, 4}, {11, 3}, {13, 2}, {17,2}, {19, 2}, {23, 1}, {29, 1}, {31, 1}, {37, 1}, {41, 1}, {43,1}, {47, 1}, {53, 1}, {59, 1}}
In[77]:= example = mergeFactorLists[g[59], g[19], g[11], g[7], g[5, 2], g[4, 3], g[2, 7]] Out[77]= {{2, 16}, {3, 9}, {5, 6}, {7, 4}, {11, 3}, {13, 2}, {17, 2}, {19, 2}, {23, 1}, {29, 1}, {31, 1}, {37, 1}, {41, 1}, {43, 1}, {47, 1}, {53, 1}, {59, 1}}
In[83]:= divSigma[2, 例子]出[83]= 8309625653259163198663074449058595410045270294408417958734031\0136565010401600000000
In[83]:= divSigma[2, example] Out[83]= 8309625653259163198663074449058595410045270294408417958734031\ 0136565010401600000000
In[92]:= eulerPhi[示例]出[92]= 30117106786279162451552137484697600000000
In[92]:= eulerPhi[example] Out[92]= 30117106786279162451552137484697600000000
In[95]:=examplenumber = Times@@Map[#[[1]]^#[[2]] &,example]出[95]= 225123336762006539948611826706656256000000
In[95]:= examplenumber = Times @@ Map[#[[1]]^#[[2]] &, example] Out[95]= 225123336762006539948611826706656256000000
In[99]:= DivisorSigma[2, examplenumber]出[99]= 8309625653259163198663074449058595410045270294408417958734031\0136565010401600000000
In[99]:= DivisorSigma[2, examplenumber] Out[99]= 8309625653259163198663074449058595410045270294408417958734031\ 0136565010401600000000
In[100]:= EulerPhi[examplenumber]出[100]= 30117106786279162451552137484697600000000
In[100]:= EulerPhi[examplenumber] Out[100]= 30117106786279162451552137484697600000000
--结束更新--
丹尼尔·利希布劳Wolfram 研究
Daniel Lichtblau Wolfram Research
这篇关于在 Mathematica 中更快地计算大整数的算术函数?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
