目录
理论基础
拉格朗日四平方数和定理
高斯恒等式
操作步骤
分解质因数
求解四平方数
应用高斯恒等式
小结
高斯恒等式输出代码
输出结果
运行结果
怎么把一个大数分解成四个小数的平方和呢？
理论基础
拉格朗日四平方数和定理
每个正整数都可以表示为至多四个正整数的平方和。或者是，每个正整数都可以表示为四个整数的平方和形式。
2 = 0^2 + 0^2 + 1^2 + 1^2
10 = 0^2 + 0^2 +1^2 + 3^2
1234 = 30^2 + 18^2 + 3^2 + 1^2
在历史上也出现过二平方和定理和三平方和定理。证明方式是二次剩余法，或二次互反律，这里不加以证明。
在已知了拉格朗日四平方和定理的基础上，可以通过枚举，累加的方式通过四层循环来寻找满足条件的四个数。
但是在数字比较大的时候，四层循环是很费劲的，十分耗时。那么我们可以用分解质因数的方式，先把大数分解成若干个质因数的乘积，求出质因数的四平方和，然后再通过质因数的四数求出大数大的四数。
高斯恒等式
如果有正整数x和y,
x = x1^2+x2^2+x3^2+x4^2
y = y1^2+y2^2+y3^2+y4^2
那么
x*y = z1^2+z2^2+z3^2+z4^2
z1 = x1*y1+x2*y1+x3*y1+x4*y1
z1 = x1*y2-x2*y1+x3*y4-x4*y3
z1 = x1*y3-x2*y4-x3*y1+x4*y2
z1 = x1*y4+x2*y3-x3*y2-x4*y1
矩阵形式
1 1 1 1
2 -1 4 -3
3 -4 -1 2
4 3 -2 -1
学过高数的同学可以用矩阵乘法推一下。
在有了高斯恒等式之后，我们就可以按照步骤来做了。
操作步骤
第一步分解质因数
分解质因数
分解质因数的方式已经在上一篇文章，讲欧拉定理的时候讲过了。感兴趣的同学可以去看一看。还是一样的方法，先高斯筛质数，然后通过能不能整除，来判断质因数。
%计算一个数的质因数
function e = EULER_E(n)
isnot_prime = zeros(1 , n);%判断是否是合数
prime = [];%质数集
primen = [];%质因数集
if(n == 1)
e = 1;
return;
end
for i = 2:n
if (isnot_prime(i) == 0)
prime(end + 1) = i;
end
for j = 1 : length(prime)
if((i*(prime(j))>n))
break;
end
isnot_prime(i*(prime(j))) = 1;
if(mod(i , prime(j)) == 0)
break;
end
end
end
for k = 1:length(prime)
if(mod(n , prime(k)) == 0)
primen(end + 1) = prime(k);
end
end
e = primen;
下一步，求解质因数的四平方数。
求解四平方数
用四层循环嵌套的方式来求。通过暴力破解的累加法。同时为了节省时间，可以只循环到sqrt(n)，即根号n,同时也可以用字典排列的方式来节省时间，比如第一层循环为i = 1：n,第二层可以是j = 1:i;依此类推。
代码如下：
%计算一个数的四个平方和数（拉格朗日四平方数定理）
function f = FOUR_number_Square_and(n)
for i = 0 : floor(sqrt(n))
for j = 0 : floor(sqrt(n))
for k = 0 : floor(sqrt(n))
for l = 0 : floor(sqrt(n))
if((i^2 + j^2 + k^2 + l^2) == n)
f = [i , j , k , l];
end
end
end
end
end
第三步，应用高斯恒等式。
应用高斯恒等式
这一步内容比较多。首先要判断n = q1^a1+q2^a2+...+qn^an ，各个质因数的次数。
然后，用高斯恒等式，分别先将每个质因数都乘一次，然后再去乘以多余的次数。
代码如下：
%拉格朗日四平方数定理实验
clc;clear;
w = input('n = ');FNSA = [];FNSA_X = [];n = w;prime = [];
primen = EULER_E(w);count = zeros(1 ,length(primen));
for i = 1 : length(primen)
k = 0;
while (mod(n , primen(i)) == 0)
n = n/primen(i);
k = k + 1;
end
count(i) = k;
end
for i = 1 : length(primen)
temp = FOUR_number_Square_and(primen(i));
for j = 1 : 4
FNSA(end+1) = temp(j);
end
end
for i = 1 : length(primen)
if(i == 1)
temp = primen(1);
for k = 1:4
FNSA_X(end+1) = FNSA(k);
tempn = FNSA_X;
end
continue;
end
temp = primen(i)*temp;
FNSA_X(1) = abs(FNSA(4*i-3)*tempn(1)+FNSA(4*i-2)*tempn(2)+FNSA(4*i-1)*tempn(3)+FNSA(4*i)*tempn(4));
FNSA_X(2) = abs(FNSA(4*i-3)*tempn(2)-FNSA(4*i-2)*tempn(1)+FNSA(4*i-1)*tempn(4)-FNSA(4*i)*tempn(3));
FNSA_X(3) = abs(FNSA(4*i-3)*tempn(3)-FNSA(4*i-2)*tempn(4)-FNSA(4*i-1)*tempn(1)+FNSA(4*i)*tempn(2));
FNSA_X(4) = abs(FNSA(4*i-3)*tempn(4)+FNSA(4*i-2)*tempn(3)-FNSA(4*i-1)*tempn(2)-FNSA(4*i)*tempn(1));
tempn = FNSA_X;
end
for j = 1 : length(count)
while(count(j)>1)
count(j) = count(j) - 1;
FNSA_X(1) = abs(FNSA(4*j-3)*tempn(1)+FNSA(4*j-2)*tempn(2)+FNSA(4*j-1)*tempn(3)+FNSA(4*j)*tempn(4));
FNSA_X(2) = abs(FNSA(4*j-3)*tempn(2)-FNSA(4*j-2)*tempn(1)+FNSA(4*j-1)*tempn(4)-FNSA(4*j)*tempn(3));
FNSA_X(3) = abs(FNSA(4*j-3)*tempn(3)-FNSA(4*j-2)*tempn(4)-FNSA(4*j-1)*tempn(1)+FNSA(4*j)*tempn(2));
FNSA_X(4) = abs(FNSA(4*j-3)*tempn(4)+FNSA(4*j-2)*tempn(3)-FNSA(4*j-1)*tempn(2)-FNSA(4*j)*tempn(1));
tempn = FNSA_X;
end
end
temp = 0;
for j = 1:4
temp = temp + FNSA_X(j)^2;
end
if(temp == w)
fprintf('r %d %d %d %d',FNSA_X(1),FNSA_X(2),FNSA_X(3),FNSA_X(4));
else
fprintf('w --%d --%d %d %d &d',temp ,FNSA_X(1),FNSA_X(2),FNSA_X(3),FNSA_X(4));
end
小结
对于高斯恒等式在代码里通过矩阵表示会简洁一些。
完整代码就是上面三个代码。两个函数，一个主程序。
另外在编辑的过程中，写高斯恒等式来费劲啦！！！小编就写了一个程序来输出高斯恒等式，为了减少没及时保存而重打一次的工作量。算是没用的小代码吧。
高斯恒等式输出代码
%打印高斯恒等式
clc;clear;
A = [1,1,1,1;2,-1,4,-3;3,-4,-1,2;4,3,-2,-1];
disp(A);
fprintf('x = ');
for i = 1 : 4
fprintf('x%d^2',i);
if(i ~= 4)
fprintf('+');
else
fprintf('\n');
end
end
fprintf('y = ');
for j = 1 : 4
fprintf('y%d^2',j);
if(j ~= 4)
fprintf('+');
else
fprintf('\n');
end
end
fprintf('x*y = ');
for k = 1 : 4
fprintf('z%d^2',k);
if(k ~= 4)
fprintf('+');
end
end
for i = 1 : 4
for j = 1 : 4
if (j == 1)
fprintf('\n z%d = ',j);
end
if(j ~= 1)
if(A(i , j)>0)
fprintf('+');
else
fprintf('-');
end
end
fprintf('x%d*y%d',j,abs(A(i , j)));
end
end
输出结果
1 1 1 1
2 -1 4 -3
3 -4 -1 2
4 3 -2 -1
x = x1^2+x2^2+x3^2+x4^2
y = y1^2+y2^2+y3^2+y4^2
x*y = z1^2+z2^2+z3^2+z4^2
z1 = x1*y1+x2*y1+x3*y1+x4*y1
z1 = x1*y2-x2*y1+x3*y4-x4*y3
z1 = x1*y3-x2*y4-x3*y1+x4*y2
z1 = x1*y4+x2*y3-x3*y2-x4*y1
对了，补充一下源程序的运行结果：
运行结果
n = 12345678
r 2422 348 2507 271
希望各位有所收获。瑞斯拜。如果有不懂的欢迎留言评论。
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)
