pellEqu  펠 방정식 x2 - M y2 = c 의 최소해를 얻는다.
         다만, M 는 평방수가 아닌 정정수, c = 1, -1.

int pellEqu(long *x1, long *y1, long *x2, long *y2, int m);

x1, y1  (출력) 펠 방정식  x2 - M y 2 = 1 에 대한 최소해
x2, y2  (출력) 펠 방정식  x2 - M y 2 = -1 에 대한 최소해
m       (입력) 펠 방정식의 계수 M (1,4,9,16...이외의 비평방수)

c = 1 또는 -1 의 양쪽 모두에 대한 최소해를 얻을 수 있었을 때는 1,
c = -1 에 대한 최소해가 존재하지 않았을 때는 0.
c = 1 에 대한 최소해도 c = 1 에 대한 최소해도 얻을 수 없었을 때는 -1.

int pellEqu(long *x1, long *y1, long *x2, long *y2, int m)
{
    long p0, q0, p1, q1, p, q;
    int  a0, a;
    int  s, t, k;
    
    if (m <= 1) return -1;

    *x1 = *y1 = *x2 = *y2 = 0;
    p0 = 1, q0 = 0;

    /* get m^(1/2) */
    a0 = m, k = 1;               
    while (k < a0) k <<= 1, a0 >>= 1;
    do a0 = k, k = (m / k + k) >> 1; 
    while (k < a0);
    if (a0 * a0 == m) return -1;

    s = 0, t = 1, k = 0;
    p = a = a0, q = 1;
    for ( ; ; ) {
        s = a * t - s;
        t = (m - s*s) /t;
        if (t == 1) break;
        k++;
        a = (a0 + s) / t;
        p1 = p, q1 = q;
        p = a * p + p0;
        q = a * q + q0;
        if (p < 0 || q < 0) return -1;
        p0 = p1, q0 = q1;
    }
    if (k & 1) {
        *x1 = p, *y1 = q;
        return 0;
    } else {
        *x1 = p * p + m * q * q, *y1 = 2 * p * q;
        *x2 = p, *y2 = q;
        return 1;
    }
}

펠 방정식의 이름은 나-가 잘못해 명명한 것으로 되어 있다. 펠은 실재의 수학자이지만, 이 방정식에 대한 연구는 전혀 없었다.

펠 방정식

　　　　x² - M y² = c 　(M 는 평방수가 아닌 정정수, c = +1, -1)

에 두어, 자명한 해 x = 1, y = 0 이외의 x, y 가 가장 작은 정의 정수이다 해를, 그 펠 방정식의 최소해라고 한다. 방정식의 우변 c = -1 로 했을 때에는 해가 없는 경우가 많지만, 그러나 -1 에 대 하는 최소해가 존재하는 경우에는, 반드시 그 최소해 쪽이 +1 에 대한 최소해 보다 작다.

예를 들면, M = 2 의 경우의 방정식에 대해서는, x, y 의 페어 (1,1), (3,2), (7,5), (17,12), (41,29), (99,70), (239,169), (577,408)... 하지만 그 해이지만, -1 에 대한 최소해는 x=1, y=1 이며,+1 에 대한 최소해는 x=2, y=3 로 .

펠 방정식의 일반해는 다음의 점화식에서 요구된다.

　　　　x_n+1 = x₁x_n + M y₁y_n
　　　　y_n+1 = x₁y_n + y₁x_n

다만 x₁, y₁은 최소해를 나타낸다. 따라서 펠 방정식 (을)를 풀려면 그 최소해를 요구하면 좋다.

M 의 값에 의해, 최소해가 의외로 큰 수가 되는 일이 있다. 예를 들면, M=61, c=1 에 대한 최소해는 x=1766319049, y=226153980 로도 된다.

최소해는 다음과 같게 점화적으로 요구된다.

　　　　s₀ = 0　　t₀ = 1
　　　　a₀ = M^1/2 의 정수부
　　　　x_-1 = 1　　y_-1 = 0　　x₀ = a₀　　y₀ = 1
　　　　s_k+1 = a_kt_k - s_k　　 t_k+1 = (M - s_k+1²) / t_k
　　　　a_k+1 = (a₀ + s_k+1) / t_k+1의 정수부
　　　　x_k+1 = a_k+1 x_k + x_k-1 　　y_k+1 = a_k+1 y_k + y_k-1

t_k+1의 값이 1 이 되었을 때의, x_k, y_k의 값 하지만 펠 방정식의 최소치가 된다.