Commits

Lars Viklund committed ae06db4

Solved #66.

  • Participants
  • Parent commits e24f44e

Comments (0)

Files changed (3)

File Windows/Windows.vcxproj

     <ClInclude Include="..\src\060.h" />
     <ClInclude Include="..\src\063.h" />
     <ClInclude Include="..\src\064.h" />
+    <ClInclude Include="..\src\066.h" />
     <ClInclude Include="..\src\068.h" />
     <ClInclude Include="..\src\069.h" />
     <ClInclude Include="..\src\070.h" />

File Windows/Windows.vcxproj.filters

     <ClInclude Include="..\src\094.h">
       <Filter>Header Files\Problems</Filter>
     </ClInclude>
+    <ClInclude Include="..\src\066.h">
+      <Filter>Header Files\Problems</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="..\src\main.cc">
+#include "problem.h"
+#include "primes.h"
+
+struct Context
+{
+	mpz_class a, p, q, P, Q;
+};
+
+template <>
+void Problem<66>()
+{
+	mpz_class bestx  = 0;
+	int bestD = 0;
+	SievedPrimes sp(1000000);
+	for (int D = 2; D <= 1000; ++D)
+	{
+		int S = D;
+		int a0 = (int)sqrt((float)S);
+		if (a0*a0 == S)
+			continue;
+		Context c0 = { a0, a0, 1, 0, 1 };
+		
+		Context ctx = {};
+		ctx.P = a0;
+		ctx.Q = D - a0*a0;
+		ctx.a = (a0 + ctx.P) / ctx.Q;
+		ctx.p = a0*ctx.a + 1;
+		ctx.q = ctx.a;
+
+		std::vector<Context> history(1, c0);
+		Context last = {};
+		while (ctx.p*ctx.p - D * ctx.q * ctx.q != 1)
+		{
+			history.push_back(ctx);
+			auto const& n1 = history[history.size()-1];
+			auto const& n2 = history[history.size()-2];
+
+			ctx.P = n1.a*n1.Q - n1.P;
+			ctx.Q = (D - ctx.P*ctx.P) / n1.Q;
+			ctx.a = (a0 + ctx.P) / ctx.Q;
+			ctx.p = ctx.a*n1.p + n2.p;
+			ctx.q = ctx.a*n1.q + n2.q;
+		}
+
+		if (ctx.p > bestx)
+		{
+			bestx = ctx.p;
+			bestD = D;
+		}
+	}
+	std::cout << "The maximum X = " << bestx << " occurs at D = " << bestD << "." << std::endl;
+}