计算从 1 到 k 的 lcm(i, j)、i 和 j 之和的最有效方法是什么？
What is the most performant method to calculate the sum of lcm(i, j), i and j ranging from 1 to k?
问题：
我的解决方案：
这个解决方案有效，但不幸的是，当 k 很大时它真的太慢了（例如 k 可以是一个像 69533550916004 这样大的数字）
var k = 14;
var res = 0;
for(let i = 1; i <= k; i++){
for(let j = 1; j <= k; j++){
res += lcm(i, j);
}
}
console.log(res)
function lcm(a,b){
return (a * b) / gcd(a,b);
}
function gcd(a,b){
if(!b){
return a;
}
return gcd(b, a % b);
}
谢谢
首先简单优化：
如前所述lcm(i,j)与lcm(j,i)
相同
所以你计算了几乎所有总和的热值两次。相反，您可以只计算一次重复热，然后将总和结果加倍。
剩余的非重复热 lcm(i,j) 其中 i==j
lcm(i,i)=i 所以它们形成一个数列 1+2+3+...+k 可以用等式计算： k*(k+1)/2
素数分解 2D LUT 使用 SoE
之前的两种优化结合起来的速度略高于 2 倍（但是复杂度仍然是二次的，所以对于更大的数字来说仍然不够快！！！
对于更高级和更快的方法，您需要足够的内存，以便您可以执行基于 SoE（埃拉托色尼筛法）的素数分解以加快速度显着提高 GCD 和 LCM 操作（bignums 的复杂性更高，常数时间更小）。这个想法是让 k 个列表保存每个数字的质数分解，直到 k。
然后要获得 GCD，您只需比较两个列表并将两个列表中存在的所有素数（及其幂）相乘即可。
LCM 是相似的，但是您将几乎所有素数相乘，您只是使用那些只出现在两个列表中一次的素数。这可以通过不使用列表中的一个数字的质数分解并乘以该数字来进一步增强。
两个列表的搜索都可以在 O(m) 中完成，因为 SoE 将创建已经排序的列表，其中 m 是素数分解 i 或 j ...这比它们的实际值小得多。
如果您想预分配单个 i-th 列表，您可以将不超过 sqrt(i) 的素数作为上限。可以通过以下方式计算出约 10% 的误差：
sqrt(i)/ln(sqrt(i));
所以我用这个只是为了确保我有足够的：
ceil(1.2*sqrt(i)/ln(sqrt(i)))+1;
对于整数，您可以使用 log2（或位数）代替 + 足够大的余量。
然而像这样的 2D LUT 将是大 k ！！！它不适合标准计算机的内存，所以这不适合 bignum k.
这种方法会给你带来更大的提升（测量显示速度比小 k 上的 #1+#2 优化快 ~10 倍 运行 并且随着 k).这种方法会影响 bigints 的复杂性（更好），但它依赖于 bigint 实现。我估计它在我的 bigints 上略低于二次方，但没有尝试 k 离乘法阈值很远，因为这对 bigints 来说仍然太慢了。
使用 SoE
最多 k 1D LUT
这将加速GCD和LCM操作，但没有以前的LUT那么快。从好的方面来说，这需要更少的内存。但是，您当前的限制仍然需要 ~15.8 TB 的存储空间，这在未来可能会实现。然而，这仍然会很慢...
使用 PSLQ
这可能是快速完成此任务的最有前途的方法。对使用
形成的级数进行PSLQ分析
k={ 1,2,3,... }
A(k) = sum_of_lcm_up_to(k)
A = 1,7,28,72,177,303,604,948,1497,...
这可能return你更直接的等式。
IIRC BPP formula for computing Pi 是这样找到的。直到那时，我们还不知道我们可以计算 Pi 的各个数字而不需要以前的数字，也不知道它可以在没有大数算术的情况下计算...
但是我从来没有执行过也没有完全理解 PSLQ，因为它与我的专业领域相去甚远，我天真的理解是，这就像 PCA 对数字系列的各个数字所做的那样，试图找到系列迭代之间的数字关系
这里 C++ 结合了 #1,#2,#3 的代码可以用 bigints 编译，也可以不编译（如果不使用任何第 3 个party libs）你只需用你想使用的任何东西重写 bigint typedef。但要注意 32 位无符号整数（没有 bigints）这会溢出小数字所以限制只有 k<=390 !!!
//---------------------------------------------------------------------------
#include <math.h>
//---------------------------------------------------------------------------
typedef unsigned int bigint; // rewrite int to whatewer bigint datatype you got
//---------------------------------------------------------------------------
// prime decomposition LUT
int pdec_n=0,pdec_siz=0,*pdec_num=NULL,**pdec_LUT=NULL;;
void pdec_exit() // free memory
{
if (pdec_num) delete[] pdec_num;
if (pdec_LUT)
{
if (pdec_LUT[0]) delete[] pdec_LUT[0];
delete[] pdec_LUT;
}
pdec_n=0; pdec_siz=0; pdec_num=NULL; pdec_LUT=NULL;
}
void pdec_init(int n) // pdec_LUT[i][j] holds prime decomposition of i where j={ 0,1,2...pd_num[i]-1 }
{
double x;
int i,j,a;
pdec_exit();
pdec_n=n;
pdec_num=new int[n+1]; // primes in decomposition
pdec_LUT=new int*[n+1]; // decompositions
// compute decompositions size
pdec_num[0]=0;
pdec_num[1]=0;
for (j=0,i=2;i<=n;i++)
{
x=i;
x=sqrt(x); // primes up to sqrt(x)
x=x*log(2.7182818284590452353602874713527)/log(x);
x=ceil(1.2*x)+1; // +20% and round up +1 as a margin for error
a=x;
pdec_num[i]=a;
j+=a;
}
pdec_siz=j+10;
// compute starts of each decomposition
pdec_LUT[0]=new int[j]; // allocate buffer for all decompositions once
pdec_LUT[1]=pdec_LUT[0];
for (j=1,i=2;i<=n;j=i,i++)
{
pdec_LUT[i]=pdec_LUT[j]+pdec_num[j];
pdec_num[j]=0;
}
pdec_num[n]=0;
// SoE prime decomposition
for (i=2;i<=n;i++)
if (pdec_num[i]==0) // is prime?
for (j=i;j<=n;j+=i) // add it to all its multiples
for (a=j;a%i==0;a/=i){ pdec_LUT[j][pdec_num[j]]=i; pdec_num[j]++; } // add also powers
}
//---------------------------------------------------------------------------
bigint pdec_gcd(int a,int b) // fast GCD using pdec_LUT
{
bigint _gcd;
int i,j,ni,nj,*pi,*pj;
pi=pdec_LUT[a]; ni=pdec_num[a];
pj=pdec_LUT[b]; nj=pdec_num[b];
for (_gcd=1,i=0,j=0;(i<ni)&&(j<nj);)
{
a=pi[i];
b=pj[j];
if (a==b){ _gcd*=a; i++; j++; continue; }
if (a<b) i++; else j++;
}
return _gcd;
}
//---------------------------------------------------------------------------
bigint pdec_lcm(int a,int b) // fast LCM using pdec_LUT
{
bigint _lcm;
int i,j,ni,nj,*pi,*pj;
if (pdec_num[a]<pdec_num[b]){ i=a; a=b; b=i; } // b has less multiplications
pi=pdec_LUT[a]; ni=pdec_num[a];
pj=pdec_LUT[b]; nj=pdec_num[b];
for (_lcm=a,i=0,j=0;(i<ni)&&(j<nj);) // a is fully multiplied
{
a=pi[i];
b=pj[j];
if (a==b){ i++; j++; continue; }
if (a< b){ i++; }
if (a> b){ _lcm*=b; j++; } // use primes from b only if not present in a
}
for (;j<nj;j++) _lcm*=pj[j];
return _lcm;
}
//---------------------------------------------------------------------------
int gcd(int a,int b) // slow euclid GCD
{
if(!b) return a;
return gcd(b, a % b);
}
//---------------------------------------------------------------------------
bigint lcm(int a,int b) // slow LCM
{
bigint _lcm;
_lcm=a;
_lcm*=b;
_lcm/=gcd(a,b);
return _lcm;
}
//---------------------------------------------------------------------------
bigint sum_of_lcm0(int k)
{
int i,j;
bigint s;
// sum of lcm(i,j) up to k
for (s=0,i=1;i<=k;i++)
for (j=i+1;j<=k;j++) // do not compute the same terms twice
s+=lcm(i,j);
s+=s; // double the repeating therms
s+=k*(k+1)>>1; // add sum of series: 1+2+3+...+k
return s;
}
//---------------------------------------------------------------------------
bigint sum_of_lcm1(int k)
{
int i,j;
bigint s;
// sum of lcm(i,j) up to k
for (s=0,i=1;i<=k;i++)
for (j=i+1;j<=k;j++)
s+=pdec_lcm(i,j);
s+=s; // double the repeating therms
s+=k*(k+1)>>1; // add sum of series: 1+2+3+...+k
return s;
}
//---------------------------------------------------------------------------
用法：
int k;
bigint s;
k=390; // cant get above this on 32bit uints
pdec_init(k);
s=sum_of_lcm0(k); // #1+#2 k limited by bigint
s=sum_of_lcm1(k); // #1+#2+#3 k limited by bigint and pdec_init
pdec_exit();
这里使用 pdec_init(1000) 和 64 位 bigints 测量更快的方法：
[ 0.507 ms] pdec init 38 KBytes
[ 1.222 ms] k= 100 sum= 18446224
[ 4.613 ms] k= 200 sum= 295280840
[ 10.734 ms] k= 300 sum= 1485491616
[ 21.983 ms] k= 400 sum= 4692090228
[ 38.277 ms] k= 500 sum= 11454503268
[ 52.168 ms] k= 600 sum= 23723529128
[ 67.608 ms] k= 700 sum= 43946379280
[ 83.128 ms] k= 800 sum= 74987181444
[ 106.942 ms] k= 900 sum=120019988208
[ 132.031 ms] k=1000 sum=183011304660
这里是 David Eisenstat 的方法，使用 2D 素数分解 LUT，使用无符号数学并移植到 C++
bigint sum_of_lcm2(int k)
{
int g,*p,n,a,a0,i,sign;
bigint s,f;
// sum of lcm(i,j) up to k
for (s=0,g=1;g<=k;g++)
{
f=k/g;
f=((f+1)*f)>>1;
f*=f;
f*=g;
p=pdec_LUT[g];
n=pdec_num[g];
for (sign=0,a=0,i=0;i<n;i++)
{
a0=a; a=p[i];
if (a0==a) continue; // ignore powers
f*=(a-1);
sign=!sign;
}
if (sign) s-=f; else s+=f;
}
return s;
}
和三围：
[ 62.085 ms] pdec init 18888 KBytes
[ 0.685 ms] k= 1000 sum= 183011304660
[ 1.373 ms] k= 2000 sum= 2926354006228
[ 2.072 ms] k= 3000 sum= 14805042797184
[ 2.781 ms] k= 4000 sum= 46781822900688
[ 3.539 ms] k= 5000 sum= 114218125116052
[ 4.756 ms] k= 6000 sum= 236808736848248
[ 5.223 ms] k= 7000 sum= 438721136165016
[ 6.806 ms] k= 8000 sum= 748454504589824
[ 7.645 ms] k= 9000 sum= 1198740116138576
[ 8.424 ms] k= 10000 sum= 1827127167830060
[ 9.303 ms] k= 11000 sum= 2675130127257140
[ 10.292 ms] k= 12000 sum= 3788610508835184
[ 10.220 ms] k= 13000 sum= 5218119174379260
[ 10.736 ms] k= 14000 sum= 7019284675919704
[ 10.782 ms] k= 15000 sum= 9249622447024312
[ 13.937 ms] k= 16000 sum= 11973761092807672
[ 14.610 ms] k= 17000 sum= 15260078673441888
[ 14.535 ms] k= 18000 sum= 19179420697925096
[ 15.512 ms] k= 19000 sum= 23809411309067844
[ 15.049 ms] k= 20000 sum= 29232792766430776
[ 15.454 ms] k= 21000 sum= 35530729593528460
[ 19.366 ms] k= 22000 sum= 42798909693505916
[ 19.164 ms] k= 23000 sum= 51128789507207824
[ 17.958 ms] k= 24000 sum= 60615343952225808
[ 18.378 ms] k= 25000 sum= 71367596207492300
[ 22.220 ms] k= 26000 sum= 83490183930172416
[ 22.565 ms] k= 27000 sum= 97093514169469316
[ 21.821 ms] k= 28000 sum= 112296054697946652
[ 24.109 ms] k= 29000 sum= 129222336308688284
[ 24.484 ms] k= 30000 sum= 147985909450776556
[ 23.518 ms] k= 31000 sum= 168726314083539724
[ 24.437 ms] k= 32000 sum= 191573248859171168
[ 25.810 ms] k= 33000 sum= 216665462512674040
[ 26.562 ms] k= 34000 sum= 244144733826449780
[ 25.961 ms] k= 35000 sum= 274165636445073084
[ 26.890 ms] k= 36000 sum= 306861019309280956
[ 27.347 ms] k= 37000 sum= 342403256303560068
[ 28.250 ms] k= 38000 sum= 380951178712974600
[ 28.792 ms] k= 39000 sum= 422654544640465536
[ 32.039 ms] k= 40000 sum= 467702202175553120
[ 31.689 ms] k= 41000 sum= 516261770711667188
[ 32.103 ms] k= 42000 sum= 568493102008765720
[ 33.605 ms] k= 43000 sum= 624601491354014924
[ 33.251 ms] k= 44000 sum= 684766511634347904
[ 38.647 ms] k= 45000 sum= 749168090843424836
[ 34.710 ms] k= 46000 sum= 818010181586572760
[ 36.381 ms] k= 47000 sum= 891505184775229460
[ 35.670 ms] k= 48000 sum= 969815208006099344
[ 36.745 ms] k= 49000 sum= 1053194632366715464
[ 39.265 ms] k= 50000 sum= 1141860626298697932
[ 41.430 ms] k= 51000 sum= 1235961937510193376
[ 39.108 ms] k= 52000 sum= 1335793856714706980
[ 39.395 ms] k= 53000 sum= 1441568696460268148
[ 41.918 ms] k= 54000 sum= 1553459252290960280
[ 40.703 ms] k= 55000 sum= 1671762162244218552
[ 41.705 ms] k= 56000 sum= 1796724219557144384
[ 43.305 ms] k= 57000 sum= 1928541735647080528
[ 49.717 ms] k= 58000 sum= 2067474716400847492
[ 43.599 ms] k= 59000 sum= 2213789303836196316
[ 48.271 ms] k= 60000 sum= 2367722439679998864
[ 45.474 ms] k= 61000 sum= 2529554127958058476
[ 47.711 ms] k= 62000 sum= 2699556160225194044
[ 50.947 ms] k= 63000 sum= 2877966575602455352
[ 53.505 ms] k= 64000 sum= 3065104550449828668
[ 54.559 ms] k= 65000 sum= 3261217251285430196
[ 53.066 ms] k= 66000 sum= 3466580213396622284
[ 51.004 ms] k= 67000 sum= 3681481683988726580
[ 52.176 ms] k= 68000 sum= 3906300655499965688
[ 52.597 ms] k= 69000 sum= 4141140842108345608
[ 54.037 ms] k= 70000 sum= 4386474829016834456
[ 55.015 ms] k= 71000 sum= 4642634007103668512
[ 56.036 ms] k= 72000 sum= 4909694650694395656
[ 56.753 ms] k= 73000 sum= 5188180625623422844
[ 56.944 ms] k= 74000 sum= 5478395966863778240
[ 58.519 ms] k= 75000 sum= 5780503390808063112
[ 58.851 ms] k= 76000 sum= 6095030347012208488
[ 59.481 ms] k= 77000 sum= 6422238822380824088
[ 59.821 ms] k= 78000 sum= 6762431127210826000
[ 60.474 ms] k= 79000 sum= 7115946997843513388
[ 62.483 ms] k= 80000 sum= 7483237599224918568
[ 62.773 ms] k= 81000 sum= 7864354728203024788
[ 63.597 ms] k= 82000 sum= 8259936085084243856
[ 63.737 ms] k= 83000 sum= 8670356746307209408
[ 65.268 ms] k= 84000 sum= 9095726603620478688
[ 66.408 ms] k= 85000 sum= 9536717580904610040
[ 66.726 ms] k= 86000 sum= 9993536125103095688
[ 67.575 ms] k= 87000 sum= 10466432130003173312
[ 68.633 ms] k= 88000 sum= 10956004885713459408
[ 68.506 ms] k= 89000 sum= 11462647299842391920
[ 69.760 ms] k= 90000 sum= 11986403043182154584
[ 70.185 ms] k= 91000 sum= 12528192303586301056
[ 71.192 ms] k= 92000 sum= 13088114114334229944
[ 72.118 ms] k= 93000 sum= 13666442859056921956
[ 72.514 ms] k= 94000 sum= 14263784093996812372
[ 73.450 ms] k= 95000 sum= 14880480246958258764
[ 75.154 ms] k= 96000 sum= 15516955556268238668
[ 75.907 ms] k= 97000 sum= 16173675941972720104
[ 76.416 ms] k= 98000 sum= 16851041998775635308
[ 76.931 ms] k= 99000 sum= 17549420431062508128
[ 77.511 ms] k=100000 sum= 18269345553999897648
可以通过从 LUT 中删除素数的幂来加快速度。然而对于大 k 这将不适合内存，所以在这种情况下你需要 1D 素数 LUT 并在 运行 上进行分解，这要慢得多但不需要那么多内存。 IIRC 那里有一些快速的素数分解方法，当使用 FFT 或 NTT 知道最多 sqrt(k) 的素数列表时（除非我在内存中错过匹配它与一些不同的问题）
k 的值很大 -- 二次方时间行不通。这是一个 O(k log log k) 时间算法，与埃拉托色尼筛法相同。
首先，一些数论。我们知道lcm(i, j) = (i*j) / gcd(i, j)。如果我们试图评估
k k
sum sum (i*j),
i=1 j=1
然后我们可以使用分配律重写：
k k k 2
sum sum (i*j) = (sum i) = ((k+1)*k / 2)^2.
i=1 j=1 i=1
我们可以尝试通过对可能的最大公约数求和来挽救这个想法 g:
k k k [k/g] 2
sum sum lcm(i, j) =? sum (1/g)*( sum (g*i) ),
i=1 j=1 g=1 i=1
但这当然是错误的，因为我们最终（例如）lcm(2, 2) 被表示了两次，一次是 g=1，一次是 g=2。为了理顺它，我们求助于莫比乌斯反演，它给了我们正确的公式
k k k [k/g] 2
sum sum lcm(i, j) = sum f(g)*(1/g)*( sum (g*i) ),
i=1 j=1 g=1 i=1
k 2
= sum g*f(g)*(([k/g]+1)*[k/g] / 2) ,
g=1
其中f(g)的定义是
f(g) = product (1-p).
prime p|g
我们可以通过修改埃拉托色尼筛法来批量评估 f。 Python下面3个：
def fast_method(k):
f = [1] * (k + 1)
for i in range(2, k + 1):
if f[i] != 1:
continue
for j in range(i, k + 1, i):
f[j] *= 1 - i
return sum(g * f[g] * ((k // g + 1) * (k // g) // 2) ** 2 for g in range(1, k + 1))
import math
def slow_method(k):
return sum(math.lcm(i, j) for i in range(1, k + 1) for j in range(k + 1))
def test():
for k in range(100):
assert fast_method(k) == slow_method(k), k
test()
问题：
我的解决方案：
这个解决方案有效，但不幸的是，当 k 很大时它真的太慢了（例如 k 可以是一个像 69533550916004 这样大的数字）
var k = 14;
var res = 0;
for(let i = 1; i <= k; i++){
for(let j = 1; j <= k; j++){
res += lcm(i, j);
}
}
console.log(res)
function lcm(a,b){
return (a * b) / gcd(a,b);
}
function gcd(a,b){
if(!b){
return a;
}
return gcd(b, a % b);
}
谢谢
首先简单优化：
如前所述
相同lcm(i,j)与lcm(j,i)所以你计算了几乎所有总和的热值两次。相反，您可以只计算一次重复热，然后将总和结果加倍。
剩余的非重复热
lcm(i,j)其中i==jlcm(i,i)=i所以它们形成一个数列1+2+3+...+k可以用等式计算：k*(k+1)/2素数分解 2D LUT 使用 SoE
之前的两种优化结合起来的速度略高于 2 倍（但是复杂度仍然是二次的，所以对于更大的数字来说仍然不够快！！！
对于更高级和更快的方法，您需要足够的内存，以便您可以执行基于 SoE（埃拉托色尼筛法）的素数分解以加快速度显着提高 GCD 和 LCM 操作（bignums 的复杂性更高，常数时间更小）。这个想法是让
k个列表保存每个数字的质数分解，直到k。然后要获得 GCD，您只需比较两个列表并将两个列表中存在的所有素数（及其幂）相乘即可。
LCM 是相似的，但是您将几乎所有素数相乘，您只是使用那些只出现在两个列表中一次的素数。这可以通过不使用列表中的一个数字的质数分解并乘以该数字来进一步增强。
两个列表的搜索都可以在
O(m)中完成，因为 SoE 将创建已经排序的列表，其中m是素数分解i或j...这比它们的实际值小得多。如果您想预分配单个
i-th列表，您可以将不超过sqrt(i)的素数作为上限。可以通过以下方式计算出约 10% 的误差：sqrt(i)/ln(sqrt(i));所以我用这个只是为了确保我有足够的：
ceil(1.2*sqrt(i)/ln(sqrt(i)))+1;对于整数，您可以使用 log2（或位数）代替 + 足够大的余量。
然而像这样的 2D LUT 将是大
k！！！它不适合标准计算机的内存，所以这不适合 bignumk.这种方法会给你带来更大的提升（测量显示速度比小
k上的 #1+#2 优化快 ~10 倍 运行 并且随着k).这种方法会影响 bigints 的复杂性（更好），但它依赖于 bigint 实现。我估计它在我的 bigints 上略低于二次方，但没有尝试k离乘法阈值很远，因为这对 bigints 来说仍然太慢了。使用 SoE
最多k1D LUT这将加速GCD和LCM操作，但没有以前的LUT那么快。从好的方面来说，这需要更少的内存。但是，您当前的限制仍然需要 ~15.8 TB 的存储空间，这在未来可能会实现。然而，这仍然会很慢...
使用 PSLQ
这可能是快速完成此任务的最有前途的方法。对使用
形成的级数进行PSLQ分析k={ 1,2,3,... } A(k) = sum_of_lcm_up_to(k) A = 1,7,28,72,177,303,604,948,1497,...这可能return你更直接的等式。
IIRC BPP formula for computing Pi 是这样找到的。直到那时，我们还不知道我们可以计算 Pi 的各个数字而不需要以前的数字，也不知道它可以在没有大数算术的情况下计算...
但是我从来没有执行过也没有完全理解 PSLQ，因为它与我的专业领域相去甚远，我天真的理解是，这就像 PCA 对数字系列的各个数字所做的那样，试图找到系列迭代之间的数字关系
这里 C++ 结合了 #1,#2,#3 的代码可以用 bigints 编译，也可以不编译（如果不使用任何第 3 个party libs）你只需用你想使用的任何东西重写 bigint typedef。但要注意 32 位无符号整数（没有 bigints）这会溢出小数字所以限制只有 k<=390 !!!
//---------------------------------------------------------------------------
#include <math.h>
//---------------------------------------------------------------------------
typedef unsigned int bigint; // rewrite int to whatewer bigint datatype you got
//---------------------------------------------------------------------------
// prime decomposition LUT
int pdec_n=0,pdec_siz=0,*pdec_num=NULL,**pdec_LUT=NULL;;
void pdec_exit() // free memory
{
if (pdec_num) delete[] pdec_num;
if (pdec_LUT)
{
if (pdec_LUT[0]) delete[] pdec_LUT[0];
delete[] pdec_LUT;
}
pdec_n=0; pdec_siz=0; pdec_num=NULL; pdec_LUT=NULL;
}
void pdec_init(int n) // pdec_LUT[i][j] holds prime decomposition of i where j={ 0,1,2...pd_num[i]-1 }
{
double x;
int i,j,a;
pdec_exit();
pdec_n=n;
pdec_num=new int[n+1]; // primes in decomposition
pdec_LUT=new int*[n+1]; // decompositions
// compute decompositions size
pdec_num[0]=0;
pdec_num[1]=0;
for (j=0,i=2;i<=n;i++)
{
x=i;
x=sqrt(x); // primes up to sqrt(x)
x=x*log(2.7182818284590452353602874713527)/log(x);
x=ceil(1.2*x)+1; // +20% and round up +1 as a margin for error
a=x;
pdec_num[i]=a;
j+=a;
}
pdec_siz=j+10;
// compute starts of each decomposition
pdec_LUT[0]=new int[j]; // allocate buffer for all decompositions once
pdec_LUT[1]=pdec_LUT[0];
for (j=1,i=2;i<=n;j=i,i++)
{
pdec_LUT[i]=pdec_LUT[j]+pdec_num[j];
pdec_num[j]=0;
}
pdec_num[n]=0;
// SoE prime decomposition
for (i=2;i<=n;i++)
if (pdec_num[i]==0) // is prime?
for (j=i;j<=n;j+=i) // add it to all its multiples
for (a=j;a%i==0;a/=i){ pdec_LUT[j][pdec_num[j]]=i; pdec_num[j]++; } // add also powers
}
//---------------------------------------------------------------------------
bigint pdec_gcd(int a,int b) // fast GCD using pdec_LUT
{
bigint _gcd;
int i,j,ni,nj,*pi,*pj;
pi=pdec_LUT[a]; ni=pdec_num[a];
pj=pdec_LUT[b]; nj=pdec_num[b];
for (_gcd=1,i=0,j=0;(i<ni)&&(j<nj);)
{
a=pi[i];
b=pj[j];
if (a==b){ _gcd*=a; i++; j++; continue; }
if (a<b) i++; else j++;
}
return _gcd;
}
//---------------------------------------------------------------------------
bigint pdec_lcm(int a,int b) // fast LCM using pdec_LUT
{
bigint _lcm;
int i,j,ni,nj,*pi,*pj;
if (pdec_num[a]<pdec_num[b]){ i=a; a=b; b=i; } // b has less multiplications
pi=pdec_LUT[a]; ni=pdec_num[a];
pj=pdec_LUT[b]; nj=pdec_num[b];
for (_lcm=a,i=0,j=0;(i<ni)&&(j<nj);) // a is fully multiplied
{
a=pi[i];
b=pj[j];
if (a==b){ i++; j++; continue; }
if (a< b){ i++; }
if (a> b){ _lcm*=b; j++; } // use primes from b only if not present in a
}
for (;j<nj;j++) _lcm*=pj[j];
return _lcm;
}
//---------------------------------------------------------------------------
int gcd(int a,int b) // slow euclid GCD
{
if(!b) return a;
return gcd(b, a % b);
}
//---------------------------------------------------------------------------
bigint lcm(int a,int b) // slow LCM
{
bigint _lcm;
_lcm=a;
_lcm*=b;
_lcm/=gcd(a,b);
return _lcm;
}
//---------------------------------------------------------------------------
bigint sum_of_lcm0(int k)
{
int i,j;
bigint s;
// sum of lcm(i,j) up to k
for (s=0,i=1;i<=k;i++)
for (j=i+1;j<=k;j++) // do not compute the same terms twice
s+=lcm(i,j);
s+=s; // double the repeating therms
s+=k*(k+1)>>1; // add sum of series: 1+2+3+...+k
return s;
}
//---------------------------------------------------------------------------
bigint sum_of_lcm1(int k)
{
int i,j;
bigint s;
// sum of lcm(i,j) up to k
for (s=0,i=1;i<=k;i++)
for (j=i+1;j<=k;j++)
s+=pdec_lcm(i,j);
s+=s; // double the repeating therms
s+=k*(k+1)>>1; // add sum of series: 1+2+3+...+k
return s;
}
//---------------------------------------------------------------------------
用法：
int k;
bigint s;
k=390; // cant get above this on 32bit uints
pdec_init(k);
s=sum_of_lcm0(k); // #1+#2 k limited by bigint
s=sum_of_lcm1(k); // #1+#2+#3 k limited by bigint and pdec_init
pdec_exit();
这里使用 pdec_init(1000) 和 64 位 bigints 测量更快的方法：
[ 0.507 ms] pdec init 38 KBytes
[ 1.222 ms] k= 100 sum= 18446224
[ 4.613 ms] k= 200 sum= 295280840
[ 10.734 ms] k= 300 sum= 1485491616
[ 21.983 ms] k= 400 sum= 4692090228
[ 38.277 ms] k= 500 sum= 11454503268
[ 52.168 ms] k= 600 sum= 23723529128
[ 67.608 ms] k= 700 sum= 43946379280
[ 83.128 ms] k= 800 sum= 74987181444
[ 106.942 ms] k= 900 sum=120019988208
[ 132.031 ms] k=1000 sum=183011304660
这里是 David Eisenstat 的方法，使用 2D 素数分解 LUT，使用无符号数学并移植到 C++
bigint sum_of_lcm2(int k)
{
int g,*p,n,a,a0,i,sign;
bigint s,f;
// sum of lcm(i,j) up to k
for (s=0,g=1;g<=k;g++)
{
f=k/g;
f=((f+1)*f)>>1;
f*=f;
f*=g;
p=pdec_LUT[g];
n=pdec_num[g];
for (sign=0,a=0,i=0;i<n;i++)
{
a0=a; a=p[i];
if (a0==a) continue; // ignore powers
f*=(a-1);
sign=!sign;
}
if (sign) s-=f; else s+=f;
}
return s;
}
和三围：
[ 62.085 ms] pdec init 18888 KBytes
[ 0.685 ms] k= 1000 sum= 183011304660
[ 1.373 ms] k= 2000 sum= 2926354006228
[ 2.072 ms] k= 3000 sum= 14805042797184
[ 2.781 ms] k= 4000 sum= 46781822900688
[ 3.539 ms] k= 5000 sum= 114218125116052
[ 4.756 ms] k= 6000 sum= 236808736848248
[ 5.223 ms] k= 7000 sum= 438721136165016
[ 6.806 ms] k= 8000 sum= 748454504589824
[ 7.645 ms] k= 9000 sum= 1198740116138576
[ 8.424 ms] k= 10000 sum= 1827127167830060
[ 9.303 ms] k= 11000 sum= 2675130127257140
[ 10.292 ms] k= 12000 sum= 3788610508835184
[ 10.220 ms] k= 13000 sum= 5218119174379260
[ 10.736 ms] k= 14000 sum= 7019284675919704
[ 10.782 ms] k= 15000 sum= 9249622447024312
[ 13.937 ms] k= 16000 sum= 11973761092807672
[ 14.610 ms] k= 17000 sum= 15260078673441888
[ 14.535 ms] k= 18000 sum= 19179420697925096
[ 15.512 ms] k= 19000 sum= 23809411309067844
[ 15.049 ms] k= 20000 sum= 29232792766430776
[ 15.454 ms] k= 21000 sum= 35530729593528460
[ 19.366 ms] k= 22000 sum= 42798909693505916
[ 19.164 ms] k= 23000 sum= 51128789507207824
[ 17.958 ms] k= 24000 sum= 60615343952225808
[ 18.378 ms] k= 25000 sum= 71367596207492300
[ 22.220 ms] k= 26000 sum= 83490183930172416
[ 22.565 ms] k= 27000 sum= 97093514169469316
[ 21.821 ms] k= 28000 sum= 112296054697946652
[ 24.109 ms] k= 29000 sum= 129222336308688284
[ 24.484 ms] k= 30000 sum= 147985909450776556
[ 23.518 ms] k= 31000 sum= 168726314083539724
[ 24.437 ms] k= 32000 sum= 191573248859171168
[ 25.810 ms] k= 33000 sum= 216665462512674040
[ 26.562 ms] k= 34000 sum= 244144733826449780
[ 25.961 ms] k= 35000 sum= 274165636445073084
[ 26.890 ms] k= 36000 sum= 306861019309280956
[ 27.347 ms] k= 37000 sum= 342403256303560068
[ 28.250 ms] k= 38000 sum= 380951178712974600
[ 28.792 ms] k= 39000 sum= 422654544640465536
[ 32.039 ms] k= 40000 sum= 467702202175553120
[ 31.689 ms] k= 41000 sum= 516261770711667188
[ 32.103 ms] k= 42000 sum= 568493102008765720
[ 33.605 ms] k= 43000 sum= 624601491354014924
[ 33.251 ms] k= 44000 sum= 684766511634347904
[ 38.647 ms] k= 45000 sum= 749168090843424836
[ 34.710 ms] k= 46000 sum= 818010181586572760
[ 36.381 ms] k= 47000 sum= 891505184775229460
[ 35.670 ms] k= 48000 sum= 969815208006099344
[ 36.745 ms] k= 49000 sum= 1053194632366715464
[ 39.265 ms] k= 50000 sum= 1141860626298697932
[ 41.430 ms] k= 51000 sum= 1235961937510193376
[ 39.108 ms] k= 52000 sum= 1335793856714706980
[ 39.395 ms] k= 53000 sum= 1441568696460268148
[ 41.918 ms] k= 54000 sum= 1553459252290960280
[ 40.703 ms] k= 55000 sum= 1671762162244218552
[ 41.705 ms] k= 56000 sum= 1796724219557144384
[ 43.305 ms] k= 57000 sum= 1928541735647080528
[ 49.717 ms] k= 58000 sum= 2067474716400847492
[ 43.599 ms] k= 59000 sum= 2213789303836196316
[ 48.271 ms] k= 60000 sum= 2367722439679998864
[ 45.474 ms] k= 61000 sum= 2529554127958058476
[ 47.711 ms] k= 62000 sum= 2699556160225194044
[ 50.947 ms] k= 63000 sum= 2877966575602455352
[ 53.505 ms] k= 64000 sum= 3065104550449828668
[ 54.559 ms] k= 65000 sum= 3261217251285430196
[ 53.066 ms] k= 66000 sum= 3466580213396622284
[ 51.004 ms] k= 67000 sum= 3681481683988726580
[ 52.176 ms] k= 68000 sum= 3906300655499965688
[ 52.597 ms] k= 69000 sum= 4141140842108345608
[ 54.037 ms] k= 70000 sum= 4386474829016834456
[ 55.015 ms] k= 71000 sum= 4642634007103668512
[ 56.036 ms] k= 72000 sum= 4909694650694395656
[ 56.753 ms] k= 73000 sum= 5188180625623422844
[ 56.944 ms] k= 74000 sum= 5478395966863778240
[ 58.519 ms] k= 75000 sum= 5780503390808063112
[ 58.851 ms] k= 76000 sum= 6095030347012208488
[ 59.481 ms] k= 77000 sum= 6422238822380824088
[ 59.821 ms] k= 78000 sum= 6762431127210826000
[ 60.474 ms] k= 79000 sum= 7115946997843513388
[ 62.483 ms] k= 80000 sum= 7483237599224918568
[ 62.773 ms] k= 81000 sum= 7864354728203024788
[ 63.597 ms] k= 82000 sum= 8259936085084243856
[ 63.737 ms] k= 83000 sum= 8670356746307209408
[ 65.268 ms] k= 84000 sum= 9095726603620478688
[ 66.408 ms] k= 85000 sum= 9536717580904610040
[ 66.726 ms] k= 86000 sum= 9993536125103095688
[ 67.575 ms] k= 87000 sum= 10466432130003173312
[ 68.633 ms] k= 88000 sum= 10956004885713459408
[ 68.506 ms] k= 89000 sum= 11462647299842391920
[ 69.760 ms] k= 90000 sum= 11986403043182154584
[ 70.185 ms] k= 91000 sum= 12528192303586301056
[ 71.192 ms] k= 92000 sum= 13088114114334229944
[ 72.118 ms] k= 93000 sum= 13666442859056921956
[ 72.514 ms] k= 94000 sum= 14263784093996812372
[ 73.450 ms] k= 95000 sum= 14880480246958258764
[ 75.154 ms] k= 96000 sum= 15516955556268238668
[ 75.907 ms] k= 97000 sum= 16173675941972720104
[ 76.416 ms] k= 98000 sum= 16851041998775635308
[ 76.931 ms] k= 99000 sum= 17549420431062508128
[ 77.511 ms] k=100000 sum= 18269345553999897648
可以通过从 LUT 中删除素数的幂来加快速度。然而对于大 k 这将不适合内存，所以在这种情况下你需要 1D 素数 LUT 并在 运行 上进行分解，这要慢得多但不需要那么多内存。 IIRC 那里有一些快速的素数分解方法，当使用 FFT 或 NTT 知道最多 sqrt(k) 的素数列表时（除非我在内存中错过匹配它与一些不同的问题）
k 的值很大 -- 二次方时间行不通。这是一个 O(k log log k) 时间算法，与埃拉托色尼筛法相同。
首先，一些数论。我们知道lcm(i, j) = (i*j) / gcd(i, j)。如果我们试图评估
k k
sum sum (i*j),
i=1 j=1
然后我们可以使用分配律重写：
k k k 2
sum sum (i*j) = (sum i) = ((k+1)*k / 2)^2.
i=1 j=1 i=1
我们可以尝试通过对可能的最大公约数求和来挽救这个想法 g:
k k k [k/g] 2
sum sum lcm(i, j) =? sum (1/g)*( sum (g*i) ),
i=1 j=1 g=1 i=1
但这当然是错误的，因为我们最终（例如）lcm(2, 2) 被表示了两次，一次是 g=1，一次是 g=2。为了理顺它，我们求助于莫比乌斯反演，它给了我们正确的公式
k k k [k/g] 2
sum sum lcm(i, j) = sum f(g)*(1/g)*( sum (g*i) ),
i=1 j=1 g=1 i=1
k 2
= sum g*f(g)*(([k/g]+1)*[k/g] / 2) ,
g=1
其中f(g)的定义是
f(g) = product (1-p).
prime p|g
我们可以通过修改埃拉托色尼筛法来批量评估 f。 Python下面3个：
def fast_method(k):
f = [1] * (k + 1)
for i in range(2, k + 1):
if f[i] != 1:
continue
for j in range(i, k + 1, i):
f[j] *= 1 - i
return sum(g * f[g] * ((k // g + 1) * (k // g) // 2) ** 2 for g in range(1, k + 1))
import math
def slow_method(k):
return sum(math.lcm(i, j) for i in range(1, k + 1) for j in range(k + 1))
def test():
for k in range(100):
assert fast_method(k) == slow_method(k), k
test()