关于将一个数分解成四个数平方和的算法matlab
目录
理论基础
拉格朗日四平方数和定理
高斯恒等式
操作步骤
分解质因数
求解四平方数
应用高斯恒等式
小结
高斯恒等式输出代码
输出结果
运行结果
怎么把一个大数分解成四个小数的平方和呢？
理论基础
拉格朗日四平方数和定理
每个正整数都可以表示为至多四个正整数的平方和。或者是，每个正整数都可以表示为四个整数的平方和形式。
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:nif (isnot_prime(i) == 0)prime(end + 1) = i;endfor j = 1 : length(prime)if((i*(prime(j))>n))break;endisnot_prime(i*(prime(j))) = 1;if(mod(i , prime(j)) == 0)break;endend
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];endendendend
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;endcount(i) = k;
end
for i = 1 : length(primen)temp = FOUR_number_Square_and(primen(i));for j = 1 : 4FNSA(end+1) = temp(j);end
end
for i = 1 : length(primen)if(i == 1)temp = primen(1);for k = 1:4FNSA_X(end+1) = FNSA(k);tempn = FNSA_X;endcontinue;endtemp = 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:4temp = 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));
elsefprintf('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 : 4fprintf('x%d^2',i);if(i ~= 4)fprintf('+');elsefprintf('\n');end
end
fprintf('y = ');
for j = 1 : 4fprintf('y%d^2',j);if(j ~= 4)fprintf('+');elsefprintf('\n');end
endfprintf('x*y = ');for k = 1 : 4fprintf('z%d^2',k);if(k ~= 4)fprintf('+');endend
for i = 1 : 4for j = 1 : 4if (j == 1)fprintf('\n z%d = ',j);endif(j ~= 1)if(A(i , j)>0)fprintf('+');elsefprintf('-');endendfprintf('x%d*y%d',j,abs(A(i , j)));end
end输出结果
1 1 1 12 -1 4 -33 -4 -1 24 3 -2 -1x = 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希望各位有所收获。瑞斯拜。如果有不懂的欢迎留言评论。
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！