a^2+b^2=N

wayne

第一问：

1. MAX = 1000; Union@Flatten@Table[ii^2 + jj^2, {ii, Sqrt[MAX/2]}, {jj, ii, Sqrt[MAX - ii^2]}]

复制代码

第二问： 对于大于2的情况：

1. f[n_] := If[Sqrt[n/2] \[Element] Integers, f[n/2] + 1, Quotient[Table[If[Mod[ii[[1]], 4] == 1, ii[[2]] + 1, 1], {ii, FactorInteger[n]}] /. List -> Times, 2]]

复制代码

〇〇

这些函数式语言真的很简洁

葡萄糖

wayne 发表于 2012-1-20 12:30
第一问：
MAX = 1000; Union@Flatten@Table[ii^2 + jj^2, {ii, Sqrt[MAX/2]}, {jj, ii, Sqrt[MAX - ii^2]}]

第二问：
对于大于2的情况：
f[n_] := If[Sqrt[n/2] \[Element] Integers, f[n/2] + 1, Quotient[Table[If[Mod[ii[[1]], 4] == 1, ii[[2]] + 1, 1], {ii,  FactorInteger[n]}] /. List -> Times, 2]]

第二问好像和OEIS上的结果有点不一样呀
https://oeis.org/A025426

f[n_] := If[Sqrt[n/2] \[Element] Integers, f[n/2] + 1,
   Quotient[
    Table[If[Mod[ii[[1]], 4] == 1, ii[[2]] + 1, 1], {ii,
       FactorInteger[n]}] /. List -> Times, 2]];
Blst = Table[f[k], {k, 100}];
Clst = Table[
   Length[Solve[x^2 + y^2 == k && 0 < x <= y, {x, y}, Integers]],
   {k, 100}];
Transpose[{Blst, Clst}] /. {a_, b_} /; a != b :> Style[{a, b}, Red]

葡萄糖

wayne 发表于 2012-1-19 09:57
已知丢番图方程  a^2+b^2=N，a>=b>0

1）如果只有一组正整数解，这样的N有哪些，给出所有的N
2）对于确定的N，解有几组。可以给出计数方法吗

1. fun[n_Integer] :=
   
2. Module[{decom, \[Delta], \[Epsilon], \[Xi]},
   
3.   decom = Table[
   
4.     Select[FactorInteger[n],
   
5.       Mod[#[[1]], 4] == i && PrimeQ[#[[1]]] &][[All,
   
6.      2]], {i, {2, 1, 3}}];
   
7.   \[Delta] = (Times @[url=home.php?mod=space&uid=6175]@[/url] (1 + decom[[2]]))*(Times @@ ((1 + (-1)^#)/2 & /@
   
8.          decom[[3]]));
   
9.   \[Epsilon] = (Times @@ (((-1)^# + 1)/2 & /@
   
10.         decom[[1]])) (Times @@ ((1 + (-1)^#)/2 & /@
    
11.         decom[[2]])) (Times @@ ((1 + (-1)^#)/2 & /@ decom[[3]]));
    
12.   \[Xi] = (Times @@ ((-1)^# & /@
    
13.         decom[[1]])) (((-1)^Times @@ (1 + decom[[2]]) - 1)/
    
14.       2) (Times @@ ((1 + (-1)^#)/2 & /@ decom[[3]]));
    
15.   Thread[{"\[Delta]", "\[Epsilon]",
    
16.      "\[Xi]"} -> {\[Delta], \[Epsilon], \[Xi]}]]

复制代码

参考链接：
https://mathworld.wolfram.com/SumofSquaresFunction.html

wayne

葡萄糖 发表于 2020-8-18 16:23
第二问好像和OEIS上的结果有点不一样呀
https://oeis.org/A025426

我的计算有误。你的链接很不错，我顺带着找到了这个：
https://mathworld.wolfram.com/GausssCircleProblem.html

wayne

根据链接 https://mathworld.wolfram.com/SumofSquaresFunction.html，得到两种方法：

1. 4 Table[Sum[Sin[(Pi d)/2], {d, Divisors[n]}], {n, 1000}]

复制代码

1. SquaresR[2, Range[1000]]

复制代码