2011-12-19 16:55:33 +00:00
|
|
|
|
-- Consider quadratic Diophantine equations of the form:
|
|
|
|
|
--
|
|
|
|
|
-- x^2 – Dy^2 = 1
|
|
|
|
|
--
|
|
|
|
|
-- For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.
|
|
|
|
|
--
|
|
|
|
|
-- It can be assumed that there are no solutions in positive integers when D is square.
|
|
|
|
|
--
|
|
|
|
|
-- By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:
|
|
|
|
|
--
|
|
|
|
|
-- 32 – 2×22 = 1
|
|
|
|
|
-- 22 – 3×12 = 1
|
|
|
|
|
-- 92 – 5×42 = 1
|
|
|
|
|
-- 52 – 6×22 = 1
|
|
|
|
|
-- 82 – 7×32 = 1
|
|
|
|
|
--
|
|
|
|
|
-- Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.
|
|
|
|
|
--
|
|
|
|
|
-- Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.
|
|
|
|
|
|
|
|
|
|
import Control.Monad
|
|
|
|
|
import Data.List
|
|
|
|
|
import Data.Maybe
|
|
|
|
|
|
2011-12-19 21:02:08 +00:00
|
|
|
|
-- import Debug.Trace
|
|
|
|
|
|
|
|
|
|
-- -- Brute force search
|
|
|
|
|
-- isqrt = floor . sqrt . fromIntegral
|
|
|
|
|
--
|
|
|
|
|
-- isSquare d = (isqrt d)^2 == d
|
|
|
|
|
--
|
|
|
|
|
--
|
|
|
|
|
--
|
|
|
|
|
-- hasSolution d x =
|
|
|
|
|
-- let y = isqrt (( x^2 - 1 ) `div` d)
|
|
|
|
|
-- in x^2 - (y^2)*d == 1
|
|
|
|
|
--
|
|
|
|
|
-- phead [] = -2
|
|
|
|
|
-- phead (x:xs) = x
|
|
|
|
|
--
|
|
|
|
|
-- sol d = if isSquare d
|
|
|
|
|
-- then -1
|
|
|
|
|
-- else phead $ filter (hasSolution d) [2..10000000]
|
|
|
|
|
|
|
|
|
|
-- Using Pell equation formulae
|
|
|
|
|
-- http://fr.wikipedia.org/wiki/Équation_de_Pell-Fermat#Cas_x.C2.B2-ny.C2.B2_.3D_1
|
2011-12-19 16:55:33 +00:00
|
|
|
|
|
2011-12-19 21:02:08 +00:00
|
|
|
|
isqrt = floor . sqrt . fromIntegral
|
2011-12-19 16:55:33 +00:00
|
|
|
|
isSquare d = (isqrt d)^2 == d
|
|
|
|
|
|
2011-12-19 21:02:08 +00:00
|
|
|
|
ratFracSqrt :: Integer -> Integer -> Integer -> [(Integer,(Integer,Integer,Integer))]
|
|
|
|
|
ratFracSqrt a b c =
|
|
|
|
|
let m = a - (b^2)
|
|
|
|
|
fa = fromInteger a :: Float
|
|
|
|
|
fb = fromInteger b :: Float
|
|
|
|
|
fc = fromInteger c :: Float
|
|
|
|
|
fm = fromInteger m :: Float
|
|
|
|
|
x = floor $ (fc * ((sqrt fa) - fb))/ (fm)
|
|
|
|
|
h = c
|
|
|
|
|
i = -1 * ( c*b + x*m )
|
|
|
|
|
j = m
|
|
|
|
|
in
|
|
|
|
|
if ( x^2 == a )
|
|
|
|
|
then (0,(0,0,0)):ratFracSqrt a b c
|
|
|
|
|
else if (m `rem` h /= 0)
|
|
|
|
|
then error "Bad division"
|
|
|
|
|
else (x,(a,b,c)):ratFracSqrt a (i `div` h) (j `div` h)
|
|
|
|
|
|
|
|
|
|
ratFracPeriodLength s =
|
|
|
|
|
let
|
|
|
|
|
sig = head $ tail s
|
|
|
|
|
index = fromJust $ elemIndex sig $( tail $ tail s )
|
|
|
|
|
in if head s == (0,(0,0,0))
|
|
|
|
|
then 0
|
|
|
|
|
else index + 1
|
2011-12-19 16:55:33 +00:00
|
|
|
|
|
2011-12-19 21:02:08 +00:00
|
|
|
|
-- The "reduite" is the number p/q = a_0 + ( 1/ (a_1 + ...) )
|
|
|
|
|
reduite (x:[]) = (x,1)
|
|
|
|
|
reduite (x:xs) = let
|
|
|
|
|
precalc = reduite xs
|
|
|
|
|
num = fst precalc
|
|
|
|
|
den = snd precalc
|
|
|
|
|
in
|
|
|
|
|
(x * num + den, num)
|
2011-12-19 16:55:33 +00:00
|
|
|
|
|
2011-12-19 21:02:08 +00:00
|
|
|
|
-- Solution for x^2 - dy^2 = 1
|
|
|
|
|
sol d = let
|
|
|
|
|
rfs = ratFracSqrt d 0 d
|
|
|
|
|
len = ratFracPeriodLength rfs
|
|
|
|
|
usedLen = if len `rem` 2 == 0 then (len) else (2*len)
|
|
|
|
|
in
|
|
|
|
|
if isSquare d
|
|
|
|
|
then -1
|
|
|
|
|
else fst $ reduite $ take usedLen $ map fst rfs
|
|
|
|
|
-- else fst $ reduite $ trace (show $ take usedLen $ map fst rfs) (take usedLen $ map fst rfs)
|
2011-12-19 16:55:33 +00:00
|
|
|
|
|
|
|
|
|
main = let
|
2011-12-19 21:02:08 +00:00
|
|
|
|
m = 1000
|
2011-12-19 17:02:19 +00:00
|
|
|
|
lst = map sol [2..m]
|
2011-12-19 16:55:33 +00:00
|
|
|
|
maxelem = foldl' (\x y -> if x>y then x else y) 0 lst
|
|
|
|
|
in do
|
2011-12-19 21:02:08 +00:00
|
|
|
|
forM_ [2..m] (\d -> do putStr $ (show d) ++ " -> "; print $ lst !! (fromIntegral d-2))
|
2011-12-19 17:02:19 +00:00
|
|
|
|
putStr "Maximum value: "
|
2011-12-19 16:55:33 +00:00
|
|
|
|
print maxelem
|
2011-12-19 17:02:19 +00:00
|
|
|
|
putStr "Maximum D: "
|
|
|
|
|
print $ 2 + (fromJust $ elemIndex maxelem lst)
|