只求和的情况下，可以提速一些：f:=(c=0;Do);If,c+=r3(1+k+r)],{r3,1,n/4},{k,Select]&&#>=4r3&]/r3}];c);f@1234
本帖最后由 王守恩 于 2025-3-8 17:45 编辑
基本数量关系这样就可以。\(\D(r_1+r_2)^2-(r_1-r_2)^2=\big(\sqrt{(r_1+r_3)^2-(r_1-r_3)^2\ \ \ \ }+\sqrt{(r_2+r_3)^2-(r_2-r_3)^2\ \ \ \ }\ \big)^2\)
三个圆的半径都是整数时,三个半径之和——由上可得——Union@Flatten@Table——这个算式不好。还是请 northwolves 指正。
三个圆的半径都是整数时,三个半径之和是这样一串数——{9, 18, 27, 36, 45, 49, 54, 63, 72, 81, 98, 147, 169, 196, 245, 294, 338, 343, 392, 441, 507, 676, 845, 882, 961,
1014, 1183, 1323, 1352, 1521, 1764, 1849, 1922, 2205, 2646, 2883, 3087, 3249, 3528, 3698, 3844, 3969, 4805, 5329, 5547, 5766, 6498,——OEIS没有。
补充内容 (2025-3-10 10:33):
这样可能好一些()。Union@Flatten@Table
观察比对了一会儿半径比值的规律，得到一个相对高效的代码
a:=(b=SelectFirst;Range^2/b);
f:=(c=0;Do);If,c+=r3 (1+k+r)],{r3,1,n/4},{k,a}];c);
f
4800044591
这样也可计算大一点的数字了：f=346532068711
根据关系式\(r_1r_2\)是完全平方数，可以设\(r_1=da^2,r_2=db^2\), 其中\((a,b)=1,a\le b\)
带入关系得到\(r_3=\frac{da^2b^2}{(a+b)^2}\),由此得到\((a+b)^2|d\),于是我们可以设\(d=(a+b)^2h\)
所以代入得到\(r_1=(a+b)^2a^2h,r_2=(a+b)^2b^2h,r_3=a^2b^2h\)
这说明\(r_2\)因子分解必然含平方因子。
题目说明的不清楚，假设题目是给定n时，求最大的S(n)
于是我们需要\((a+b)^2b^2h\le n\)时求\(((a+b)^2(a^2+b^2)+a^2b^2)h\)的最大值。
现在如果假设有一个解满足\(b\ge 1\), 我们选择\(a'=(a+b)b-1,b'=1\),h不变，那么同样满足约束条件，对应的S(n)值变为
于是\(a'+b'=(a+b)b\ge a+b, a'b'=ab+b^2-1\gt ab, a'^2+b'^2=(ab+b^2-1)^2+1\gt (2a-1+b^2)^2\ge (a+b^2)^2\ge(a+b)^2\gt a^2+b^2\).
所以S(n)必然变大，也就是我们总可以选择b=1
所以我们需要寻找满足
\((a+1)^2 h\le n\)时，求\( ((a+1)^2(a^2+1)+a^2)h\)的最大值。
记\(c=a+1\ge 2\), 我们要求\(c^2h\le n\)时求\((c^2-c+1)^2 h\)的最大值。
同样容易证明h中含有平方因子时，可以将平方因子转移到c里面结果会更大。所以最大情况h必然不含平方因子。
如果时求所有满足条件的情况之和，那么我们可以记
\(c=(a+b)b\ge 2\)
于是我们可以按c来分类
\(2\le c\le \sqrt{n}\), 在此条件下可以得到\(1\le h\le h_c=\lfloor \frac{n}{c^2}\rfloor\)
然后我们可以对每个范围内的c将它表示为两个不同整数a+b,b的乘积，而且要满足\(b\ge a\),也就是\(b\lt a+b \le 2b\),也就是这两个因子相差不能太远
对每个可能的a,b计算\((a+b)^2(a^2+b^2)+a^2b^2=(a^2+ab+b^2)^2\)并且求和。然后对和再乘上\(\frac{h_c(h_c+1)}2\)。
最后对于不同的c累加上面的结果即可。
13楼代码改进版：
a:=(b=Sqrt@SelectFirst;Range/b);
f:=(c=0;Do,c+=r3 (1+k^2+r^2)],{r3,1,n/4},{k,a}];c);
f//Timing
{0.328125, 4800044591}
参照mathe版主15楼的思路，提速很可观：
f:=Total@Flatten@Table[{(a+b)^2*(a^2)*h,(a+b)^2(b^2)*h,a^2b^2*h},{a,1,n^(1/4)},{b,Select==1&]},{h,1,n/((a+b)a)^2}];f//Timing
{0.25, 480386841332}
本帖最后由 northwolves 于 2025-3-8 19:31 编辑
mathe版主一语道破天机
f:=Sum;1/2 (a^2+a b+b^2)^2 h (1+h),{a,1,n^(1/4)},{b,Select==1&]}];f//Timing
{1.04688, 480583650997741120228449478}