[C++] 22695번 Spirograph

https://www.acmicpc.net/problem/22695

---

 

25/10/21

 

 

수치 적분 문제인데, 그냥 하면 된다.

---

 

문제 접근 방식:

 

 

일단 나는 Hypotrochoid에 대한 식을 위키에서 찾아봤음을 먼저 이야기해야 한다.(때문에 난이도에 대한 평가가 정확하지 않을 수 있다.)

 

$$ \begin{align} x(\theta) &= (R-r)\cos \theta + d\cos\left( \frac{R-r}{r} \theta\right) \\ y(\theta) &= (R-r)\sin \theta - d\sin\left( \frac{R-r}{r}\theta \right) \end{align} $$

 

우리는 해당 식을 통해 Arc length를 구할 수 있다.

 

Arc length를 구하기 위해 $\theta$에 대해 적분을 해보자.

 

$$ \begin{align} \int \sqrt{ \left( \frac{dx}{d\theta} \right)^{2} + \left( \frac{dy}{d\theta} \right)^{2} }d\theta &= \int \sqrt{ (R-r)^{2} + \frac{d^{2}(R-r)^{2}}{r^{2}} - \frac{2d(R-r)^{2}}{r}\cos\left( \frac{R}{r}\theta \right) }d\theta \\ &= \int \sqrt{ a-b\cos(c\theta) }d\theta \end{align} $$

 

이때, 각각의 계수는 $a, b, c$로 치환하였다.

 

해당 모양의 적분은 타원 적분으로, Closed Form이 없다는 것이 자명하게 밝혀져있다. 따라서 수치적분을 해야하는데, 범위를 살펴볼 필요가 있다.

 

범위는 $0$부터 $2\pi \times \frac{\text{LCM}(R, r)}{R}$이다. 이때, $d = 0$인 경우, 그냥 원이 나오므로 범위가 $2\pi$까지이다.

 

다양한 수치적분 방법이 있지만, 나의 경우는 Trapezoidal Rule을 사용하여 해당 적분 방법으로 구현하였다.

---

아래는 내가 위의 접근 방식과 같이 작성한 C++ 코드이다. 더보기를 누르면 확인할 수 있다.

// 22695번 Spirograph
// 미적분학, 수치 해석
#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;
#define fastio ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
#define endl '\n'
#define PI 3.14159265358979323846
typedef long double ld;

ld iters = 1000000.0;

int gcd(int a, int b){
    if (b == 0) return a;
    return gcd(b, a % b);
}

int lcm(int a, int b){
    return a*b/gcd(a, b);
}

ld f(ld a, ld b, ld c, ld x){
    ld k = a - b*cosl(c*x);
    if (k < 0) return 0.0L;
    return sqrtl(k);
}

ld trapzoid(ld R, ld r, ld d){
    ld a, b, c;
    a = (R-r)*(R-r) + d*d*(R-r)*(R-r)/(r*r);
    b = 2*d*(R-r)*(R-r)/r;
    c = R/r;

    ld ans = 0;
    ld total = 2*PI;
    if (d != 0) total *= lcm((int)R, (int)r)/(int)R;
    ld x = 0, add = total/iters;
    ans += f(a, b, c, x);
    for (int i = 1; i <= (int)(iters-1); ++i){
        x += add;
        ans += 2*f(a, b, c, x);
    }
    ans += f(a, b, c, x);
    ans *= (add/2);

    return ans;
}

int main(void){
    fastio

    ld R, r, d;
    cout << fixed << setprecision(20);
    while (1){
        cin >> R >> r >> d;
        if (R == 0 && r == 0 && d == 0) return 0;
        cout << trapzoid(R, r, d) << endl;
    }
    return 0;
}