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の値 がペル方程式の最小値となる。