From 341ea46a46f8d672061db1ad9d22a92acdc49261 Mon Sep 17 00:00:00 2001
From: Jesse McDonald <nybble41@gmail.com>
Date: Fri, 18 Sep 2015 23:26:45 -0500
Subject: [PATCH] WIP

---
 Problem066.hs | 47 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 Problem066.hs

diff --git a/Problem066.hs b/Problem066.hs
new file mode 100644
index 0000000..93cb5ea
--- /dev/null
+++ b/Problem066.hs
@@ -0,0 +1,47 @@
+-- Consider quadratic Diophantine equations of the form:
+--
+--   x^2 – D*y^2 = 1
+--
+-- It can be assumed that there are no solutions in positive integers when D is square.
+--
+-- Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.
+import Euler
+import Data.Maybe
+import Debug.Trace
+
+intSqrt n = go 0 n where
+  go a b |   m^2   >  n = go a (m-1)
+         | (m+1)^2 <= n = go (m+1) b
+         | otherwise    = m
+    where m = (a + b) `div` 2
+
+-- Chakravala method:
+--   (a^2-n*b^2 = k)
+--   1. pick any solution (a,b,k) such that gcd(a,b)=1
+--   2. if k == 1, return solution (a,b)
+--   3. choose a positive integer m such that (a + b*m)/k is an integer minimizing m^2-n
+--   4. compose with (m,1,m^2-n) to get a'=(a*m+n*b)/k, b'=(a+b*m)/k, k'=(m^2-n)/k
+--   5. repeat from step 2 with (a', b', k')
+chakravala n = let a0 = (intSqrt n+1)^2 in go (a0,1,a0^2-n)
+  where
+    go (a,b,0) = Nothing
+    go (a,b,1) = Just (a, b)
+    go (a,b,k) = traceShow (a,b,k) $ go ((a*m+n*b) `div` abs k, (a+b*m) `div` abs k, (m^2-n) `div` k)
+      where
+        fI :: (Integral a, Num b) => a -> b
+        fI = fromIntegral
+        -- TODO: Find a solution for ??? that yields integer m_i
+        m_i i = ((i*(lcm b k `div` k)-???)*k-a) `div` b
+        i0 = floor $ ((fI a + fI b * sqrt (fI n))/(fI k) + ???) / fI (lcm b k `div` k)
+        m = minimumBy (compare `on` (\m -> abs (m^2-n))) $ traceShowId $ 
+              [ m' | m' <- map (m_i . (i0+)) [-3..3]
+                   , (a*m'+n*b)       /= 0
+                   , (a+m'*b)         /= 0
+              ]
+
+-- a+m*b `mod` k == 0
+-- a+m*b-i*k == 0
+-- m == (i*k-a)/b == i*k/b - a/b
+-- i*k `mod` b == a `mod` b
+
+main = print $ maximumBy (compare `on` snd) [ traceShowId (d, x) | d <- [2..1000], Just (x, y) <- [chakravala d] ]