48 lines
1.8 KiB
Haskell
48 lines
1.8 KiB
Haskell
-- 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] ]
|