ZOJ 3759 - ZOJ Monthly, March 2014 pell方程求解

原题：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=5204

配图好萌啊....

题意：给出p, A

当p == 3时，求关于n, y的方程 n(n + 1)/2 = A * y^2 的最小正整数解

当p == 5时，求关于n, y的方程 n(3n -1)/2 = A * y^2 的最小正整数解

当p == 6时，求关于n, y的方程 n(2n -1)    = A * y^2 的最小正整数解

以p == 3为例

        n (n + 1) / 2 = A * y^2

=>   n^2 + n = 2A * y^2

=>   n^2 + n - 2A * y^2 = 0

把上式看成关于n的方程，则有 dlt = b^2 - 4ac = 1 + 8A * y^2

因为 n = (-b + sqrt(dlt)) / 2a 是整数，则sqrt(dlt) 必然为整数，不妨设 sqrt(dlt) = x 

注：因为本题-b <= 1，-b - sqrt(dlt)不会为解，就不讨论了

则有 1 + 8A * y^2 = x^2

即 x^2 - 8A * y^2 = 1

这样就化成pell方程 x^2 - D * y^2 = 1的形式了

因为原题要最小解，所以这里dlt也要求最小解，即pell方程要的也是最小解

所以我们求得pell方程的一组最小解(x, y)后，这个y就是原式中的y，另有 n = (x - b) / 2a

那么现在原式存在解的条件就是 pell方程有解，且方程存在一个解满足(x - b) % 2a 要等于0

我们设pell方程的最小解为(X1, Y1), 所有解为(Xi, Yi)

这里我们在证明一个性质：在本题的3个pell方程中，若(X1 - b) % 2a != 0，则(Xi - b) % 2a也均不等于0

首先我们列出本题的三个关于n的方程以及相应的pell方程

p == 3 时：  n^2 + n - 2A * y^2 = 0, 有 a = 1, b =  1, pell方程 x^2 - 8A * y^2 = 1, 即 D = 8A

p == 5 时：3n^2 -  n - 2A * y^2 = 0, 有 a = 3, b = -1, pell方程 x^2 - 24A * y^2 = 1, 即 D = 24A

p == 6 时：2n^2 -  n -   A * y^2 = 0, 有 a = 2, b = -1, pell方程 x^2 - 8A * y^2 = 1, 即 D = 8A

这三种情况都有同一个性质，即 D % 2a == 0，这是一个重要性质

因为对于pell方程的全部解，我们有递推式

Xi+1 = X1 * Xi + D * Y1 * Yi

Yi+1 = X1 * Yi + Xi * Y1

由此我们可以推出 Xi % 2a = (X1 % 2a) ^ i

对于p == 3, 若(X1 - b) % 2a = (X1 -  1) % 2 != 0  => X1 % 2 != 1  => X1 % 2 == 0  =>  Xi % 2 == 0

对于p == 5, 若(X1 - b) % 2a = (X1 + 1) % 6 != 0  => X1 % 6 != 5

那么X1 % 6 只能为 0 ~ 4, 分别算一下就知道这些数的任意次方都不可能模6得5

对于p == 6的情况分析起来和p == 5是一样的

于是现在原式有解的条件就是pell方程存在解，且最小解(x, y)满足(x - b) % 2a == 0

Pell方程有解的条件是D不为完全平方数，在Pell方程有解的情况下，求Pell方程最小解的C++代码如下：

//
bool pell( int D, int& x, int& y ) {
    int sqrtD = sqrt(D + 0.0);
    if( sqrtD * sqrtD == D ) return false;
    int c = sqrtD, q = D - c * c, a = (c + sqrtD) / q;
    int step = 0;
    int X[] = { 1, sqrtD };
    int Y[] = { 0, 1 };
    while( true ) {
        X[step] = a * X[step^1] + X[step];
        Y[step] = a * Y[step^1] + Y[step];
        c = a * q - c;
        q = (D - c * c) / q;
        a = (c + sqrtD) / q;
        step ^= 1;
        if( c == sqrtD && q == 1 && step ) {
            x = X[0], y = Y[0];
            return true;
        }
    }
}

能够求出最小解，剩下的就好办了，判断(x - b) % 2a是否为0即可，坑爹的是需要高精

Java代码：

//
import java.math.BigInteger;
import java.util.Scanner;

public class Main {
    static class Pell {
        int D;
        BigInteger x, y;
        boolean ok;
        void solve() {
            int sqrtD = (int)Math.sqrt((double)D);
            if( sqrtD * sqrtD == D ) {
                ok = false;
                return ;
            }
            else ok = true;
            BigInteger N = BigInteger.valueOf(D);
            BigInteger SqrtD = BigInteger.valueOf(sqrtD);
            BigInteger c = SqrtD;
            BigInteger q = N.subtract(c.pow(2));
            BigInteger a = c.add(SqrtD).divide(q);
            int step = 0;
            BigInteger X[] = {BigInteger.ONE, SqrtD};
            BigInteger Y[] = {BigInteger.ZERO, BigInteger.ONE};
            while( true ) {
                X[step] = a.multiply(X[step^1]).add(X[step]);
                Y[step] = a.multiply(Y[step^1]).add(Y[step]);
                c = a.multiply(q).subtract(c);
                q = N.subtract(c.pow(2)).divide(q);
                a = SqrtD.add(c).divide(q);
                step ^= 1;
                if( c.equals(SqrtD) && q.equals(BigInteger.ONE) && step != 0 )
                    break;
            }
            x = X[0];
            y = Y[0];
        }
    }
    public static void main(String args[]) {
        Scanner cin = new Scanner(System.in);
        int p, A, D, a = 0, b = 0;
        while( (p = cin.nextInt()) != 0 ) {
            A = cin.nextInt();
            Pell pell = new Pell();
            if( p == 3 ) {
                pell.D = 8 * A;
                a = 1;
                b = 1;
            }
            if( p == 5 ) {
                pell.D = 24 * A;
                a = 3;
                b = -1;
            }
            if( p == 6 ) {
                pell.D = 8 * A;
                a = 2;
                b = -1;
            }
            pell.solve();
            if( pell.ok ) {
                BigInteger X = pell.x.subtract(BigInteger.valueOf(b));
                if( X.mod(BigInteger.valueOf(2 * a)).equals(BigInteger.ZERO) ) {
                    X = X.divide(BigInteger.valueOf(2 * a));
                    System.out.println( X + " " + pell.y );
                }
                else System.out.println( "Impossible!" );
            }
            else System.out.println( "Impossible!" );
        }
    }
}