106 lines
3.2 KiB
Haskell
106 lines
3.2 KiB
Haskell
-- 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
|
||
|
||
-- 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
|
||
|
||
isqrt = floor . sqrt . fromIntegral
|
||
isSquare d = (isqrt d)^2 == d
|
||
|
||
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
|
||
|
||
-- 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)
|
||
|
||
-- 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)
|
||
|
||
main = let
|
||
m = 1000
|
||
lst = map sol [2..m]
|
||
maxelem = foldl' (\x y -> if x>y then x else y) 0 lst
|
||
in do
|
||
forM_ [2..m] (\d -> do putStr $ (show d) ++ " -> "; print $ lst !! (fromIntegral d-2))
|
||
putStr "Maximum value: "
|
||
print maxelem
|
||
putStr "Maximum D: "
|
||
print $ 2 + (fromJust $ elemIndex maxelem lst)
|