zoukankan      html  css  js  c++  java
  • 佩尔方程

    关于佩尔方程

    佩尔方程是具有(x^2-ny^2=1)形式的丢番图方程(不定方程)

    (n)为完全平方数的时候,这个方程只有平凡解((pm1,0)),对于其他情况拉格朗日证明了佩尔方程总有非平凡解。而这些解都可以由(sqrt{n})的连分数求出

    关于连分数怎么求,比如(sqrt{7}=2+frac{1}{r_1},a_0=2)

    (r_1=frac{sqrt{7}+2}{3}=1+frac{1}{r_2},a_1=1)

    (r_2=frac{sqrt{7}+1}{2}=1+frac{1}{r_2},a_2=1)

    (r_3=frac{sqrt{7}+1}{3}=1+frac{1}{r_3},a_3=1)

    (r_4=sqrt{7}+2=4+frac{1}{r_4},a_4=4)

    从这开始就是([4,1,1,1])这样的周期了

    我们取第一个周期

    我们得到了([2:1,1,1,1])

    这时我们的连分数为(frac{8}{3}),(h_3=8,k_3=3),发现(8^2-7*3^2=1)

    这就是佩尔方程的最小解了

    其实对于(frac{h_n}{k_n})(h_n,k_n)也有递推式

    (h_n=a_nh_{n-1}+h_{n-2})

    (k_n=a_nk_{n-1}+k_{n-2})

    我们由最小解可以推得通解

    此时的最小解为(x_1=8,y_1=3)

    (x_{i+1}=x_1x_i+ny_1y_i)

    (y_{i+1}=x_1y_i+y_1x_i)

    可以推出

    (egin{bmatrix}x_k\y_kend{bmatrix}=egin{bmatrix}x_1~~~~ny_1\y_1~~~~x_1end{bmatrix}^{k-1}egin{bmatrix}x_1\y_1end{bmatrix})

    可以由矩阵快速幂得到第k项通解

    例题 hdu2281-Square Number

    传送门

    解题思路:

    (x^2=frac{n(2n+1)(n+1)}{6n}=frac{(2n+1)(n+1)}{6}=frac{1}{3}(n^2+frac{3}{2}n+frac{9}{16}-frac{1}{16})=frac{1}{3}(n+frac{3}{4})^2-frac{1}{48})

    所以((4n+3)^2-48x^2=1)

    这是显然的佩尔方程形式,而这个最小解为(4n+3=7,x=1)(n=1,x=1)

    可以用最小解去递推通解,还要判一下是否(n)满足为整数

    打了个表看了下最多跑到第十六项 暴力来就好了

    #include <bits/stdc++.h>
    #include <ext/pb_ds/assoc_container.hpp>
    #include <ext/pb_ds/hash_policy.hpp>
    #include <ext/pb_ds/tree_policy.hpp>
    #include <ext/pb_ds/trie_policy.hpp>
    using namespace __gnu_pbds;
    using namespace std;
    // freopen("k.in", "r", stdin);
    // freopen("k.out", "w", stdout);
    // clock_t c1 = clock();
    // std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;
    //#pragma comment(linker, "/STACK:1024000000,1024000000")
    mt19937 rnd(time(NULL));
    #define de(a) cout << #a << " = " << a << endl
    #define rep(i, a, n) for (int i = a; i <= n; i++)
    #define per(i, a, n) for (int i = n; i >= a; i--)
    #define ls ((x) << 1)
    #define rs ((x) << 1 | 1)
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int, int> PII;
    typedef pair<double, double> PDD;
    typedef pair<char, char> PCC;
    typedef pair<ll, ll> PLL;
    typedef vector<int> VI;
    #define inf 0x3f3f3f3f
    const ll INF = 0x3f3f3f3f3f3f3f3f;
    const ll MAXN = 1e5 + 7;
    const ll MAXM = 4e5 + 7;
    const ll MOD = 1e9 + 7;
    const double eps = 1e-7;
    ll x[105], y[105];
    ll ansn[105];
    int main()
    {
        ll n = 48;
        x[1] = 7, y[1] = 1;
        ansn[1] = 1;
        for (int i = 2; i <= 20; i++)
        {
            x[i] = x[i - 1] * x[1] + n * y[1] * y[i - 1];
            if ((x[i] - 3) % 4 == 0)
                ansn[i] = (x[i] - 3) / 4;
            else
                ansn[i] = -1;
            y[i] = x[1] * y[i - 1] + y[1] * x[i - 1];
        }
        ll N;
        while (~scanf("%lld", &N) && N)
        {
            ll ans = -inf;
            for (int i = 1; i <= 16; i++)
            {
                if (ansn[i] == -1)
                    continue;
                if (N < ansn[i])
                    break;
                ans = i;
            }
            printf("%lld %lld
    ", ansn[ans], y[ans]);
        }
        return 0;
    }
    

    例题 Street Numbers

    传送门

    知道题意之后就跟上面一模一样的过程了

    代码

    #include <iostream>
    #include <stdio.h>
    #include <cmath>
    #include <algorithm>
    #include <string>
    #include <cstring>
    #include <string.h>
    #include <map>
    #include <queue>
    #include <stack>
    #include <list>
    #include <cctype>
    #include <fstream>
    #include <sstream>
    #include <iomanip>
    #include <set>
    #include <vector>
    #include <cstdlib>
    #include <time.h>
    using namespace std;
    /* freopen("k.in", "r", stdin);
    freopen("k.out", "w", stdout); */
    // clock_t c1 = clock();
    // std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;
    //#pragma comment(linker, "/STACK:1024000000,1024000000")
    #define de(a) cout << #a << " = " << a << endl
    #define rep(i, a, n) for (int i = a; i <= n; i++)
    #define per(i, a, n) for (int i = n; i >= a; i--)
    #define ls ((x) << 1)
    #define rs ((x) << 1 | 1)
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int, int> PII;
    typedef pair<double, double> PDD;
    typedef pair<ll, ll> PLL;
    typedef vector<int, int> VII;
    #define inf 0x3f3f3f3f
    const ll INF = 0x3f3f3f3f3f3f3f3f;
    const ll MAXN = 2e6 + 7;
    const ll MAXM = 4e5 + 7;
    const ll MOD = 1e9 + 7;
    const double eps = 1e-7;
    const double pi = acos(-1.0);
    ll x[20], y[20];
    int main()
    {
        vector<pair<ll, ll> > ans;
        x[1] = 3, y[1] = 1;
        ll n = 8;
        for (int i = 2; i <= 20; i++)
        {
            x[i] = x[i - 1] * x[1] + n * y[1] * y[i - 1];
            y[i] = x[1] * y[i - 1] + y[1] * x[i - 1];
            if ((x[i] - 1) % 2 == 0)
                ans.push_back({(x[i] - 1) / 2, y[i]});
            if (ans.size() == 10)
                break;
        }
        for (int i = 0; i < 10; i++)
            printf("%10lld%10lld
    ", ans[i].second, ans[i].first);
        return 0;
    }
    

    推佩尔方程的板子

    连分数版

    #include <bits/stdc++.h>
    #include <ext/pb_ds/assoc_container.hpp>
    #include <ext/pb_ds/hash_policy.hpp>
    #include <ext/pb_ds/tree_policy.hpp>
    #include <ext/pb_ds/trie_policy.hpp>
    using namespace __gnu_pbds;
    using namespace std;
    // freopen("k.in", "r", stdin);
    // freopen("k.out", "w", stdout);
    // clock_t c1 = clock();
    // std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;
    //#pragma comment(linker, "/STACK:1024000000,1024000000")
    mt19937 rnd(time(NULL));
    #define de(a) cout << #a << " = " << a << endl
    #define rep(i, a, n) for (int i = a; i <= n; i++)
    #define per(i, a, n) for (int i = n; i >= a; i--)
    #define ls ((x) << 1)
    #define rs ((x) << 1 | 1)
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int, int> PII;
    typedef pair<double, double> PDD;
    typedef pair<char, char> PCC;
    typedef pair<ll, ll> PLL;
    typedef vector<int> VI;
    #define inf 0x3f3f3f3f
    const ll INF = 0x3f3f3f3f3f3f3f3f;
    const ll MAXN = 1e5 + 7;
    const ll MAXM = 4e5 + 7;
    const ll MOD = 1e9 + 7;
    const double eps = 1e-7;
    ll a[20000];
    bool pell_minimum_solution(ll n, ll &x0, ll &y0)
    {
        ll m = (ll)sqrt((double)n);
        double sq = sqrt(n);
        int i = 0;
        if (m * m == n)
            return false; //当n是完全平方数则佩尔方程无解
        a[i++] = m;
        ll b = m, c = 1;
        double tmp;
        do
        {
            c = (n - b * b) / c;
            tmp = (sq + b) / c;
            a[i++] = (ll)(floor(tmp));
            b = a[i - 1] * c - b;
            //printf("%lld %lld %lld
    ",a[i-1],b,c);
        } while (a[i - 1] != 2 * a[0]);
        ll p = 1, q = 0;
        for (int j = i - 2; j >= 0; j--)
        {
            ll t = p;
            p = q + p * a[j];
            q = t;
            //printf("a[%d]=%lld %lld %lld
    ",j,a[j],p,q);
        }
        if ((i - 1) % 2 == 0)
        {
            x0 = p;
            y0 = q;
        }
        else
        {
            x0 = 2 * p * p + 1;
            y0 = 2 * p * q;
        }
        return true;
    }
    int main()
    {
        ll n, x, y;
        //x^2-ny^2=1
        while (~scanf("%lld", &n))
        {
            if (pell_minimum_solution(n, x, y))
            {
                printf("%lld^2-%lld*%lld^2=1	", x, n, y);
                printf("%lld-%lld=1
    ", x * x, n * y * y);
            }
        }
    }
    

    暴力

    int ax[MAXN], ay[MAXN];
    void pell(int &x, int &y, int n)
    { //暴力寻找pell方程最小解
        y = 1;
        while (true)
        {
            x = (ll)sqrt(n * y * y + 1);
            if (x * x - n * y * y == 1)
                break;
            y++;
        }
    }
    int main()
    {
        int n;
        while (scanf("%d", &n) != EOF)
        {
            int m = (int)sqrt((double)n);
            if (m * m == n)
            { //d不能为完全平方数
                cout << "No Solution" << endl;
                continue;
            }
            int x = 0, y = 0;
            pell(x, y, n); //暴力找到最小解
            printf("%d %d
    ", x, y);
        }
        return 0;
    }
    

    递推式

    for (int i = 2; i <= 20; i++)
    {
        x[i] = x[i - 1] * x[1] + n * y[1] * y[i - 1];
        y[i] = x[1] * y[i - 1] + y[1] * x[i - 1];
    }
    

    矩阵递推求第k项

    (egin{bmatrix}x_k\y_kend{bmatrix}=egin{bmatrix}x_1~~~~ny_1\y_1~~~~x_1end{bmatrix}^{k-1}egin{bmatrix}x_1\y_1end{bmatrix})

    用这个去算就好了

  • 相关阅读:
    Angular Universal 学习笔记
    SAP Spartacus 如何获得当前渲染页面的 CMS 元数据
    Angular 服务器端渲染的学习笔记(二)
    Angular 服务器端渲染的学习笔记(一)
    第三方外部 Saas提供商如何跟使用 SAP 系统的客户进行对接接口集成
    如何从 SAP Spartacus Product Detail 页面,找到其 Angular 实现 Component 的位置
    具备自动刷新功能的 SAP ABAP ALV 报表
    C++学习目录
    c--条件编译
    c--文件读写--二进制
  • 原文地址:https://www.cnblogs.com/graytido/p/13927747.html
Copyright © 2011-2022 走看看