如何高效产生本原勾股数组？

northwolves

Reap+Sow还可以再快上一点点：

1. Clear["Global`*"];(*清除所有变量*)
   
2. Timing[max=2*10^6;
   
3. SortBy[Sort/@Reap[Do[If[GCD[m,n]==1&&OddQ@(n+m),Sow@{m^2-n^2,2m n,m^2+n^2}],{n,(Sqrt[2max-1]-1)/2},{m,Range[n+1,Sqrt[max-n^2]]}]][[2,1]],Last]]

复制代码

nyy

northwolves 发表于 2023-12-9 15:34
这个还可以快上几倍，勾股数组的斜边的通项公式

ClearAll["Global`*"];(*清除所有变量*)

我发现你比我牛逼。
sortby比sort排序要快很多！
我今天才知道

nyy

northwolves 发表于 2023-12-9 15:34
这个还可以快上几倍，勾股数组的斜边的通项公式

ClearAll["Global`*"];(*清除所有变量*)

1. Clear["Global`*"];(*清除所有变量*)
   
2. Timing[
   
3.     max=2*10^6;(*边长的最大值*)
   
4.     jmax=-1+Sqrt[-1+max];(*由j^2+(j+2)^2<=2*max计算出j的最大值*)
   
5.     aa=Flatten[#,1]&@Table[Sort@{(i*j),(i^2-j^2)/2,(i^2+j^2)/2},(*生成勾股数组,并排序*)
   
6.                            {j,1,jmax,2},(*指定j的范围*)
   
7.                            {i,Select[Range[j+2,Sqrt[2*max-j^2],2],GCD[j,#]==1&]}(*指定i的范围*)
   
8.                           ];
   
9.     bb=SortBy[aa,Last](*按照第三列升序排列,SortBy排序比Sort快很多,原因不知道*)
   
10. ]
    

复制代码

比你的代码快一点点，零点几秒吧

nyy

northwolves 发表于 2023-12-9 15:45
Reap+Sow还可以再快上一点点：

不知道能不能用Compile，这个也许能编译成c语言，然后更快

chyanog

1. cf = Compile[{{max, _Integer}},
   
2.    Module[{jmax, bag},
   
3.     bag = Internal`Bag[Most@{0}];
   
4.     jmax = Sqrt[max - 1] - 1;
   
5.     Do[If[GCD[j, i] == 1,
   
6.       Internal`StuffBag[bag, Floor@{i j, (i^2 - j^2)/2, (i^2 + j^2)/2}, 1]],
   
7.     {j, 1, jmax, 2},
   
8.     {i, j + 2, Sqrt[2 max - j^2], 2}];
   
9.     Partition[Internal`BagPart[bag, All], 3]
   
10.     ]
    
11.    ];
    
12. ans = SortBy[Sort /@ cf[2 10^6], Last]; // AbsoluteTiming
    
13. Length[ans]

复制代码

{0.305319, Null}

318320

nyy

chyanog 发表于 2023-12-14 16:49

Internal`Bag是什么意思？

nyy

chyanog 发表于 2023-12-14 16:49

1. Clear["Global`*"];(*清除所有变量*)
   
2. cf=Compile[{{max, _Integer}},
   
3.     Module[{jmax,aa,i,j},(*指定局部变量*)
   
4.         jmax=-1+Sqrt[-1+max];(*由j^2+(j+2)^2<=2*max计算出j的最大值*)
   
5.         aa=Flatten[#,1]&@Table[Sort@{(i*j),(i^2-j^2)/2,(i^2+j^2)/2},(*生成勾股数组,并排序*)
   
6.                                {j,1,jmax,2},(*指定j的范围*)
   
7.                                {i,Select[Range[j+2,Sqrt[2*max-j^2],2],GCD[j,#]==1&]}(*指定i的范围*)
   
8.                               ]
   
9.     ]
   
10. ]
    
11. bb=SortBy[cf[2*10^6],Last]// AbsoluteTiming;(*按照第三列升序排列,SortBy排序比Sort快很多,原因不知道*)
    

复制代码

我的代码错在什么地方？

nyy

nyy 发表于 2023-12-8 13:42
把你的代码肢解一下，并且添加上注释，这样更容易看明白、更好维护

...

{93,476,485},{44,483,485}
这两组都是本原勾股数，但是排序有问题。
应该是{44,483,485}在{93,476,485}的前面