图形学问题--内摆线及格点问题（困惑很久的题目！求高手支招！！）

浪淘沙

内摆线是让一小圆在一大圆内滚动，小圆上的一点所描绘出的轨迹。
令大圆的圆心在原点，并假设起始点在最右端，则其参数式可表为
x(t) = (R-r) cos(t) + r cos((R-r)t/r)
y(y) = (R-r) sin(t) - r sin((R-r)t/r)
其中R和r分别为大小圆的半径。
令C(R,r)为由R和r定义出的内摆线中，其坐标为整数、且存在t值使sin(t)和cos(t)皆为
有理数的点的集合。
令S(R,r) = Σ|x|+|y|对所有(x,y)∈C(R,r)的和。
令T(N) = ΣS(R,r)对所有3≦R≦N，1≦r≦R/2的和。
已知：
C(3, 1) = {(3, 0), (-1, 2), (-1,0), (-1,-2)}。
C(2500, 1000) = {(2500, 0), (772, 2376), (772, -2376), (516, 1792),
              (516, -1792), (500, 0), (68, 504), (68, -504),　(-1356, 1088),
             (-1356, -1088), (-1500, 1000), (-1500, -1000)}。
注意：(-625, 0)并不是集合C(2500, 1000)的一个元素因为对应的sin(t)值不为有理数。
S(3, 1) = (|3| + |0|) + (|-1| + |2|) + (|-1| + |0|) + (|-1| + |-2|) = 10
T(3) = 10
T(10) = 524
T(100) = 580442
T(10^3) = 583108600
请求出T(10^6)。

wayne

参数方程是这个吗？
\[x(t) = (R-r) \cos t + r \cos\frac{(R-r)t}{r}\]
\[y(t) = (R-r) \sin t - r \sin\frac{(R-r)t}{r}\]

浪淘沙

是的，参数方程式这个没错

浪淘沙

内摆线的一些介绍可能有用，大家可以参考下
http://mathworld.wolfram.com/Hypocycloid.html

wayne

根据参数方程，得到\(1):    x^2+y^2 = (R-r)^2 +r^2\)
另外可以得到 \((x-(R-r)\cos t)^2 + (y-(R-r)\sin t)^2 =r^2\)，稍作变换，即 \(2):     R-r = x\cos t +y\sin t\)

因为$R,r$是整数，不妨设$r_1=r,r_2=R-r$,
因为正余弦是有理数，不妨设$\sin t =\frac{2 a b}{a^2+b^2}, \cos t =(a^2 - b^2)/(a^2 + b^2)$,$a,b$为正整数 ,  代入$1),2)$,求解得到：

\[\left\{x\to \frac{a^2 r_2-2 a b r_1-b^2 r_2}{a^2+b^2},y\to \frac{a^2 r_1+2 a b r_2-b^2 r_1}{a^2+b^2}\right\},\left\{x\to \frac{a^2 r_2+2 a b r_1-b^2 r_2}{a^2+b^2},y\to \frac{-a^2 r_1+2 a b r_2+b^2 r_1}{a^2+b^2}\right\}\]
要使x,y是整数，就要保证上面的解可以 化简掉分母。
(未完待续)

浪淘沙

这推导看起来很有滋味，继续学习ing!

但上面的a,b不一定都要正整数吧，sin(t)可以是负值的

浪淘沙

我找一个很久前做过这道题的人帮忙，他翻出个旧文档，说可能有些参数小改了。大家看看这些代码有没参考价值？
被这题堵着真心烦呀。。。

1. #include <vector>
   
2. #include <set>
   
3. #include <numeric>
   
4. #include <iostream>
   
5. 
   
6. using namespace std;
   
7. 
   
8. #define SZ(x) (int)((x).size())
   
9. 
   
10. typedef long long int64;
    
11. 
    
12. template <typename T> T gcd(T x, T y) {for (T t; x; ) t = x, x = y % x, y = t; return y; }
    
13. 
    
14. 
    
15. const int N = 3.6e4;//N要改
    
16. 
    
17. struct node {
    
18.     int64 a, b, c;
    
19. 
    
20.     node(int A, int B, int C) : a(A), b(B), c(C) {}
    
21. } ;
    
22. 
    
23. bool operator < (const node &a, const node &b) {
    
24.     return a.a == b.a ? (a.b < b.b) : (a.a < b.a);
    
25. }
    
26. 
    
27. node operator * (const node &a, const node &b) {
    
28.     return node(a.a * b.a - a.b * b.b, a.a * b.b + a.b * b.a, a.c * b.c);
    
29. }
    
30. 
    
31. int cnt;
    
32. set<node> S;
    
33. 
    
34. void test(int a, int b, int c) {
    
35.     auto p = node(a, b, c);
    
36.     if (S.count(p)) return ;
    
37.     vector<node> v;
    
38.     for (auto q = p; q.c < N; q = q * p) v.push_back(q), S.insert(q);
    
39. 
    
40.     for (int i = 0; i < SZ(v); ++i) {
    
41.         for (int j = i + 1; j < SZ(v); ++j) {
    
42.             int x = gcd(i + 1, j + 1);
    
43.             int pr = N * (i + 1) / (i + j + 2) / v[j].c;
    
44.             for (int t = 1; t <= pr; ++t) {
    
45.                 int64 r = t * v[j].c * (i + 1) / x, R = t * v[j].c * (j + 1) / x;
    
46.                 if (r == 500 && R == 1500)
    
47.                     cout << r << " " << R << " " << (R * v[i].a / v[i].c + r * v[j].a / v[j].c) << " " << (R * v[i].b / v[i].c - r * v[j].b / v[j].c) << endl;
    
48.             }
    
49.         }
    
50.     }
    
51. }
    
52. 
    
53. void gen(int a, int b, int c) {
    
54.     if (c > N) return ;
    
55.     if (a && b) test(a, b, c);
    
56.     ++cnt;
    
57.     int g = (a + b + c) << 1, d, e, f;
    
58.     gen(d = g - a, e = g - b, f = g + c);
    
59.     if (a) g = (a <<= 1) << 1, gen(d - a, e - g, f - g);
    
60.     if (b) g = (b <<= 1) << 1, gen(d - g, e - b, f - g);
    
61. }
    
62. 
    
63. int main() {
    
64.     ios_base::sync_with_stdio(false);
    
65.     gen(3,4,5);//正常使用要改参数？
    
66.     gen(4,3,5);//正常使用要改参数？
    
67.     cout << cnt << endl;
    
68. 
    
69.     return 0;
    
70. }

复制代码

wayne

wayne 发表于 2014-7-8 07:55
根据参数方程，得到\(1):    x^2+y^2 = (R-r)^2 +r^2\)
另外可以得到 \((x-(R-r)\cos t)^2 + (y-(R-r)\sin ...

5# 继续：
$r_1=r,r_2=R-r$, $\sin t =\frac{2 a b}{a^2+b^2}, \cos t =(a^2 - b^2)/(a^2 + b^2)$,$a,b$为整数，且$gcd(|a|,|b|)=1$ ,  代入$1),2)$,求解得到：

\[\left\{x\to \frac{a^2 r_2-2 a b r_1-b^2 r_2}{a^2+b^2},y\to \frac{a^2 r_1+2 a b r_2-b^2 r_1}{a^2+b^2}\right\},\left\{x\to \frac{a^2 r_2+2 a b r_1-b^2 r_2}{a^2+b^2},y\to \frac{-a^2 r_1+2 a b r_2+b^2 r_1}{a^2+b^2}\right\}\]
既然$x,y$为整数，那么$ r_1 x+r_2 y,r_1 x-r_2 y,r_2 x+r_1 y,r_2 x-r_1 y$ 也是整数。于是$a^2+b^2|r_1^2+r_2^2$

mathe

我们知道$x^2+y^2=r_2^2+r_1^2+2r_1r_2\cos((r_1+r_2)h)$其中$r_1h=t$

mathe

我们知道如果cos(u),sin(u),cos(v),sin(v)都是有理数，那么cos(u-v),sin(u-v)也是有理数。由此利用辗转相除法，我们可以得出设$(r_1,r_2)=d$,那么cos(dh),sin(dh)也是有理数。于是对于$(r_1,r_2)=d$的情况，可以通过分析它们互素的情况然后全部乘公共倍数即可