AOJ 2116 Subdividing a Land (JAG 冬合宿 2007 E)

Pell 方程式！！！

• 問題へのリンク
• editorial (day2-E)

問題概要 (意訳)

正の整数 が与えられる。 を正の整数として、

の値が最小となるような を 1 つ求めよ。

制約

• 

考えたこと

が平方数のときは として、

• 

が答えとなる (最適値は 0)。 が平方数でないときは、Pell 方程式と呼ばれ、

を満たす整数 が存在する。これが最適解となる。具体的には、連分数展開によって求められる (以下の kirika さんの整数論テクニック集参照)。

github.com

コード

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

// Pell's Equation
pair<long long, long long> Pell(long long D, int c = 1) {
    if (D == 1) return make_pair(-1, -1);
    long long a = D, b = 1;
    while (a > b) {
        a /= 2;
        b *= 2;
    }
    a = b + 1;
    while (a > b) {
        a = b; 
        b = (b + D/b)/2;
    }
    if (a*a == D) return make_pair(-1, -1);
    
    long long a0 = a;
    bool parity = false;
    b = 1;
    long long x2 = 1, x = a, y2 = 0, y = 1, q;
    while (true) {
        b = (D - a*a) / b;
        q = (a0 + a) / b;
        a = q * b - a;
        parity = !parity;
        if (b == 1) break;
        long long tx = x, tx2 = x2, ty = y, ty2 = y2;
        x2 = tx, x = tx * q + tx2; 
        y2 = ty, y = ty * q + ty2;
    }
    long long x0 = x, y0 = y;
    if (!parity) {
        if (c == 1) return make_pair(x, y);
        else return make_pair(-1, -1);
    }
    else if (c == -1) return make_pair(x, y);
    
    long long tx = x, ty = y;
    x = x0 * tx + D * y0 * ty, y = tx * y0 + x0 * ty;
    return make_pair(x, y);
}

pair<long long, long long> PellMul(pair<long long, long long> p, pair<long long, long long> q, long long D) {
    long long f = p.first * q.first + D * p.second * q.second;
    long long s = p.first * q.second + p.second * q.first;
    return make_pair(f, s);
}


void AOJ2116() {
    long long N;
    int id = 0;
    while (cin >> N, N) {
        N *= 2;
        long long sq = (long long)(sqrt(N) + 0.5);
        if (sq * sq == N) cout << "Case " << ++id << ": " << sq << " " << 1 << endl;
        else {
            auto res = Pell(N);
            cout << "Case " << ++id << ": " << res.first << " " << res.second << endl;
        }
    }
}

int main() {
    AOJ2116();
}