基于Lagrange四平方定理的整数求解及r4(n)绘图（Matlab/Mathematica）

嘿，作为刚上手Matlab和Mathematica的新手，我来帮你一步步搞定这两个数论问题～下面分别针对两个需求给出两种工具的实现方案，代码都加了注释，好懂易上手！

问题1：找出自然数n的所有四平方和整数表示 #

根据Lagrange四平方定理，每个自然数n都能写成四个整数的平方和：$n = a^2 + b^2 + c^2 + d^2$（a、b、c、d为整数，包括0和负数）。下面是两种工具的实现：

Mathematica 实现 #

Mathematica自带的数论函数能快速帮我们找出所有解，这里用范围限制+筛选的方式，还可以根据需求去重：

n = 10;  (* 替换成你要计算的自然数 *)
maxVal = Floor[Sqrt[n]];  (* 单个平方数最大不超过n，整数范围到sqrt(n)的整数部分 *)

(* 生成所有可能的整数组合，包含正负和0 *)
allCombos = Tuples[Range[-maxVal, maxVal], 4];
(* 筛选满足平方和等于n的组合 *)
validCombos = Select[allCombos, Total[#^2] == n &];

(* 可选：去重——把正负、顺序不同但本质相同的解合并，比如(a,b,c,d)和(-a,b,c,d)视为同一组 *)
uniqueCombos = DeleteDuplicates[Sort /@ validCombos, (Sort[Abs[#1]] == Sort[Abs[#2]]) &];

(* 输出结果 *)
Print["所有符合条件的四平方和组合："];
uniqueCombos

如果只需要非负整数解，把Range[-maxVal, maxVal]改成Range[0, maxVal]即可，去重规则也能简化。

Matlab 实现 #

Matlab里可以通过数组广播生成所有组合，再筛选有效解：

n = 10;  % 替换成目标自然数
max_val = floor(sqrt(n));

% 生成所有可能的整数组合（包含正负和0）
[a, b, c, d] = ndgrid(-max_val:max_val, -max_val:max_val, -max_val:max_val, -max_val:max_val);
% 计算平方和并筛选符合条件的索引
idx = (a.^2 + b.^2 + c.^2 + d.^2) == n;
% 提取有效解
valid_solutions = [a(idx), b(idx), c(idx), d(idx)];

% 可选：去重——将每个解的绝对值排序后去重
sorted_abs = sort(abs(valid_solutions), 2);
[~, unique_idx] = unique(sorted_abs, 'rows');
unique_solutions = valid_solutions(unique_idx, :);

% 输出结果
disp('所有符合条件的四平方和组合：');
disp(unique_solutions);

同样，若只需非负解，把-max_val:max_val改成0:max_val就行。

问题2：绘制r4(n)函数并分析（基于Jacobi定理） #

Jacobi定理明确了四平方和的表示方式数$r_4(n)$的计算公式：

$r_4(n) = 8 \times \sum_{\substack{d|n \ 4 \nmid d}} d$
这里的“表示方式数”是考虑顺序和正负的，比如(1,0,0,0)和(-1,0,0,0)算不同方式，(0,1,0,0)也算另一种，和问题1的去重逻辑不同哦。

Mathematica 实现 #

先实现Jacobi定理的公式，再计算一系列n的$r_4(n)$并绘图：

(* 定义r4(n)函数：基于Jacobi定理 *)
r4[n_] := 8*Total[Select[Divisors[n], Mod[#, 4] != 0 &]];

(* 计算n从1到50的r4(n)值 *)
nRange = Range[1, 50];
r4Values = r4 /@ nRange;

(* 绘制折线图，直观展示规律 *)
ListLinePlot[r4Values, 
  PlotLabel -> "四平方和表示方式数r4(n)",
  AxesLabel -> {"n", "r4(n)"},
  GridLines -> Automatic,
  PlotStyle -> Blue,
  Markers -> Automatic
]

(* 简单规律验证 *)
Print("当n是4的倍数时，r4(n)=r4(n/4)：");
Print("r4(4)=", r4(4), "，r4(1)=", r4(1));
Print("r4(8)=", r4(8), "，r4(2)=", r4(2));

从图里能观察到：$r_4(n)$是积性函数（当n和m互质时，$r_4(nm)=r_4(n)r_4(m)$）；当n是4的倍数时，$r_4(n)=r_4(n/4)$，这些都是Jacobi定理带来的典型规律。

Matlab 实现 #

同样先实现公式，再计算绘图：

% 定义r4(n)函数
function count = r4(n)
    divisors = divisors(n);  % 获取n的所有正约数
    valid_divisors = divisors(mod(divisors, 4) ~= 0);  % 筛选不被4整除的约数
    count = 8 * sum(valid_divisors);
end

% 计算n从1到50的r4(n)值
n_range = 1:50;
r4_values = arrayfun(@r4, n_range);

% 绘制折线图
figure;
plot(n_range, r4_values, 'b-o', 'LineWidth', 1.5);
title('四平方和表示方式数r4(n)');
xlabel('n');
ylabel('r4(n)');
grid on;
set(gca, 'FontSize', 12);

% 简单规律验证
disp('验证当n是4的倍数时，r4(n)=r4(n/4)：');
disp(['r4(4)=', num2str(r4(4)), '，r4(1)=', num2str(r4(1))]);
disp(['r4(8)=', num2str(r4(8)), '，r4(2)=', num2str(r4(2))]);

运行后能看到，$r_4(n)$的增长和n的约数分布直接相关：比如奇质数p的$r_4(p)=8*(1+p)$，$r_4(2)=24$，这些都能从公式推导出来。

内容的提问来源于stack exchange，提问作者Paola Tiranti