Project Euler Problem 66 (C++和Python)

Problem 66 : Diophantine equation

Consider quadratic Diophantine equations of the form:

x² – Dy² = 1

For example, when D=13, the minimal solution in x is 649² – 13×180² = 1.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:

3² – 2×2² = 1
2² – 3×1² = 1
9² – 5×4² = 1
5² – 6×2² = 1
8² – 7×3² = 1

Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.

Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.

C++ 代码

#include <iostream>
#include <vector>
#include <iterator>
#include <cmath>
#include <ctime>
#include <cassert>

using namespace std;

#define UNIT_TEST
// #define PE0066_DEBUG

class PE0066
{
private:
    vector<int> getContinuedFraction(int N);
    vector<int> getFraction_x(int D, vector<int>& b);
    void displayFraction(int D,vector<int>& n_digit_vec,vector<int>& d_digit_vec);
    void displayFraction_x(int D, vector<int>& x_digit_vec);
    bool compareVectors(vector<int>& a, vector<int>& b);

public:
    int get_D_to_obtained_largest_x(int max_D);
};

vector<int> PE0066::getContinuedFraction(int N)
{
    vector<int> fractions_vec;
    int period = 0;
    int m = 0, d = 1;
    double root = sqrt((double)N);
    int a0 = int(root);
    int a  = a0;

    fractions_vec.push_back(a0);

    while (a != 2*a0)
    {
        m = a*d - m;
        d = (N - m*m) / d;
        a = int((root + m) / d);

        period++;
        fractions_vec.push_back(a);
    }

#ifdef PE0066_DEBUG
    cout << "√" << N << " = [" << fractions_vec[0] << "; (" ;
    copy(fractions_vec.begin()+1, fractions_vec.end()-1, 
         ostream_iterator<int>(cout, ","));
    int size=fractions_vec.size();
    cout << fractions_vec[size-1] << ")], period = " << period << endl;
#endif

    return fractions_vec;
}
    
void PE0066::displayFraction(int D, vector<int>& n_digit_vec, vector<int>& d_digit_vec)
{
    cout << "D = " << D << ":  x =  "  ;
    copy(n_digit_vec.rbegin(), n_digit_vec.rend(), ostream_iterator<int>(cout));
    cout << "(" << n_digit_vec.size() << "), y = ";
    copy(d_digit_vec.rbegin(), d_digit_vec.rend(), ostream_iterator<int>(cout));
    cout << "(" << d_digit_vec.size() << ")" << endl;
}

void PE0066::displayFraction_x(int D, vector<int>& x_digit_vec)
{
    cout << "when D = " << D << ", largest x = "  ;
    copy(x_digit_vec.begin(), x_digit_vec.end(), ostream_iterator<int>(cout));
    cout << endl;
}

vector<int> PE0066::getFraction_x(int D, vector<int>& b)  
{                                         // vector b is continued fraction of D
    vector<int> tmp = b;
    int period = b.size() - 1;

    if (period % 2 == 0) 
    {
        b.pop_back();
    }
    else
    {
        for(int i=1; i<period; i++)
        {
            b.push_back(tmp[i]);
        }
    }

    vector<int> n_digit_vec;  // the digits of the numerator
    vector<int> d_digit_vec;  // the digits of the denominator 

    for(unsigned k=1; k<b.size(); k++)
    {
        n_digit_vec.clear();
        d_digit_vec.clear();

        n_digit_vec.push_back(1);
        d_digit_vec.push_back(0);
        
        for (int i=k; i>=0; i--)
        {
             swap(n_digit_vec,d_digit_vec);

            if (n_digit_vec.size() < d_digit_vec.size())
            {
                n_digit_vec.push_back(0); 
            }

            for(unsigned int m=0; m<d_digit_vec.size(); m++)
            {
                n_digit_vec[m] += b[i]*d_digit_vec[m];

                if (n_digit_vec[m] / 10 > 0)
                {
                    if (m+1<n_digit_vec.size())
                    {
                        n_digit_vec[m+1] += n_digit_vec[m]/10;
                    }
                    else
                    {
                        n_digit_vec.push_back(n_digit_vec[m]/10);
                    }
                    n_digit_vec[m] %= 10;
                }
            }
        }
    }

#ifdef PE0066_DEBUG
    displayFraction(D, n_digit_vec, d_digit_vec);
#endif

    vector<int> x_digit_vec;
    x_digit_vec.assign(n_digit_vec.rbegin(), n_digit_vec.rend());

    return x_digit_vec;
}

bool PE0066::compareVectors(vector<int>& a, vector<int>& b)
{
    if (a.size() > b.size())
    {
        return true;
    }
    else if (a.size() == b.size())
    {
        for(unsigned int i=0; i<a.size();i++)
        {
            if (a[i] == b[i]) 
            {
                continue;
            }
            else if (a[i] >= b[i])
            {
                return true;
            }
            else
            {
                return false;
            }
        }
    } 
    return false;
}

int PE0066::get_D_to_obtained_largest_x(int max_D)
{
    vector<int> largest_x_vec;
    int D_value = 0;

    for (int D=2; D<=max_D; D++)
    {
        int root = (int)sqrt(double(D));;
        if (root*root == D)
        {
            continue;
        }
        else
        {
            vector<int> fractions_vec = getContinuedFraction(D);
            vector<int> x_digits_vec  = getFraction_x(D, fractions_vec);

            if (true == compareVectors(x_digits_vec, largest_x_vec))
            {
                largest_x_vec = x_digits_vec;
                D_value       = D;
#ifdef UNIT_TEST
                displayFraction_x(D, largest_x_vec);
#endif
            }

        }
    }

    return D_value;
}

int main()
{
    PE0066 pe0066;

    assert ( 5 == pe0066.get_D_to_obtained_largest_x(7));
    assert (13 == pe0066.get_D_to_obtained_largest_x(20));

    int D = pe0066.get_D_to_obtained_largest_x(1000);

    cout << "for D <= 1000, the largest x is obtained when D = " << D << endl;

    return 0;
}

Python代码

import math

def getContinuedFraction(N):
    period, m, d = 0, 0, 1
    root = int(math.sqrt(N))
    a0, a = root, root
    fractions_list = [ a0 ]

    while a != 2*a0:
        m = a*d - m
        d = (N - m*m) / d
        a = int((root + m) / d)

        period += 1
        fractions_list += [ a ]

    return fractions_list

def getFraction_x(D, b): # vector b is continued fraction of D
    tmp, period = b, len(b)-1

    if period % 2 == 0:
        b.pop()
    else:
        for i in range(1, period):
            b.append(tmp[i])

    for k in range(1, len(b)):
        n_digit_list = [ 1 ]    # n_digit_list: the digits of the numerator
        d_digit_list = [ 0 ]    # d_digit_list: the digits of the denominator
            
        for i in range(k, -1, -1):
            n_digit_list, d_digit_list = d_digit_list, n_digit_list 

            if len(n_digit_list) < len(d_digit_list):
                n_digit_list.append(0) 

            for m in range(0, len(d_digit_list)):
                n_digit_list[m] += b[i]*d_digit_list[m];

                if n_digit_list[m] // 10 > 0:
                    if m+1 < len(n_digit_list):
                        n_digit_list[m+1] += n_digit_list[m]//10
                    else:
                        n_digit_list.append(n_digit_list[m]//10)
                    n_digit_list[m] %= 10

    x_digit_list = n_digit_list[::-1]

    return x_digit_list

def getValueOfDigitsList(x_digits_list):
    """
    convert digits list to integer
    """
    digits_str_list = [ str(i) for i in x_digits_list ]
    x_str = ''.join(digits_str_list)
    if x_str == '':
        return 0
    else:
        return int(x_str)

def displayFraction_x(D, largest_x_digits_list):
    """
    print debug information when Enable_Debug is True
    """
    print("when D = ",D,", largest x = ",getValueOfDigitsList(largest_x_digits_list))

def compareTwoDigitsLists(list_a, list_b):
    """
    list_a, list_b are two x_digits_list
    compare list_a and list_b, True if list_a > list_b, otherwise False
    """
    return getValueOfDigitsList(list_a) > getValueOfDigitsList(list_b)

def get_D_to_obtained_largest_x(max_D):
    D_value, largest_x_digits_list = 0, []

    for D in range(2, max_D+1):
        if math.sqrt(D) % 1 == 0.0:
            continue
        else:
            fractions_list = getContinuedFraction(D)
            x_digits_list  = getFraction_x(D, fractions_list)

            if True == compareTwoDigitsLists(x_digits_list, largest_x_digits_list):
                largest_x_digits_list, D_value = x_digits_list, D
                if True == Enable_Debug:
                    displayFraction_x(D, largest_x_digits_list)
    return D_value

def main():
    global Enable_Debug
    Enable_Debug = False

    assert  5 == get_D_to_obtained_largest_x(7)
    assert 13 == get_D_to_obtained_largest_x(20)

    Enable_Debug = True

    D = get_D_to_obtained_largest_x(1000)
    print("for D <= 1000, the largest x is obtained when D = %d" % D)

if  __name__ == '__main__':
    main()