Merge branch 'master' of github:yogsototh/euler

066.hs

@@ -22,27 +22,84 @@ import Control.Monad
 import Data.List
 import Data.Maybe
 
-isqrt = floor . sqrt . fromIntegral
+-- 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
 
-hasSolution d x =
-    let y = isqrt (( x^2 - 1 ) `div` d)
-    in x^2 - (y^2)*d == 1
+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)
 
-phead [] = -2 
-phead (x:xs) = x
+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
 
-sol d = if isSquare d
+-- 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 phead $ filter (hasSolution d) [2..1000000]
+               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 = 100
+            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 !! (d-2))
+           forM_ [2..m] (\d -> do putStr $ (show d) ++ " -> "; print $ lst !! (fromIntegral d-2))
            putStr "Maximum value: "
            print maxelem
            putStr "Maximum D: "