Project Euler_Problem 140_Modified Fibonacci Golden Nuggets_生成函数+广义佩尔方程

原题目：

题目大意：

解题思路：

这时候就可以去暴力搜索答案了

代码：

import math

maxn = 211345365
x= -1
p=1;
k=2

while(p<=30):
    u=5*k*k+14*k+1;
    j=int(math.sqrt(u));
    if(u==j*j):
        print(p,k);
        p=p+1;
    k=k+1

但是很遗憾，这个题，数据量要求有点大了，我Python跑了半个小时也只能跑到第21个

考虑进一步优化，转化成求二次方程解的问题：

此时该程序已经变得非常强大。

代码:

#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>  //cout输出精度控制需要，fixed、setpercison（）
#include <time.h> 
#include <map> 
using namespace std;
#define ll long long



ll A[2000000], B[2000000], C[2000000],D[30000000],D1[30000000];
ll E[1000][1000];
ll R[100][50][50];
bool bA[10000000];

map <ll, bool> Map4;


ll ans1 = 0, ans2 = 0, flag;


struct Math{   //自用函数模板

    const double Pi = acos(-1);
    const double E = 2.718281828459;

    struct STL{
        
    }Stl;

    struct Num_Theory {

        /*struct Pythagorean_triangle {
            ll a, b, c;
        }Ptg_r[10000000];
        ll Ptg_r_cnt=0;
        void get_Ptg_r(ll tar) {
            ll i, j, k,p,q;
            ll x, y, z;
            for (i = 1; i * i < tar; i++) {
                if (i & 1)j = 2;
                else j = 1;
                for (; j < i; j+=2) {
                    x = 2 * i * j;
                    y = i * i - j * j;
                    z = i * i + j * j;
                    if (x + y + z > tar)break;
                    if (x > y) {
                        k = y;
                        y = x;
                        x = k;
                    }

                    if (y - x != 1)continue;

                    if (gcd(i, j) == 1) {
                       // Ptg_r[++Ptg_r_cnt].a = x; Ptg_r[Ptg_r_cnt].b = y; Ptg_r[Ptg_r_cnt].c = z;
                        Ptg_r[0].a = x; Ptg_r[0].b = y; Ptg_r[0].c = z;
                        if (M.NT.Ptg_r[0].c % (M.NT.Ptg_r[0].b - M.NT.Ptg_r[0].a) == 0) {
                            p = M.NT.Ptg_r[0].a + M.NT.Ptg_r[0].b + M.NT.Ptg_r[0].c;
                            for (q = 1; q * p < tar; q++) {
                                ans1++;
                            }


                        }


                    }
                }
            }
        }*/


        ll Rad[1000000];
        void get_rad(ll x) {
            ll i, j,k;

            for (i = 1; i <= x; i++) {
                Rad[i] = 1;
                k = i;
                for (j = 1; prime[j] <= i; j++) {
                    
                    if (k % prime[j] == 0) {
                        Rad[i] *= prime[j];

                        while (k % prime[j] == 0) {
                            k = k / prime[j];
                        }
                    }
                }
            }

        }

        ll prime_mul(ll x, ll y, ll mod) {
            ll ans = 0;

            while (y) {
                if (y & 1) {
                    ans = (ans + x) % mod;
                }
                x = (x + x) % mod;
                y >>= 1;
            }
            return ans;

        }


        ll prime_pow(ll x,ll y,ll mod) {
            ll ans = 1,u;
            while (y) {
                if (y & 1) {
                   ans = (prime_mul(ans, x, mod)) % mod; 
                }
                   x = (prime_mul(x, x, mod)) % mod;
                y >>= 1;
            }
            return ans % mod;
        }


        ll prime[6000000],phi[50000000], pcnt;
        bool prime_vis[100000000];
        bool is_prime_MR(ll x,ll r) {
            if (x == 2)return 1;
            if (!(x%2))return 0;
            ll i, j,k,t,X,Y;

            k = x-1;
            t = 0;
            while (!(k&1)) {
                k >>= 1;
                t++;
            }
            srand(time(NULL));
            for (i = 1; i <= r; i++) {

                ll a = rand() % (x - 1) + 1;
                X = prime_pow(a, k, x);

                for (j = 0; j < t; j++) {
                    Y = prime_mul(X, X, x);

                    if (Y == 1 && X != 1 && X != x - 1) {
                        return 0;
                    }
                    X = Y;
                }

                if (X != 1) {
                    return 0;
                }

            }
            return 1;

        }

        bool is_prime_Bf(ll x) {
            ll i, j,k,cnt=0;

            if (x == 2)return 1;
            if (x % 2 == 0)return 0;
            for (i = 2; i <= sqrt(x); i++) {
                if (x % i == 0)return 0;
            }
            return 1;

        }

        void get_prime_Euler(ll x) {
            pcnt = 0; memset(prime_vis, 0, sizeof(prime_vis));
            ll i, j;

            for (i = 2; i <= x; i++) {
                if (prime_vis[i] == 0) {
                    prime[++pcnt]=i;
                }
                for (j = 1; j <= pcnt; j++) {
                    if (i * prime[j] > x)break;
                    prime_vis[i * prime[j]] = 1;
                    if (i % prime[j] == 0)break;

                }
            }

        }

        void get_prime_phi(ll x) {
            ll i, j; pcnt = 0; memset(prime_vis, 0, sizeof(prime_vis));
            for (i = 2; i <= x; i++) {
                if (prime_vis[i] == 0) {
                    prime[++pcnt] = i;
                    phi[i] = i - 1;
                }
                for (j = 1; j <= pcnt; j++) {
                    if (i * prime[j] > x)break;
                    prime_vis[prime[j] * i] = 1;

                    if (i % prime[j] == 0) {
                        phi[i * prime[j]] = prime[j] * phi[i];
                        break;
                    }
                    phi[i * prime[j]] = (prime[j] - 1) * phi[i];


                }



            }
        }

        ll gcd(ll a, ll b) {
            if (a < b) {
                ll c = a;
                a = b; b = c;
            }
            if (a % b == 0)return b;
            else return gcd(b, a % b);

        }


        ll Divisible_Checker(ll x, ll y) {
            x = x % y;
            if(x==0)return 1;
            ll i, j, z;
            z = gcd(x, y);
            while (z != 1) {
                x = x / z;
                y = y / z;
                z = gcd(x, y);
            }

            while (y % 2 == 0) {
                y = y / 2;
            }
            while (y % 5 == 0) {
                y = y / 5;
            }

            if (y != 1) {
                return 0;
            }

            return 1;
        }


        ll Pell_Base_solve[100][2]; 
        ll Pell_solve[1000],PS_Len;
        void Gener_Pell_search(ll P,ll Q,ll D,ll len,ll tar) {

            ll i, j,k, u, v, x, y,flag;
            PS_Len = 0;
            while (PS_Len <tar) {
                k = 1;
                flag = 1;
                for (i = 1; i <= len; i++) {
                    x = Pell_Base_solve[i][0];
                    y = Pell_Base_solve[i][1];
                    if ((x - 7) % 5 != 0) {
                         flag = 0;
                         break;
                    }

                }

                if (flag == 1) {
                    u = Pell_Base_solve[1][0];
                    v = Pell_Base_solve[1][1];
                    for (i = 1; i <= len; i++) {
                        x = Pell_Base_solve[i][0];
                        y = Pell_Base_solve[i][1];

                        if (x < u) {
                            u = x; k = i;
                        }

                    }

                }
                if (flag == 1) {
                    i = k;
                    ll t= (Pell_Base_solve[i][0] - 7) / 5;
                    if (t != Pell_solve[PS_Len]) {
                        Pell_solve[++PS_Len] = t;
                    }
                }

                x = Pell_Base_solve[i][0];
                y = Pell_Base_solve[i][1];

                u = P * x + Q * D * y;
                v = P * y + Q * x;

                x = u, y = v;
                Pell_Base_solve[i][0] = x; Pell_Base_solve[i][1] = y;


            }


        }

    }NT;

    struct Num_Analysis {
        double x[100], y[100];

        double Lag_in(ll X, ll K) {
            int i, j;
            double ans = 0;
            double li;

            for (i = 1; i <= K; i++) {
                li = 1;
                for (j = 1; j <= K; j++) {
                    if (i == j)continue;
                    li = li * ((X - x[j]) / (x[i] - x[j]));
                }
                ans = ans + (li) * (y[i]);
            }

            return ans;
        }


    }NA;
    
    struct ComPA {
        struct Complex {
            ll x, y;
            double a, b;
        };

        Complex complex_Add(Complex A, Complex B) {
            Complex ans;
            ans.x = A.x + B.x;
            ans.y = A.y + B.y;
            ans.a=A.a + B.a;
            ans.b = A.b + B.b;
            return ans;
        }

        Complex complex_Sud(Complex A, Complex B) {
            Complex ans;
            ans.x = A.x - B.x;
            ans.y = A.y - B.y;
            ans.a = A.a - B.a;
            ans.b = A.b - B.b;
            return ans;
        }


        Complex complex_Mul(Complex A, Complex B) {
            Complex ans;
            ans.x = A.x*B.x-A.y*B.y;
            ans.y = A.x*B.y + A.y*B.x;

            ans.a = A.a * B.a - A.b * B.b;
            ans.b = A.a * B.b + A.b * B.a;
            return ans;
        }


    }CPA;

    struct Lin_Algebra {

        struct Complex {
            double a, b;
        }Cpl[100];
        Complex complex_Add(Complex A, Complex B) {
            Complex ans;
            ans.a = A.a + B.a;
            ans.b = A.b + B.b;
            return ans;
        }
        Complex complex_Sud(Complex A, Complex B) {
            Complex ans;
            ans.a = A.a - B.a;
            ans.b = A.b - B.b;
            return ans;
        }
        Complex complex_Mul(Complex A, Complex B) {
            Complex ans;
            ans.a = A.a * B.a - A.b * B.b;
            ans.b = A.a * B.b + A.b * B.a;
            return ans;
        }



        Complex FFT1[1000];
        void FFT_BF(Complex A[],ll n) {
            if (n == 1) {
                return;
            }
            ll i, j;
            Complex a[100], b[100];

            for (i = 0,j=0; i < n; i = i + 2) {
                a[j++] = A[i];
            }
            for (i = 1, j = 0; i < n; i += 2) {
                b[j++] = A[i];
            }
            
            FFT_BF(a, n / 2); FFT_BF(b, n / 2);

            Complex wn;

            for (i = 0; i < (n>>1); i++) {
                wn.a = cos((2.0 * i * M.Pi) / (double)n);
                wn.b = sin((2.0 * i * M.Pi) / (double)n);

                A[i].a = a[i].a+ complex_Mul(wn,b[i]).a;
                A[i].b = a[i].b + complex_Mul(wn, b[i]).b;

                A[i+n/2].a = a[i].a - complex_Mul(wn, b[i]).a;
                A[i+n/2].b = a[i].b - complex_Mul(wn, b[i]).b;

            }

        }



    }LA;

}M;


void serch1(ll x,ll y) {

}

void solve() {
    ll i, j,k,x,y,z,p,q,u,v;
    ll N = 30,NN=6000000;
    k = 2;
    u = 0;
    v = 0;
    z = 1;
    double a,b,c;
    x = 7, y = 1;
    
    M.NT.Pell_Base_solve[1][0] = 7;
    M.NT.Pell_Base_solve[1][1] = 1;
    M.NT.Pell_Base_solve[2][0] = 8;
    M.NT.Pell_Base_solve[2][1] = 2;
    M.NT.Pell_Base_solve[3][0] = 13;
    M.NT.Pell_Base_solve[3][1] = 5;
    M.NT.Pell_Base_solve[4][0] = 17;
    M.NT.Pell_Base_solve[4][1] = 7;
    M.NT.Pell_Base_solve[5][0] = 32;
    M.NT.Pell_Base_solve[5][1] = 14;
    M.NT.Pell_Base_solve[6][0] = 43;
    M.NT.Pell_Base_solve[6][1] = 19;




    M.NT.Gener_Pell_search(9, 4, 5, 6, N);

    for (i = 1; i <= N; i++) {
        ans1 = ans1 + M.NT.Pell_solve[i];
        printf("%lld\n", M.NT.Pell_solve[i]);
    }

    printf("%lld\n",ans1);
}

int main()
{

    long test;
    clock_t start, finish;  // clock_t为时钟计时单元数
    start = clock();        // clock()返回此时CPU时钟计时单元数
    solve();
    finish = clock();
    cout << "代码运行花费时间为：" << fixed << setprecision(8)//控制时间输出精度8位，右侧自动补0。
        << double(finish - start) / CLOCKS_PER_SEC << "s" << endl;  //时间计算过程
    
}