packages feed

pell (empty) → 0.1.0.0

raw patch · 13 files changed

+512/−0 lines, 13 filesdep +Cabaldep +QuickCheckdep +arithmoisetup-changed

Dependencies added: Cabal, QuickCheck, arithmoi, base, cabal-test-quickcheck, containers, primes

Files

+ LICENSE view
@@ -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.+
+ Math/NumberTheory/Moduli/SquareRoots.hs view
@@ -0,0 +1,50 @@+-- |+-- Module:      Math.NumberTheory.Moduli.SquareRoots+-- Copyright:   (c) 2015 by Dr. Lars Brünjes+-- Licence:     MIT+-- Maintainer:  Dr. Lars Brünjes <lbrunjes@gmx.de>+-- 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)
+ Math/NumberTheory/Moduli/SquareRoots/Test.hs view
@@ -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))
+ Math/NumberTheory/Pell.hs view
@@ -0,0 +1,115 @@+-- |+-- Module:      Math.NumberTheory.Pell+-- Copyright:   (c) 2015 by Dr. Lars Brünjes+-- Licence:     MIT+-- Maintainer:  Dr. Lars Brünjes <lbrunjes@gmx.de>+-- Stability:   Provisional+-- Portability: portable+--+-- This module provides a function to solve generalized Pell Equations, +-- using the "LMM Algorithm" described by John P. Robertson in+-- <http://http://www.jpr2718.org/pell.pdf>.+-- 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))
+ Math/NumberTheory/Pell/PQa.hs view
@@ -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"
+ Math/NumberTheory/Pell/Test.hs view
@@ -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)+--
+ Math/NumberTheory/Pell/Test/Reduced.hs view
@@ -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+
+ Math/NumberTheory/Pell/Test/Solve.hs view
@@ -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
+ Math/NumberTheory/Pell/Test/Utils.hs view
@@ -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
+ README.md view
@@ -0,0 +1,2 @@+# pell+Haskell Package to solve the Generalized Pell Equation
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ dist/build/test-pellStub/test-pellStub-tmp/test-pellStub.hs view
@@ -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
+ pell.cabal view
@@ -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