[NBUT 1224 Happiness Hotel 佩尔方程最小正整数解]连分数法解Pell方程

题意：求方程x²-Dy²=1的最小正整数解
思路：用连分数法解佩尔方程，关键是找出√d的连分数表示的循环节。具体过程参见：http://m.blog.csdn.net/blog/wh2124335/8871535
当d为完全平方数时无解
将√d表示成连分数的形式，例如：
当d不为完全平方数时，√d为无理数，那么√d总可以表示成：
记
当n为偶数时，x0=p，y0=q；当n为奇数时，x0=2p²+1，y0=2pq
求d在1000以内佩尔方程的最小正整数解的c++打表程序（正常跑比较慢，这个题需要离线打表）：
#pragma comment(linker, "/STACK:10240000")
#include <map>
#include <set>
#include <cmath>
#include <ctime>
#include <deque>
#include <queue>
#include <stack>
#include <vector>
#include <cstdio>
#include <string>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>

using namespace std;

#define X                   first
#define Y                   second
#define pb                  push_back
#define mp                  make_pair
#define all(a)              (a).begin(), (a).end()
#define fillchar(a, x)      memset(a, x, sizeof(a))
#define copy(a, b)          memcpy(a, b, sizeof(a))

typedef long long ll;
typedef pair<int, int> pii;
typedef unsigned long long ull;

#ifndef ONLINE_JUDGE
void RI(vector<int>&a,int n){a.resize(n);for(int i=0;i<n;i++)scanf("%d",&a[i]);}
void RI(){}void RI(int&X){scanf("%d",&X);}template<typename...R>
void RI(int&f,R&...r){RI(f);RI(r...);}void RI(int*p,int*q){int d=p<q?1:-1;
while(p!=q){scanf("%d",p);p+=d;}}void print(){cout<<endl;}template<typename T>
void print(const T t){cout<<t<<endl;}template<typename F,typename...R>
void print(const F f,const R...r){cout<<f<<", ";print(r...);}template<typename T>
void print(T*p, T*q){int d=p<q?1:-1;while(p!=q){cout<<*p<<", ";p+=d;}cout<<endl;}
#endif
template<typename T>bool umax(T&a, const T&b){return b<=a?false:(a=b,true);}
template<typename T>bool umin(T&a, const T&b){return b>=a?false:(a=b,true);}

const double PI = acos(-1.0);
const int INF = 1e9 + 7;
const double EPS = 1e-12;

/* -------------------------------------------------------------------------------- */

struct BigInt {
    const static int maxI = 1e8;
    const static int Len = 8;
    typedef vector<int> vi;
    typedef long long LL;
    vi num;
    bool symbol;

    BigInt() {
        num.clear();
        symbol = 0;
    }
    BigInt(int x) {
        symbol = 0;
        if (x < 0) {
            symbol = 1;
            x = -x;
        }
        num.push_back(x % maxI);
        if (x >= maxI) num.push_back(x / maxI);
    }
    BigInt(bool s, vi x) {
        symbol = s;
        num = x;
    }
    BigInt(char s[]) {
        int len = strlen(s), x = 1, sum = 0, p = s[0] == '-';
        symbol = p;
        for (int i = len - 1; i >= p; i--) {
            sum += (s[i] - '0') * x;
            x *= 10;
            if (x == 1e8 || i == p) {
                num.push_back(sum);
                sum = 0;
                x = 1;
            }
        }
        while (num.back() == 0 && num.size() > 1) num.pop_back();
    }

    void push(int x) {
        num.push_back(x);
    }

    BigInt abs() const {
        return BigInt(false, num);
    }

    bool smaller(const vi &a, const vi &b) const {
        if (a.size() != b.size()) return a.size() < b.size();
        for (int i = a.size() - 1; i >= 0; i--) {
            if (a[i] != b[i]) return a[i] < b[i];
        }
        return 0;
    }

    bool operator < (const BigInt &p) const {
        if (symbol && !p.symbol) return true;
        if (!symbol && p.symbol) return false;
        if (symbol && p.symbol) return smaller(p.num, num);
        return smaller(num, p.num);
    }

    bool operator > (const BigInt &p) const {
        return p < *this;
    }

    bool operator == (const BigInt &p) const {
        return !(p < *this) && !(*this < p);
    }

    bool operator != (const BigInt &p) const {
        return *this < p || p < *this;
    }

    bool operator >= (const BigInt &p) const {
        return !(*this < p);
    }

    bool operator <= (const BigInt &p) const {
        return !(p < *this);
    }

    vi add(const vi &a, const vi &b) const {
        vi c;
        c.clear();
        int x = 0;
        for (int i = 0; i < a.size(); i++) {
            x += a[i];
            if (i < b.size()) x += b[i];
            c.push_back(x % maxI);
            x /= maxI;
        }
        for (int i = a.size(); i < b.size(); i++) {
            x += b[i];
            c.push_back(x % maxI);
            x /= maxI;
        }
        if (x) c.push_back(x);
        while (c.back() == 0 && c.size() > 1) c.pop_back();
        return c;
    }

    vi sub(const vi &a, const vi &b) const {
        vi c;
        c.clear();
        int x = 1;
        for (int i = 0; i < b.size(); i++) {
            x += maxI + a[i] - b[i] - 1;
            c.push_back(x % maxI);
            x /= maxI;
        }
        for (int i = b.size(); i < a.size(); i++) {
            x += maxI + a[i] - 1;
            c.push_back(x % maxI);
            x /= maxI;
        }
        while (c.back() == 0 && c.size() > 1) c.pop_back();
        return c;
    }

    vi mul(const vi &a, const vi &b) const {
        vi c;
        c.resize(a.size() + b.size());
        for (int i = 0; i < a.size(); i++) {
            for (int j = 0; j < b.size(); j++) {
                LL tmp = (LL)a[i] * b[j] + c[i + j];
                c[i + j + 1] += tmp / maxI;
                c[i + j] = tmp % maxI;
            }
        }
        while (c.back() == 0 && c.size() > 1) c.pop_back();
        return c;
    }

    vi div(const vi &a, const vi &b) const {
        vi c(a.size()), x(1, 0), y(1, 0), z(1, 0), t(1, 0);
        y.push_back(1);
        for (int i = a.size() - 1; i >= 0; i--) {
            z[0] = a[i];
            x = add(mul(x, y), z);
            if (smaller(x, b)) continue;
            int l = 1, r = maxI - 1;
            while (l < r) {
                int m = (l + r + 1) >> 1;
                t[0] = m;
                if (smaller(x, mul(b, t))) r = m - 1;
                else l = m;
            }
            c[i] = l;
            t[0] = l;
            x = sub(x, mul(b, t));
        }
        while (c.back() == 0 && c.size() > 1) c.pop_back();
        return c;
    }

    BigInt operator + (const BigInt &p) const {
        if (!symbol && !p.symbol) return BigInt(false, add(num, p.num));
        if (!symbol && p.symbol) {
            return *this >= p.abs() ?
            BigInt(false, sub(num, p.num)) : BigInt(true, sub(p.num, num));
        }
        if (symbol && !p.symbol) {
            return (*this).abs() > p ?
            BigInt(true, sub(num, p.num)) : BigInt(false, sub(p.num, num));
        }
        return BigInt(true, add(num, p.num));
    }

    BigInt operator - (const BigInt &p) const {
        return *this + BigInt(!p.symbol, p.num);
    }

    BigInt operator * (const BigInt &p) const {
        BigInt res(symbol ^ p.symbol, mul(num, p.num));
        if (res.symbol && res.num.size() == 1 && res.num[0] == 0)
            res.symbol = false;
        return res;
    }

    BigInt operator / (const BigInt &p) const {
        if (p == BigInt(0)) return p;
        BigInt res(symbol ^ p.symbol, div(num, p.num));
        if (res.symbol && res.num.size() == 1 && res.num[0] == 0)
            res.symbol = false;
        return res;
    }

    BigInt operator % (const BigInt &p) const {
        return *this - *this / p * p;
    }

    void show() const {
        if (symbol) putchar('-');
        printf("%d", num[num.size() - 1]);
        for (int i = num.size() - 2; i >= 0; i--) {
            printf("%08d", num[i]);
        }
        //putchar('
');
    }

    int TotalDigit() const {
        int x = num[num.size() - 1] / 10, t = 1;
        while (x) {
            x /= 10;
            t++;
        }
        return t + (num.size() - 1) * Len;
    }

};

template<typename T>
T gcd(T a, T b) {
    return b == 0? a : gcd(b, a % b);
}

template<typename T>
struct  Fraction {
    T a, b;
    Fraction(T a, T b) {
        T g = gcd(a, b);
        this->a = a / g;
        this->b = b / g;
        if (this->b < 0) {
            this->a = this->a * T(- 1);
            this->b = this->b * T(- 1);
        }
    }
    Fraction(T a) {
        this->a = a;
        this->b = 1;
    }
    Fraction() {}
    Fraction operator + (const Fraction &that) const {
        T x = a * that.b + b * that.a, y = b * that.b;
        return Fraction(x, y);
    }
    Fraction operator - (const Fraction &that) const {
        T x = a * that.b - b * that.a, y = b * that.b;
        return Fraction(x, y);
    }
    Fraction operator * (const Fraction &that) const {
        T x = a * that.a, y = b * that.b;
        return Fraction(x, y);
    }
    Fraction operator / (const Fraction &that) const {
        T x = a * that.b, y = b * that.a;
        return Fraction(x, y);
    }
    Fraction operator += (const Fraction &that)  {
        return *this = *this + that;
    }
    Fraction operator -= (const Fraction &that)  {
        return *this = *this - that;
    }
    Fraction operator *= (const Fraction &that)  {
        return *this = *this * that;
    }
    Fraction operator /= (const Fraction &that)  {
        return *this = *this / that;
    }
    Fraction operator ! () const {
        return Fraction(b, a);
    }
    bool operator == (const Fraction &that) const {
        return a == that.a && b == that.b;
    }
    bool operator != (const Fraction &that) const {
        return a != that.a || b != that.b;
    }
};

template<typename T>
T getInt(Fraction<T> a, T d, Fraction<T> b) {
    T Min = 0, Max;
    Fraction<T> buf = a * d + b;
    Max = buf.a / buf.b;
    while (Min < Max) {
        T Mid = (Min + Max + 1) / 2;
        buf = (b - Mid) * (b - Mid);
        buf = buf / a / a;
        if (buf.a <= buf.b * d) Min = Mid;
        else Max = Mid - 1;
    }
    return Min;
}

void work(int n) {
    int k = (int)sqrt(n + 0.5);
    if (k * k == n) {
        printf("no solution");
        return ;
    }
    Fraction<BigInt> a(1), b(0), aa, bb;
    BigInt d(n);
    vector<BigInt> R;
    BigInt t = getInt(a, d, b);
    aa = a / (a * a * d - (b - t) * (b - t));
    bb = (b - t) * BigInt(- 1) / (a * a * d - (b - t) * (b - t));
    a = aa;
    b = bb;
    do {
        R.pb(t);
        t = getInt(a, d, b);
        aa = a / (a * a * d - (b - t) * (b - t));
        bb = (b - t) * BigInt(- 1) / (a * a * d - (b - t) * (b - t));
        a = aa;
        b = bb;
    } while (t != R[0] * 2);
    Fraction<BigInt> ans(R[R.size() - 1]);
    for (int i = 1; i < R.size(); i ++) {
        ans = !ans + R[R.size() - i - 1];
    }
    BigInt x0 = ans.a, y0 = ans.b;
    if (R.size() & 1) {
        x0 = ans.a * ans.a * 2 + 1;
        y0 = ans.a * ans.b * 2;
    }
    x0.show();
}


int main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
#endif // ONLINE_JUDGE
    int n;
    puts("char ans[][100] = {"", ");
    for (int i = 1; i <= 1000; i ++) {
        printf(""");
        work(i);
        printf("", ");
        if (i % 20 == 0) puts("");
    }
    puts("
};");
    return 0;
}
相关阅读:
JavaScript 数组去重
 Javascript数组 reduce() 方法
 Vue事件总线(EventBus)
前端多媒体-音频
 前端多媒体-视频
 VUE3.0 总结
 el-select地区区号选择
 vue中点击获取相应元素
 Markdown 语法
 Codeforces Round #295 (Div. 2) B. Two Buttons 520B
原文地址：https://www.cnblogs.com/jklongint/p/4743778.html