ProjEuler[39]:找勾股数

jimogsh

If p is the perimeter of a right angle triangle with integral length sides, {a,b,c}, there are exactly three solutions for p = 120.

{20,48,52}, {24,45,51}, {30,40,50}

For which value of p ≤ 1000, is the number of solutions maximised?

120可以写为三组勾股数（整数）之和，分别是{20,48,52}, {24,45,51}, {30,40,50}。要求在不超过1000的数字中求出可以写为最多组勾股数之和的数。

jimogsh

1. Timing[data = Range[500]^2;
   
2. lis = {};
   
3. For[m = 2, m < 500, m++,
   
4.   For[n = 2, n <= m, n++,
   
5.    If[MemberQ[data, data[[m]] + data[[n]]],
   
6.     lis = Append[lis, {m, n, Sqrt[m^2 + n^2]}]]]];
   
7. Total /@ lis // Commonest]

复制代码

要算10秒左右。感觉不是很好。找勾股数有什么好办法？

FreddyMusic

My notebook should be faster.

JustinChen

1. f[lim_] := Module[{m = 1, n = 2, lis = {}},
   
2. While[m <= Floor[lim/3],
   
3. While[n <= Floor[(lim - m)/2],
   
4. If[IntegerQ[Sqrt[m^2 + n^2]] && m + n + Sqrt[m^2 + n^2] <= lim,
   
5. AppendTo[lis, m + n + Sqrt[m^2 + n^2]]];
   
6. n++;];
   
7. m++; n = m + 1];
   
8. lis // Commonest]

复制代码

1. 
   
2. Timing[f[1000]]
   

复制代码

试试看~

chyanog

本来想能不用Reduce就不用，由于它不够“底层”，不过的确强大。下面两种做法都用到了Reduce,其中第二个效率不错,当范围增大时尤为明显。

1. 
   
2. Timing[sol =
   
3.   Reduce[3 < x + y + z <= 1000 && x^2 + y^2 == z^2 && x > 0 && y > 0 &&
   
4.      z > 0, {x, y, z}, Integers];
   
5. x + y + z //. (Cases[sol, Except[False]] /. {And -> List, Or -> List,
   
6.        Equal -> Rule}) // Commonest]
   

复制代码

{7.156, {840}}

1. 
   
2. Timing[GGS[max_] := Module[{x, y, r,s},
   
3.    r = Table[
   
4.      Reduce[x^2 + y^2 == z^2 && 0 < x <y <= max, {x, y},
   
5.       Integers], {z, max}];
   
6.    s = {x, y} /.
   
7.       DeleteCases[r /. {And -> List, Or -> List, Equal -> Rule},
   
8.        False];
   
9.    Sort@Partition[Flatten[s], 2] /. {x_, y_} -> {x, y, Sqrt[x^2 + y^2]}
   
10.    ];
    
11. Commonest[Tr /@ GGS[500]]]
    

复制代码

{1.578, {840}}

guocong89

问题的瓶颈在于如何记录已经发现的勾股数

1. Module[
   
2.   {found, cnt, i, j, k},
   
3.   found = {};
   
4.   cnt = Table[0, {1000}];
   
5.   For[i = 1, i < 500, i++,
   
6.    For[j = 1, j < i, j++,
   
7.     If[2 i j + 2 i i > 1000, Break[]];
   
8.     For[k = 1, k < 1000, k++,
   
9.      If[k (2 i j + 2 i i) > 1000, Break[]];
   
10.      If[!
    
11.        MemberQ[found,
    
12.         Hash[{2 i j  k, k (i i - j j), k (i i + j j)}]],
    
13.       cnt[[k (2 i j + 2 i i)]]++;
    
14.       found =
    
15.        Append[found, Hash[{2 i j  k, k (i i - j j), k (i i + j j)}]]
    
16.       ]
    
17.      ]
    
18.     ]
    
19.    ];
    
20.   j = 1;
    
21.   For[i = 2, i < 1001, i++, If[cnt[[j]] < cnt[[i]], j = i]];
    
22.   j
    
23.   ] // Timing

复制代码

时间： {0.030996, 840}

qibbxxt

matlab计算的效率要高一些的

1. clear;clc;close all
   
2. tic;
   
3. [a,b]=meshgrid(1:500);
   
4. c=a.^2+b.^2;
   
5. d=sqrt(c);
   
6. f=round(d);
   
7. e=f==d;
   
8. k=find(e);
   
9. p=a(e)+b(e)+f(e);
   
10. q=accumarray(p,1);
    
11. m=find(q==max(q));
    
12. ij=find(p==m);
    
13. re=[a(k(ij)) b(k(ij)) f(k(ij))];
    
14. re=unique(sort(re,2),'rows');
    
15. fprintf('%d^2+%d^2=%d^2\n',re');
    
16. toc

复制代码

1. m =
   
2. 
   
3.    840
   
4. 
   
5. 40^2+399^2=401^2
   
6. 56^2+390^2=394^2
   
7. 105^2+360^2=375^2
   
8. 120^2+350^2=370^2
   
9. 140^2+336^2=364^2
   
10. 168^2+315^2=357^2
    
11. 210^2+280^2=350^2
    
12. 240^2+252^2=348^2
    
13. Elapsed time is 0.025446 seconds.

复制代码

waynebuaa

此题稍加分析即可笔算出来