diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,22 @@ +The MIT License (MIT) + +Copyright (c) 2015 Dr. Lars Brünjes + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + diff --git a/Math/NumberTheory/Moduli/SquareRoots.hs b/Math/NumberTheory/Moduli/SquareRoots.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/SquareRoots.hs @@ -0,0 +1,50 @@ +-- | +-- Module: Math.NumberTheory.Moduli.SquareRoots +-- Copyright: (c) 2015 by Dr. Lars Brünjes +-- Licence: MIT +-- Maintainer: Dr. Lars Brünjes +-- Stability: Provisional +-- Portability: portable +-- +-- This module provides a function to find all square roots of a number modulo another number. +module Math.NumberTheory.Moduli.SquareRoots ( + sqrts ) where + +import Data.List (sort) +import Math.NumberTheory.Primes.Factorisation (factorise) +import Math.NumberTheory.Moduli (chineseRemainder, sqrtModPPList) + +chineseRemainders :: [([Integer], Integer)] -> [Integer] +chineseRemainders = fst . go where + go :: [([Integer], Integer)] -> ([Integer], Integer) + go [] = ([0], 1) + go xsm = foldr1 f xsm where + f (xs, m) (ys, n) = (xys, lcm m n) where + xys = do + x <- xs + y <- ys + case chineseRemainder [(x, m), (y, n)] of + Just z -> return z + Nothing -> [] + +sqrtsPP :: Integer -> (Integer, Int) -> [Integer] +sqrtsPP 1 (2, 1) = [1] +sqrtsPP 0 (p, e) = takeWhile (< p ^ e) $ map (* q) [0..] where + q = p ^ f + f = (if even e then e else succ e) `div` 2 +sqrtsPP a (p, e) + | a `mod` p /= 0 = sqrtModPPList a (p, e) + | a `mod` (p * p) /= 0 = [] + | otherwise = do + x <- sqrtsPP (a `div` (p * p)) (p, e - 2) + takeWhile (< m) [p * x + i * p ^ (e - 1) | i <- [0..]] + where + m = p ^ e + +-- |@sqrts a m@ finds all square roots of @a@ modulo @m@, +-- where @a@ is an arbitrary integer and @m@ is a positive integer. +sqrts :: Integer -> Integer -> [Integer] +sqrts a m + | a < 0 = error $ "a must not be negative, but a == " ++ show a ++ " < 0." + | otherwise = sort $ chineseRemainders $ map f $ factorise m where + f (p, e) = (sqrtsPP (a `mod` p ^ e) (p, e), p ^ e) diff --git a/Math/NumberTheory/Moduli/SquareRoots/Test.hs b/Math/NumberTheory/Moduli/SquareRoots/Test.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/SquareRoots/Test.hs @@ -0,0 +1,57 @@ +module Math.NumberTheory.Moduli.SquareRoots.Test ( + prop_sqrtsPP, + prop_sqrts) where + +import Data.Numbers.Primes (primes) +import Math.NumberTheory.Moduli.SquareRoots (sqrts) +import Test.QuickCheck + +newtype PrimePower = PrimePower (Integer, Int) deriving (Show, Eq) + +evalPP :: PrimePower -> Integer +evalPP (PrimePower (p, e)) = p ^ e + +instance Arbitrary PrimePower where + + arbitrary = sized $ \size -> + flip suchThat ((<= (fromIntegral $ size + 2)) . evalPP) $ do + indexP <- choose (0, size) + e <- choose (1, 20) + return $ PrimePower (primes !! indexP, e) + + shrink (PrimePower (p, e)) = + map PrimePower $ [(p', e) | p' <- takeWhile (< p) primes] ++ [(p, e') | e' <- [1, e - 1]] + +newtype ProblemPP = ProblemPP (Integer, PrimePower) deriving (Show, Eq) + +instance Arbitrary ProblemPP where + + arbitrary = do + pp <- arbitrary + a <- choose (0, evalPP pp - 1) + return $ ProblemPP (a, pp) + + shrink (ProblemPP (a, pp)) = + map ProblemPP $ [(a', pp) | a' <- [0, a - 1]] ++ [(a, pp') | pp' <- shrink pp, evalPP pp' > a] + +newtype Problem = Problem (Integer, Integer) deriving (Show, Eq) + +instance Arbitrary Problem where + + arbitrary = do + m <- suchThat arbitrary (> 0) + a <- choose (0, m - 1) + return $ Problem (a, m) + + shrink (Problem (a, m)) = + [Problem (a', m) | a' <- shrink a, a' >= 0] ++ + [Problem (a, m') | m' <- shrink m, m' > a] + +naive :: Problem -> [Integer] +naive (Problem (a, m)) = [x | x <- [0 .. (m - 1)], (x * x - a) `mod` m == 0] + +prop_sqrts :: Problem -> Property +prop_sqrts p@(Problem (a, m)) = sqrts a m === naive p + +prop_sqrtsPP :: ProblemPP -> Property +prop_sqrtsPP (ProblemPP (a, pp)) = prop_sqrts (Problem (a, evalPP pp)) diff --git a/Math/NumberTheory/Pell.hs b/Math/NumberTheory/Pell.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Pell.hs @@ -0,0 +1,115 @@ +-- | +-- Module: Math.NumberTheory.Pell +-- Copyright: (c) 2015 by Dr. Lars Brünjes +-- Licence: MIT +-- Maintainer: Dr. Lars Brünjes +-- Stability: Provisional +-- Portability: portable +-- +-- This module provides a function to solve generalized Pell Equations, +-- using the "LMM Algorithm" described by John P. Robertson in +-- . +-- A /generalized Pell Equation/ is a diophantine equation of the form +-- @x^2 - dy^2 = n@, where @d@ is a positive integer which is not a square +-- and where @n@ is a non-zero integer. +-- We are looking for solutions @(x,y)@, where @x@ and @y@ are non-negative integers. +module Math.NumberTheory.Pell ( + Solution, + solve ) where + +import Control.Arrow ((***)) +import Data.List (sort, nub) +import Data.Ratio ((%)) +import Data.Set (toList) +import Math.NumberTheory.Moduli.SquareRoots (sqrts) +import Math.NumberTheory.Pell.PQa (PQa(..), period) +import Math.NumberTheory.Powers.Squares (isSquare, integerSquareRoot) +import Math.NumberTheory.Primes.Factorisation (divisors) + +fmzs :: Integer -> Integer -> [(Integer, Integer, [Integer])] +fmzs d n = map (\f -> let m = n `div` (f * f) in (f, m, zs m)) $ filter (\f -> (n `mod` (f * f)) == 0) $ toList $ divisors n where + + zs :: Integer -> [Integer] + zs 1 = [0] + zs (-1) = [0] + zs m = nub $ sort $ map norm $ sqrts d am where + am = abs m + am2 = floor (am % 2) + norm x = if x <= am2 then x else x - am + +getRS :: Integer -> Integer -> Integer -> Maybe (Integer, Integer) +getRS d m z = + let + (_, pqas) = period z (abs m) d + qrs = zipWith (\x y -> (q x, g y, b y)) (tail pqas) pqas + rss = map (\(_, r, s) -> (r, s)) $ filter (\(q', _, _) -> abs q' == 1) qrs + in + case rss of + [] -> Nothing + (rs : _) -> Just rs + +-- |Represents a solution to a generalized Pell Equation. +-- The first component is the value of x, +-- the second component that of y. +type Solution = (Integer, Integer) + +solvePlusOne :: Integer -> Solution +solvePlusOne d = case getRS d 1 0 of + Just (r, s) + | r * r - d * s * s == 1 -> (r, s) + | otherwise -> (r * r + d * s * s, let rs = r * s in rs + rs) + Nothing -> error "algorithm error" + +solveMinusOne :: Integer -> Maybe Solution +solveMinusOne d = case getRS d 1 0 of + Just (r, s) + | r * r - d * s * s == (-1) -> Just (r, s) + | otherwise -> Nothing + Nothing -> error "algorithm error" + +getRep :: Integer -> Integer -> Integer -> Integer -> Integer -> Solution -> Maybe Solution -> Maybe Solution +getRep _ 1 _ _ _ x _ = Just x +getRep _ (-1) _ _ _ _ y = y +getRep d _ f m z _ y = case getRS d m z of + Nothing -> Nothing + Just (r, s) + | r * r - d * s * s == m -> Just (f * r, f * s) + | otherwise -> case y of + Nothing -> Nothing + Just (t, u) -> Just (f * (r * t + s * u * d), f * (r * u + s * t)) + +mul :: Integer -> Solution -> Solution -> Solution +mul d (x, y) (r, s) = (x * r + y * s * d, x * s + y * r) + +getReps :: Integer -> Integer -> (Solution, [Solution]) +getReps d n = ((r, s), reps) where + (r, s) = solvePlusOne d + minusOne = solveMinusOne d + reps = do + (f, m, zs) <- fmzs d n + do + z <- zs + case getRep d n f m z (r, s) minusOne of + Just (x, y) -> [if x >= 0 then (x, y) else mul d (r, s) (x, y)] + Nothing -> [] + +getMinimalReps :: Integer -> Integer -> (Solution, [Solution]) +getMinimalReps d n = ((r, s), map toMinimal reps) where + ((r, s), reps) = getReps d n + toMinimal (x, y) = minimum $ map (abs *** abs) $ filter (\(x', y') -> x' * y' >= 0) [(x, y), mul d (r, s) (x, y), mul d (r, -s) (x, y)] + +-- |@solve d n@ calculates all non-negative integer solutions of the generalized Pell Equation +-- x^2 - @d@y^2 = @n@, +-- where @d@ must be a positive integer which is not a square, +-- and @n@ must be a non-zero integer. +solve :: Integer -> Integer -> [Solution] +solve d n + | d <= 0 = error $ "D must be positive, but D == " ++ show d ++ "." + | isSquare d = error $ "D must not be a square, but D == " ++ show (integerSquareRoot d) ++ "^2." + | n == 0 = error "N must not be zero." + | otherwise = case getMinimalReps d n of + (_, []) -> [] + ((r, s), xys) -> go xys where + go xys' = normalize xys' ++ go (step xys') + normalize = sort . nub + step = map (mul d (r, s)) diff --git a/Math/NumberTheory/Pell/PQa.hs b/Math/NumberTheory/Pell/PQa.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Pell/PQa.hs @@ -0,0 +1,57 @@ +module Math.NumberTheory.Pell.PQa ( + PQa(..), + pqa, + reduced, + period) where + +import Data.Ratio ((%)) +import Math.NumberTheory.Powers.Squares (isSquare, integerSquareRoot) + +data PQa = PQa { + a :: Integer, + b :: Integer, + g :: Integer, + a' :: Integer, + p :: Integer, + q :: Integer } deriving Show + +pqa :: Integer -> Integer -> Integer -> [PQa] +pqa p0 q0 d + | q0 == 0 = error "Q0 must not be zero." + | d <= 0 = error "D must be positive." + | isSquare d = error $ "D must not be a square, but D == " ++ show dd ++ "^2." + | (p0 * p0 - d) `mod` q0 /= 0 = error $ "P0^2 must be equivalent to D modulo Q0, but " + ++ show p0 ++ "^2 == " ++ show (p0 `mod` q0) ++ " /= " ++ show (d `mod` q0) + ++ " == " ++ show d ++ " (mod " ++ show q0 ++ ")" + | otherwise = go p0 q0 (PQa 0 1 (-p0) undefined undefined undefined) (PQa 1 0 q0 undefined undefined undefined) + where + dd = integerSquareRoot d + go :: Integer -> Integer -> PQa -> PQa -> [PQa] + go p' q' x y = + let + _a' = if q' > 0 then floor $ (p' + dd) % q' else floor $ (p' + dd + 1) % q' + _a = _a' * a y + a x + _b = _a' * b y + b x + _g = _a' * g y + g x + _p = _a' * q' - p' + _q = (d - _p * _p) `div` q' + z = PQa _a _b _g _a' p' q' + in + z : go _p _q y z + +reduced :: Integer -> Integer -> Integer -> Bool +reduced x y dd + | y > 0 = (dd >= y - x) && (dd < x + y) && (dd >= x) + | otherwise = (dd < y - x) && (dd >= x + y) && (dd < x) + +period :: Integer -> Integer -> Integer -> (Int, [PQa]) +period p0 q0 d = u [] 0 $ pqa p0 q0 d where + dd = integerSquareRoot d + u acc i (x : xs) + | reduced (p x) (q x) dd = v (x : acc) i (p x) (q x) xs + | otherwise = u (x : acc) (succ i) xs + u _ _ _ = error "algorithm error" + v acc i pp qq (x : xs) + | (pp == p x) && (qq == q x) = (i, reverse acc) + | otherwise = v (x : acc) i pp qq xs + v _ _ _ _ _ = error "algorithm error" diff --git a/Math/NumberTheory/Pell/Test.hs b/Math/NumberTheory/Pell/Test.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Pell/Test.hs @@ -0,0 +1,35 @@ +module Math.NumberTheory.Pell.Test where + +import Distribution.TestSuite.QuickCheck (Test, testProperty, testGroup) +import Math.NumberTheory.Moduli.SquareRoots.Test (prop_sqrtsPP, prop_sqrts) +import Math.NumberTheory.Pell.Test.Reduced (prop_reduced) +import Math.NumberTheory.Pell.Test.Solve (Problem (..), prop_solves) + +tests :: IO [Test] +tests = return + [ + testGroup "SquareRoots" + [ + testProperty "sqrtsPP" prop_sqrtsPP, + testProperty "sqrts" prop_sqrts + ], + testGroup "Pell" + [ + testProperty "reduced" prop_reduced, + testProperty "solves 7 9" (prop_solves 100 $ Problem 7 9), + testProperty "solves 5 (-4)" (prop_solves 100 $ Problem 5 (-4)), + testProperty "solves 2 (-7)" (prop_solves 100 $ Problem 2 (-7)), + testProperty "solves" (prop_solves 100000) + ] + ] + +-- main :: IO () +-- main = do +-- test prop_sqrtsPP +-- test prop_sqrts +-- test prop_reduced +-- test (prop_solves 100 $ Problem 7 9) +-- test (prop_solves 100 $ Problem 5 (-4)) +-- test (prop_solves 100 $ Problem 2 (-7)) +-- test (prop_solves 100000) +-- diff --git a/Math/NumberTheory/Pell/Test/Reduced.hs b/Math/NumberTheory/Pell/Test/Reduced.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Pell/Test/Reduced.hs @@ -0,0 +1,51 @@ +module Math.NumberTheory.Pell.Test.Reduced ( prop_reduced ) where + +import Math.NumberTheory.Pell.PQa (reduced) +import Math.NumberTheory.Powers.Squares (isSquare, integerSquareRoot) +import Test.QuickCheck + +data Triple = Triple Integer Integer Integer deriving (Show, Eq) + +isProperTriple :: Triple -> Bool +isProperTriple (Triple _ q d) = (q /= 0) && (d > 0) && not (isSquare d) + +toDouble :: Triple -> (Double, Double) +toDouble (Triple p q d) = + let + pp = fromIntegral p + qq = fromIntegral q + dd = sqrt $ fromIntegral d + x = (pp + dd) / qq + y = (pp - dd) / qq + in + (x, y) + +isReduced :: Triple -> Bool +isReduced t = let (x, y) = toDouble t in (x > 1) && (-1 < y) && (y < 0) + +genTriple :: Gen Triple +genTriple = flip suchThat isProperTriple $ do + p <- scale (* 2) arbitrary + q <- scale (* 2) arbitrary + d <- scale (* 3) arbitrary + return $ Triple p q d + +instance Arbitrary Triple where + arbitrary = oneof $ map (suchThat genTriple) [isReduced, not . isReduced] + shrink (Triple p q d) = filter isProperTriple $ + [Triple p' q d | p' <- shrink p] ++ + [Triple p q' d | q' <- shrink q] ++ + [Triple p q d' | d' <- shrink d] + +reduced' :: Triple -> Bool +reduced' (Triple p q d) = reduced p q (integerSquareRoot d) + +prop_reduced :: Triple -> Property +prop_reduced t@(Triple _ _ d) = + counterexample (show $ toDouble t) $ + classify (reduced' t) "reduced" $ + classify (not $ reduced' t) "not reduced" $ + classify (d <= 100) "d <= 100" $ + classify (d > 100) "d > 100" $ + reduced' t === isReduced t + diff --git a/Math/NumberTheory/Pell/Test/Solve.hs b/Math/NumberTheory/Pell/Test/Solve.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Pell/Test/Solve.hs @@ -0,0 +1,52 @@ +module Math.NumberTheory.Pell.Test.Solve ( + Problem (..), + prop_solves, + naive) where + +import Control.Monad (liftM2) +import Math.NumberTheory.Pell (solve) +import Math.NumberTheory.Powers.Squares (isSquare, integerSquareRoot) +import Test.QuickCheck + +data Problem = Problem Integer Integer deriving (Show, Eq) + +isProperD :: Integral a => a -> Bool +isProperD n = (n > 0) && not (isSquare n) + +genD :: Gen Integer +genD = scale (* 2) $ suchThat arbitrary isProperD + +shrinkD :: Integer -> [Integer] +shrinkD d = filter isProperD $ shrink d + +genN :: Gen Integer +genN = scale (* 2) $ sized genN' where + genN' 0 = elements [-1, 1] + genN' 1 = elements [-4, 4] + genN' _ = do + x <- suchThat arbitrary (> 0) + let y = integerSquareRoot x + elements [-y, y] + +shrinkN :: Integer -> [Integer] +shrinkN n = filter (/= 0) $ shrink n + +instance Arbitrary Problem where + arbitrary = liftM2 Problem genD genN + shrink (Problem d n) = [Problem d' n | d' <- shrinkD d] ++ [Problem d n' | n' <- shrinkN n] + +naive :: Integer -> Integer -> Integer -> [(Integer, Integer)] +naive maxY d n = [(integerSquareRoot $ n + d * y * y, y) | y <- [0..maxY], isSquare $ n + d * y * y] + +prop_solves :: Integer -> Problem -> Property +prop_solves limit (Problem d n) = + classify (n == 1) "n == 1" $ + classify (n == (-1)) "n == -1" $ + classify (n == 4) "n == 4" $ + classify (n == (-4)) "n == -4" $ + classify (abs n `notElem` [1, 4]) "|n| /= 1, 4" $ + classify (n * n < d) "n^2 < d" $ + classify (n * n > d) "n^2 > d" $ + classify (d <= 100) "d <= 100" $ + classify (d > 100) "d > 100" $ + takeWhile ((<= limit) . snd) (solve d n) === naive limit d n diff --git a/Math/NumberTheory/Pell/Test/Utils.hs b/Math/NumberTheory/Pell/Test/Utils.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Pell/Test/Utils.hs @@ -0,0 +1,9 @@ +module Math.NumberTheory.Pell.Test.Utils ( + (~~) ) where + +import Test.QuickCheck + +infix 4 ~~ + +(~~) :: Double -> Double -> Property +x ~~ y = counterexample ("expected " ++ show x ++ " ~~ " ++ show y) $ abs (x - y) < 1e-4 diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# pell +Haskell Package to solve the Generalized Pell Equation diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/dist/build/test-pellStub/test-pellStub-tmp/test-pellStub.hs b/dist/build/test-pellStub/test-pellStub-tmp/test-pellStub.hs new file mode 100644 --- /dev/null +++ b/dist/build/test-pellStub/test-pellStub-tmp/test-pellStub.hs @@ -0,0 +1,5 @@ +module Main ( main ) where +import Distribution.Simple.Test.LibV09 ( stubMain ) +import Math.NumberTheory.Pell.Test ( tests ) +main :: IO () +main = stubMain tests diff --git a/pell.cabal b/pell.cabal new file mode 100644 --- /dev/null +++ b/pell.cabal @@ -0,0 +1,55 @@ +-- Initial pell.cabal generated by cabal init. For further documentation, +-- see http://haskell.org/cabal/users-guide/ + +name: pell +version: 0.1.0.0 +synopsis: Package to solve the Generalized Pell Equation. +description: Finds all solutions of the generalized Pell Equation. +homepage: https://github.com/brunjlar/pell +license: MIT +license-file: LICENSE +author: Lars Bruenjes +maintainer: lbrunjes@gmx.de +copyright: (c) 2015 by Dr. Lars Brünjes +category: Math, Algorithms, Number Theory +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.20.0 + +library + exposed-modules: Math.NumberTheory.Pell, Math.NumberTheory.Moduli.SquareRoots + other-modules: Math.NumberTheory.Pell.PQa + -- other-extensions: + build-depends: base >=4.7 && <4.8, + arithmoi, + containers + -- hs-source-dirs: + default-language: Haskell2010 + +Test-Suite test-pell + type: detailed-0.9 + test-module: Math.NumberTheory.Pell.Test + other-modules: Math.NumberTheory.Moduli.SquareRoots, + Math.NumberTheory.Moduli.SquareRoots.Test, + Math.NumberTheory.Pell, + Math.NumberTheory.Pell.PQa, + Math.NumberTheory.Pell.Test.Reduced, + Math.NumberTheory.Pell.Test.Solve, + Math.NumberTheory.Pell.Test.Utils + build-depends: base >= 4.7 && <4.8, + arithmoi, + containers, + QuickCheck >= 2.8, + primes, + Cabal >= 1.20.0, + cabal-test-quickcheck + default-language: Haskell2010 + +source-repository head + type: git + location: https://github.com/brunjlar/pell + +source-repository this + type: git + location: https://github.com/brunjlar/pell + tag: 0.1.0.0