如何把一个4k+1的素数表达成两个整数的平方和

nyy

如何笔算将某些数分拆为两个平方数之和
https://bbs.emath.ac.cn/forum.php?mod=viewthread&tid=15803
(出处: 数学研发论坛)
https://bbs.emath.ac.cn/forum.ph ... 15803&pid=78125

这儿有例子！


作者：东城居士
链接：https://www.zhihu.com/question/3 ... 1900975429948519641
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

Don Zagier 的证明虽然非常巧妙，但并没有告诉我们如何寻找正整数 aaa 和 b.b.b. 是否有一种有效的算法求出 a,ba,ba,b 呢？对于这个问题，同样是在 1990 年 的 The American Mathematical Monthly 97 上，Stan Wagon 介绍了如下简单高效堪称天秀般的算法：首先求出 x2≡−1 (mod p)x^2\equiv-1~(mod ~p)x^2\equiv-1~(mod ~p) 的正整数解 x,x,x, 然后将 ppp 和 xxx 做辗转相除，则不超过 p\sqrt{p}\sqrt{p} 的两个最大的余数即为要求的 aaa 和 b.b.b. 比如 p=111111111111111111111113≡1 (mod 4),p=111111111111111111111113\equiv1~(mod ~4),p=111111111111111111111113\equiv1~(mod ~4), 求出 x2≡−1 (mod p)x^2\equiv-1~(mod ~p)x^2\equiv-1~(mod ~p) 的正整数解为 x=37673906709576022058672.x=37673906709576022058672.x=37673906709576022058672. 将 ppp 和 xxx 做辗转相除，可得余数分别为

这个算法很不错，虽然我不知道原理，发上来共享一下，不知道对于不是素数的整数是否成立！

nyy

1. Clear["Global`*"];(*mathematica11.2,win7(64bit)Clear all variables*)
   
2. (*4k+1素数*)
   
3. p=10^20+129
   
4. (*找到x^2=-1(mod p),随机选择一个a,使得a是非剩余,然后计算得到x*)
   
5. Do[If[JacobiSymbol[k,p]==-1,a=k;Break[]],{k,2,10^10}];
   
6. Print[a];
   
7. x=PowerMod[a,(p-1)/4,p];(*x^2=(a^((p-1)/4))^2=a^((p-1)/2)=-1(mod p)*)
   
8. (*准备辗转相除法*)
   
9. k=0;(*计数,辗转相除法中小于根号p的余数的个数*)
   
10. {a,b}={p,x}(*辗转相除法最初两个数的初始赋值*)
    
11. out={};(*用来保存小于根号p的最大的两个整数*)
    
12. While[Mod[a,b]>0&&k<2,
    
13.     r=Mod[a,b];(*求余数*)
    
14.     {a,b}={b,r};(*得到新的被除数与除数*)
    
15.     If[r<Sqrt[p],out=Append[out,r];k=k+1]
    
16. ]
    
17. Print[{out,p,Total[out^2]-p}](*打印出两个数与素数,并且检验a^2+b^2-p是否等于零*)
    

复制代码

代码写出来了，虽然是比较丑的代码，但是能用！

nyy

nyy 发表于 2025-5-21 13:21
代码写出来了，虽然是比较丑的代码，但是能用！

{{
8734299514697096415092324422262172643487651173596610058044551075223059534652741558693062858574717874,
4869498124813486982231377956355473573990401367667605534476012902371537410710789858816644285044662191
},
10^200+357,
0}

10^200+357写成两个整数的平方和

nyy

nyy 发表于 2025-5-21 13:25
{{
87342995146970964150923244222621726434876511735966100580445510752230595346527415586930628585747 ...

有价值的是上面那个算法！
那个算法比较神奇！

nyy

nyy 发表于 2025-5-21 22:11
有价值的是上面那个算法！
那个算法比较神奇！

这个算法思想很赞，你的代码可以优化，比如使用NestWhile

秀出你的代码！

northwolves

1. f[p_] := (a = NestWhile[1 + # &, 2, JacobiSymbol[#, p] != -1 &];
   
2.    b = PowerMod[a, (p - 1)/4, p];
   
3.    NestWhileList[{Mod[#[[2]], #[[1]]], #[[1]]} &,Sort[{p, b}], #[[2]]^2 > p &][[-1]]);
   
4. f[10^200 + 357]

复制代码

{4869498124813486982231377956355473573990401367667605534476012902371537410710789858816644285044662191,8734299514697096415092324422262172643487651173596610058044551075223059534652741558693062858574717874}

nyy

Let p be a prime number of the form 4k+1. Fermat's theorem asserts that p is a sum of two squares, p=x2+y2.

There are different proofs of this statement (descent, Gaussian integers,...). And recently I've learned there is the following explicit formula (due to Gauss): x=12(2kk)(modp), y=(2k)!x(modp) (|x|,|y|<p/2). But how to prove it?

https://math.stackexchange.com/q ... fermats-4k1-theorem

这个算法直接给出解了！但是非常难算，不实用！

nyy

northwolves 发表于 2025-5-23 09:26
{4869498124813486982231377956355473573990401367667605534476012902371537410710789858816644285044662 ...

1. NestWhile[(1 + #) &, 2, (JacobiSymbol[#, p] != -1) &]

复制代码

个人认为，你的代码应该写成这样。
后面的是一个函数，所以这行代码应该有两个&
按照标准写法是这样，不知道你的怎么回事。
看别人的代码费劲。
尤其是一行代码，且没注释没缩进的代码

northwolves

一句代码

1. f[p_]:=NestWhileList[{Mod[#[[2]],#[[1]]],#[[1]]}&,Sort[{p,PowerMod[NestWhile[1+#&,2,JacobiSymbol[#,p]!=-1&],(p-1)/4,p]}],#[[2]]^2>p&][[-1]];

复制代码

nyy

How to prove that all primes of the form  4k+1  can be represented by the sum of two squares in only one way regardless of the order?

https://math.stackexchange.com/q ... ented-by-the-sum-of