-- 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] ]