Euler Problem 66:Investigate the Diophantine equation x² − Dy² = 1.

Varijanta 1. Korištenjem Mathematica

Rješenje se može vrlo lako dobiti korištenjem Wolframove Mathematike:

maxX = 0;
Dmin = 0;
For[d = 2, d <= 1000, d++,
If[Mod[Length[Divisors[d]], 2] == 1, Continue[]];
solX = x /.FindInstance[x^2 - d*y^2 == 1 && x > 0 && y > 0, {x, y},Integers][[1]];
If[solX > maxX, maxX = solX; Dmin = d]]; Dmin

Varijanta 2: C# Implementacija sa BigNumber strukturom:

using System;
using System.Collections.Generic;
using System.Numerics;
namespace EulerProblem66
{
    class Program
    {
        static void Main(string[] args)
        {
            BigInteger x = 0;
            BigInteger maxx = 0;
            int maxd = 0;
            for (int d = 2; d <= 1000; d++)
            {
                if (d % Math.Sqrt(d) == 0)
                    continue;
                List<long> c = ContinedFraction(d);
                x = PellEq(c);
                if (x > maxx)
                {
                    maxx = x;
                    maxd = d;
                }
            }
            Console.WriteLine(maxd);
            ///Kraj
            Console.WriteLine("Pres any key to continue..");
            Console.Read();
        }
        static long GCD(long a, long b)
        {
            long r = a % b;
            if (r == 0)
                return b;
            else
                return GCD(b, r);
        }
        static List<long> ContinedFraction(long d)
        {
            List<long> e = new List<long>();
            long d2 = (long)Math.Sqrt(d);
            long a = d2;
            long p = 0;
            long q = 1;
            do
            {
                e.Add(a);
                p = a * q - p;
                q = (d - p * p) / q;
                a = (p + d2) / q;
            } while (q != 1);
 
            e.Add(a);
            return e;
        }
        static BigInteger PellEq(List<long> c)
        {
            int l = c.Count - 1;
            int per = l % 2 == 0 ? l - 1 : 2 * l - 1;
 
            BigInteger a = c[0];
            BigInteger a1 = 1;
            BigInteger b = 1;
            BigInteger b1 = 0;
            BigInteger t;
 
            a = c[0];
            b = c[0] * b + b1;
            for (int i = 1; i <= per; i++)
            {
                t = a;
                a = c[(i - 1) % l + 1] * a + a1;
                a1 = t;
                t = b;
                b = c[(i - 1) % l + 1] * b + b1;
                b1 = t;
            }
            return a;
        }
    }
}