cody习题：寻找圆上的整数坐标点

bainhome

向雷锋学习：替祁某人代发题目，我自己在熬夜赶结题项目报告，还没来得及想，不过直觉感觉挺有意思，我比较喜欢这种设问简单但内涵深意的代码题。
cody原题见这里。
大意就是寻找一个以原点为圆心的圆上有多少个坐标点二维坐标值都是整数（正负零都算）。
例如test suite里半径为5时，应该是12个满足要求的坐标：

1. (0, 5) and (0, -5)
   
2. (5, 0) and (-5, 0)
   
3. (2, 1) and (2, -1)
   
4. (-2, 1) and (-2, -1)
   
5. (1, 2) and (1, -2)
   
6. (-1, 2) and (-1, -2)

复制代码

但是由于一些test的suite非常大，枚举恐怕运算时间非常长...
列出test suite：

1. %% 1
   
2. assert(isequal(circle_points(1),4))
   
3. %% 2
   
4. assert(isequal(circle_points(3),4))
   
5. %% 3
   
6. assert(isequal(circle_points(5),12))
   
7. %% 4
   
8. assert(isequal(circle_points(65),36))
   
9. %% 5
   
10. assert(isequal(circle_points(64090),324))
    
11. %% 6
    
12. assert(isequal(circle_points(2417899275),20))
    
13. %% 7
    
14. assert(isequal(circle_points(31432690549),8748))
    
15. %% 8
    
16. assert(isequal(circle_points(1021090952484265),236196))
    
17. %% 9
    
18. assert(isequal(circle_points(6095127531752228),78732))
    
19. %% 10
    
20. assert(isequal(circle_points(5*circle_points(630209)),12))
    
21. %% 11
    
22. assert(isequal(circle_points(65*circle_points(236009)),36))
    
23. %% 12
    
24. filetext = fileread('circle_points.m');
    
25. assert(isempty(strfind(filetext, 'switch')))
    
26. assert(isempty(strfind(filetext, 'case')))
    

复制代码

test12貌似不让用switch-case流程。

winner245

不知道是我理解有误还是 test suite 写得有误，半径为 5 时，x^2+y^2 = 5^2，而test suite给出的是 x^2+y^2 = 5 的情形。如果是按照 x^2+y^2 = 5^2，那么12个点应该是：(3,4,5), (4,3,5) 在4个象限的镜像，故有8个点。再加上圆周与坐标轴的 4 个交点，故总共12个点。

nwcwww

抛砖引玉：

1. function ans = circle_points(r)
   
2.     factor(r);
   
3.     histc(ans, unique(ans(rem(ans-4, 4)==1)));
   
4.     prod(2*ans+1)*4;
   
5. end

复制代码

lin2014

实际上这是个求整数解的简单方程。
用maple验证一下，当r = 1000时，应为28组解；
> S := isolve(x^2+y^2 = 1000^2);
{x = -1000, y = 0}
{x = -960, y = -280}
{x = -960, y = 280}
{x = -936, y = -352}
{x = -936, y = 352}
{x = -800, y = -600}
{x = -800, y = 600}
{x = -600, y = -800}
{x = -600, y = 800}
{x = -352, y = -936}
{x = -352, y = 936}
{x = -280, y = -960}
{x = -280, y = 960}
{x = 0, y = -1000}
{x = 0, y = 1000}
{x = 280, y = -960}
{x = 280, y = 960}
{x = 352, y = -936}
{x = 352, y = 936}
{x = 600, y = -800}
{x = 600, y = 800}
{x = 800, y = -600}
{x = 800, y = 600}
{x = 936, y = -352}
{x = 936, y = 352}
{x = 960, y = -280}
{x = 960, y = 280}
{x = 1000, y = 0}
> nops([S]); #计算个数
28

当 r= 10000时，应有36组解；
         100000（10万）时有44个解。
用maple很快就可以得出答案。

lin2014

nwcwww 发表于 2014-12-6 02:22
抛砖引玉：

算法很精妙，运行得也很快，能解释一下就更好了。
把所有解都列出来时的速度怎么样？
试试 r = 10^8时运行时间（列出所有解）。

rocwoods

nwcwww 发表于 2014-12-6 02:22
抛砖引玉：

nwc老兄抛出的是玉，算法已经很难提升了，而且size也很小了，别人只能抛砖啦！我来一个accumarray的。

1. function ans = circle_points(r)
   
2. factor(r);
   
3. ans(rem(ans-4, 4)==1);
   
4. accumarray(ans(:),1,[],[],[],1);
   
5. prod(2*nonzeros(ans(5:end))+1)*4;
   
6. end

复制代码

qibbxxt

1. 
   

复制代码

该问题似乎是寻找圆上的整数点的个数，由于测试的例子中，半径比较大，因此采用枚举所有整数点的方法不可取

因此就得从算法角度去考虑，和方程的整数解，复数或者是直角三角形相关，网络上有很多相关的算法

细心点，可以发现cody原题目下面的评论

1. Tim on 28 Mar 2013
   
2. 
   
3. This is derived from the Mathematica algorithm for sequence A046080 at oeis.org (arrived at from A046109). Instead of factoring (because of the large integers) it checks for divisibility by the various primes (up to 325643, which is enough to handle the test set).

复制代码

在oeis里面输入A046080,即  http://oeis.org/A046080

可以看到公式

Let n = 2^e_2 * product_i p_i^f_i * product_j q_j^g_j where p_i == 1 mod 4, q_j == 3 mod 4; then a(n) = (1/2)*(product_i (2*f_i + 1) - 1).

剩下的就是编写程序，一些用能的函数
1. factor: 求质因数
2. histc/accumarray/bsxfun(@eq,...)/hist： 统计个数,2014b的版本可能是histcounts


大家有什么新的想法，希望一起讨论