如何求出1-n以内勾股数之和

wayne

20# G-Spider
thank you ,gelivable!

zeroieme

14# wayne

    for(b=1;b*b<n/2.0;b++)
    for(a=b+1;a*a<n;a++)
几乎是全搜索了

if(gcd(a-b,2*b)==1)
……
s=n/t
sum+=s*(s+1)*a*(a+b)
产生出互素基本勾股数 a、  b、 c，然后用高斯公式求整个勾股数族 k a、 k b、k c 的和

分析勾股数公式
$a = 2mn$
$b = m^2— n^2$
$c = m^2  + n^2$
要 a、  b、 c互素，必然 m、n互素而且一奇一偶。
以4的倍数作为a。然后将$a/2$，分解质因数，组合出m、n。

zeroieme

更直接的
给出m (m<N)。则n [1,s]，s=[$sqrt(N-m^2)$]。
Sam(a)=m s (s+1)
Sam(b+c)=2*m*m *s

复杂度可以减到O(N)


——————
应当m (m<$sqrt(N)$)
复杂度O($sqrt(N)$)

wayne

22# zeroieme
我的算法时间复杂度就是O(N)

zeroieme

for(b=1;b*b<n/2.0;b++)
    for(a=b+1;a*a<n;a++)
$(sqrt(N))^2 -> N $

wayne

25# zeroieme
  ,严格来说是 O(${N}/{2\sqrt{2}}$)

zeroieme

这样的话由m导出
$a=2mn$
$b=m^2 - n^2$
$c=m^2  + n^2$
应当是O($sqrt(N)$)

另外O()是正比于某个增长度，不应当带常数项。

zeroieme

错了，不能产生全部衍生勾股数

wayne

14# wayne
该代码还可以改进， 循环内部可以少一层if判断，在unsigned long int范围内，速度已经无可挑剔了：
如果想查看所有的勾股数，去掉注释即可， 建议用重定向的方式运行程序，比如  
a.exe 1000000>a.txt   将产生75.6M的文件。

1. #include <stdio.h>
   
2. typedef unsigned long int DW;
   
3. inline DW gcd(DW m,DW n){DW t;  while(m!=0) t=m,m=n%m,n=t;return n;}  
   
4. DW Test(DW n){
   
5.     DW t,s,ii,cnt=1,sum=0,c=1;
   
6.     DW a,b;
   
7.     for(b=1;b*b<=n/2;b++)
   
8.         for(a=b+1;t=a*a+b*b,t<=n;a++){
   
9.                         if(gcd(a-b,2*b)==1){
   
10.                                 s=n/t;       
    
11.                                 //for(ii=1;ii<=s;ii++)
    
12.                                 //printf("%5u.%-5u\t%u,%u,%u\n",c++,cnt,ii*(a*a-b*b),2*ii*a*b,ii*t);
    
13.                                 //cnt++;
    
14.                                 sum+=s*(s+1)*a*(a+b);
    
15.                           }
    
16.                 }      
    
17.     return sum;
    
18. }
    
19. 
    
20. int main(int c ,char **v)
    
21. {   if (c==2) {
    
22.                 printf("total of triples less than %s :\t%lu\n",v[1],Test(atoi(v[1])));
    
23.     }
    
24.     return 0;
    
25. }

复制代码

zeroieme

29# wayne

挑根小刺
if(gcd(a-b,2*b)==1)

既然a-b与 2*b互质，那么a、b不能同奇偶
for(a=b+1;t=a*a+b*b,t<=n;a++)
改成
for(a=b+1;t=a*a+b*b,t<=n;a+=2)
是否更快点