quadratic-irrational 0.0.6 → 0.1.0
raw patch · 9 files changed
+412/−297 lines, 9 filesdep +semigroupsdep ~arithmoidep ~mtlPVP ok
version bump matches the API change (PVP)
Dependencies added: semigroups
Dependency ranges changed: arithmoi, mtl
API changes (from Hackage documentation)
- Numeric.QuadraticIrrational: instance GHC.Classes.Ord Numeric.QuadraticIrrational.QI
- Numeric.QuadraticIrrational: qiOne :: QI
- Numeric.QuadraticIrrational: qiZero :: QI
Files
- .travis.yml +80/−60
- ChangeLog.md +7/−0
- quadratic-irrational.cabal +10/−5
- src/Numeric/QuadraticIrrational.hs +191/−204
- src/Numeric/QuadraticIrrational/CyclicList.hs +6/−1
- src/Numeric/QuadraticIrrational/Ring.hs +48/−0
- src/Numeric/QuadraticIrrational/SquareFree.hs +48/−0
- tests/QuadraticIrrational.hs +7/−16
- tests/doctests.hs +15/−11
.travis.yml view
@@ -1,114 +1,134 @@+# This Travis job script has been generated by a script via+#+# haskell-ci '--config' 'cabal.haskell-ci' '--output' '.travis.yml' 'quadratic-irrational.cabal'+#+# For more information, see https://github.com/haskell-CI/haskell-ci+#+# version: 0.2.1+# language: c-sudo: false+dist: xenial +git:+ submodules: false # whether to recursively clone submodules+ cache: directories: - $HOME/.cabal/packages - $HOME/.cabal/store before_cache:- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/build-reports.log+ - rm -fv $CABALHOME/packages/hackage.haskell.org/build-reports.log # remove files that are regenerated by 'cabal update'- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/00-index.*- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/*.json- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.cache- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar.idx+ - rm -fv $CABALHOME/packages/hackage.haskell.org/00-index.*+ - rm -fv $CABALHOME/packages/hackage.haskell.org/*.json+ - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.cache+ - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar+ - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar.idx - - rm -rfv $HOME/.cabal/packages/head.hackage+ - rm -rfv $CABALHOME/packages/head.hackage matrix: include:+ - compiler: "ghc-8.8.1"+ env: GHCHEAD=true+ addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.8.1], sources: [hvr-ghc]}}+ - compiler: "ghc-8.6.5"+ addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.6.5], sources: [hvr-ghc]}}+ - compiler: "ghc-8.4.4"+ addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.4.4], sources: [hvr-ghc]}}+ - compiler: "ghc-8.2.2"+ addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.2.2], sources: [hvr-ghc]}}+ - compiler: "ghc-8.0.2"+ addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.0.2], sources: [hvr-ghc]}} - compiler: "ghc-7.10.3" env: TEST=--disable-tests addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-7.10.3], sources: [hvr-ghc]}}- - compiler: "ghc-8.0.2"- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.0.2], sources: [hvr-ghc]}}- - compiler: "ghc-8.2.2"- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.2.2], sources: [hvr-ghc]}}- - compiler: "ghc-8.4.3"- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.4.3], sources: [hvr-ghc]}}- - compiler: "ghc-8.6.1"- env: GHCHEAD=true- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.6.1], sources: [hvr-ghc]}} - compiler: "ghc-head" env: GHCHEAD=true addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-head], sources: [hvr-ghc]}} allow_failures:+ - compiler: "ghc-8.8.1" - compiler: "ghc-head" before_install:- - HC=${CC}+ - HC=/opt/ghc/bin/${CC} - HCPKG=${HC/ghc/ghc-pkg} - unset CC+ - CABAL=/opt/ghc/bin/cabal+ - CABALHOME=$HOME/.cabal+ - export PATH="$CABALHOME/bin:$PATH" - ROOTDIR=$(pwd)- - mkdir -p $HOME/.local/bin- - "PATH=/opt/ghc/bin:/opt/ghc-ppa-tools/bin:$HOME/local/bin:$PATH" - HCNUMVER=$(( $(${HC} --numeric-version|sed -E 's/([0-9]+)\.([0-9]+)\.([0-9]+).*/\1 * 10000 + \2 * 100 + \3/') )) - echo $HCNUMVER install:- - cabal --version+ - ${CABAL} --version - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"- - BENCH=${BENCH---enable-benchmarks} - TEST=${TEST---enable-tests}- - HADDOCK=${HADDOCK-true}- - INSTALLED=${INSTALLED-true}+ - BENCH=--enable-benchmarks - GHCHEAD=${GHCHEAD-false}- - travis_retry cabal update -v- - "sed -i.bak 's/^jobs:/-- jobs:/' ${HOME}/.cabal/config"+ - travis_retry ${CABAL} update -v+ - sed -i.bak 's/^jobs:/-- jobs:/' $CABALHOME/config - rm -fv cabal.project cabal.project.local- # Overlay Hackage Package Index for GHC HEAD: https://github.com/hvr/head.hackage - | if $GHCHEAD; then- sed -i.bak 's/-- allow-newer:.*/allow-newer: *:base, *:template-haskell, *:ghc, *:Cabal, *:stm/' ${HOME}/.cabal/config-- echo 'repository head.hackage' >> ${HOME}/.cabal/config- echo ' url: http://head.hackage.haskell.org/' >> ${HOME}/.cabal/config- echo ' secure: True' >> ${HOME}/.cabal/config- echo ' root-keys: 07c59cb65787dedfaef5bd5f987ceb5f7e5ebf88b904bbd4c5cbdeb2ff71b740' >> ${HOME}/.cabal/config- echo ' 2e8555dde16ebd8df076f1a8ef13b8f14c66bad8eafefd7d9e37d0ed711821fb' >> ${HOME}/.cabal/config- echo ' 8f79fd2389ab2967354407ec852cbe73f2e8635793ac446d09461ffb99527f6e' >> ${HOME}/.cabal/config- echo ' key-threshold: 3' >> ${HOME}/.cabal.config-- cabal new-update head.hackage -v- fi- - grep -Ev -- '^\s*--' ${HOME}/.cabal/config | grep -Ev '^\s*$'- - "printf 'packages: \".\"\\n' > cabal.project"- - cat cabal.project- - if [ -f "./configure.ac" ]; then- (cd "." && autoreconf -i);+ sed -i.bak 's/-- allow-newer:.*/allow-newer: *:base, *:template-haskell, *:ghc, *:Cabal/' $CABALHOME/config fi+ - grep -Ev -- '^\s*--' $CABALHOME/config | grep -Ev '^\s*$'+ - rm -f cabal.project+ - touch cabal.project+ - "printf 'packages: \".\"\\n' >> cabal.project"+ - "printf 'write-ghc-environment-files: always\\n' >> cabal.project"+ - touch cabal.project.local+ - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | grep -vE -- '^(quadratic-irrational)$' | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"+ - cat cabal.project || true+ - cat cabal.project.local || true+ - if [ -f "./configure.ac" ]; then (cd "." && autoreconf -i); fi - rm -f cabal.project.freeze- - cabal new-build -w ${HC} ${TEST} ${BENCH} --project-file="cabal.project" --dep -j2 all- - cabal new-build -w ${HC} --disable-tests --disable-benchmarks --project-file="cabal.project" --dep -j2 all- - rm -rf "."/.ghc.environment.* "."/dist+ - ${CABAL} new-freeze -w ${HC} ${TEST} ${BENCH} --project-file="cabal.project" --dry+ - "cat \"cabal.project.freeze\" | sed -E 's/^(constraints: *| *)//' | sed 's/any.//'"+ - rm "cabal.project.freeze"+ - ${CABAL} new-build -w ${HC} ${TEST} ${BENCH} --project-file="cabal.project" --dep -j2 all+ - ${CABAL} new-build -w ${HC} --disable-tests --disable-benchmarks --project-file="cabal.project" --dep -j2 all+ - rm -rf .ghc.environment.* "."/dist - DISTDIR=$(mktemp -d /tmp/dist-test.XXXX) # Here starts the actual work to be performed for the package under test; # any command which exits with a non-zero exit code causes the build to fail. script: # test that source-distributions can be generated- - (cd "." && cabal sdist)- - mv "."/dist/quadratic-irrational-*.tar.gz ${DISTDIR}/+ - ${CABAL} new-sdist all+ - mv dist-newstyle/sdist/*.tar.gz ${DISTDIR}/ - cd ${DISTDIR} || false - find . -maxdepth 1 -name '*.tar.gz' -exec tar -xvf '{}' \;- - "printf 'packages: quadratic-irrational-*/*.cabal\\n' > cabal.project"- - cat cabal.project+ - rm -f cabal.project+ - touch cabal.project+ - "printf 'packages: \"quadratic-irrational-*/*.cabal\"\\n' >> cabal.project"+ - "printf 'write-ghc-environment-files: always\\n' >> cabal.project"+ - touch cabal.project.local+ - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | grep -vE -- '^(quadratic-irrational)$' | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"+ - cat cabal.project || true+ - cat cabal.project.local || true # this builds all libraries and executables (without tests/benchmarks)- - cabal new-build -w ${HC} --disable-tests --disable-benchmarks -f check-bounds all-- # Build with installed constraints for packages in global-db- - if $INSTALLED; then echo cabal new-build -w ${HC} --disable-tests --disable-benchmarks $(${HCPKG} list --global --simple-output --names-only | sed 's/\([a-zA-Z0-9-]\{1,\}\) */--constraint="\1 installed" /g') all | sh; else echo "Not building with installed constraints"; fi+ - ${CABAL} new-build -w ${HC} --disable-tests --disable-benchmarks all # build & run tests, build benchmarks- - cabal new-build -w ${HC} ${TEST} ${BENCH} all- - if [ "x$TEST" = "x--enable-tests" ]; then cabal new-test -w ${HC} ${TEST} ${BENCH} -f check-bounds all; fi+ - ${CABAL} new-build -w ${HC} ${TEST} ${BENCH} all+ - if [ "x$TEST" = "x--enable-tests" ]; then ${CABAL} new-test -w ${HC} ${TEST} ${BENCH} all; fi + # cabal check+ - (cd quadratic-irrational-* && ${CABAL} check)+ # haddock- - rm -rf ./dist-newstyle- - if $HADDOCK; then cabal new-haddock -w ${HC} ${TEST} ${BENCH} all; else echo "Skipping haddock generation";fi+ - ${CABAL} new-haddock -w ${HC} ${TEST} ${BENCH} all -# REGENDATA ["quadratic-irrational.cabal"]+ # Build without installed constraints for packages in global-db+ - rm -f cabal.project.local; ${CABAL} new-build -w ${HC} --disable-tests --disable-benchmarks all;++ # Constraint sets+ - rm -rf cabal.project.local++# REGENDATA ["--config","cabal.haskell-ci","--output",".travis.yml","quadratic-irrational.cabal"] # EOF
ChangeLog.md view
@@ -1,3 +1,10 @@+# 0.1.0 (2019-04-26)++* Allow imaginary square roots, e. g., `qi 1 1 (-5) 1`.+* Remove `Ord QI` instance: complex values cannot be ordered.+* Roots of 0 are reduced: `qi a b 0 d` becomes `qi a 0 2 d`.+* Remove `qiZero` and `qiOne`.+ # 0.0.6 (2018-08-29) * Support GHC up to 8.6.1.
quadratic-irrational.cabal view
@@ -1,6 +1,6 @@ name: quadratic-irrational category: Math, Algorithms, Data-version: 0.0.6+version: 0.1.0 license: MIT license-file: LICENSE author: Johan Kiviniemi <devel@johan.kiviniemi.name>@@ -51,7 +51,7 @@ > 3.5 = [3; 2] > (1+√5)/2 = [1; 1, 1, 1, …] > √2 = [1; 2, 2, 2, …]-tested-with: GHC==7.10.3, GHC==8.0.2, GHC==8.2.2, GHC==8.4.3, GHC==8.6.1+tested-with: GHC==7.10.3, GHC==8.0.2, GHC==8.2.2, GHC==8.4.3, GHC==8.6.5, GHC==8.8.1 build-type: Simple cabal-version: >= 1.10@@ -69,12 +69,16 @@ exposed-modules: Numeric.QuadraticIrrational , Numeric.QuadraticIrrational.CyclicList , Numeric.QuadraticIrrational.Internal.Lens+ other-modules: Numeric.QuadraticIrrational.SquareFree+ , Numeric.QuadraticIrrational.Ring hs-source-dirs: src build-depends: base >= 4.8 && < 5- , arithmoi >= 0.4+ , arithmoi >= 0.7 , containers >= 0.5 && < 0.7- , mtl >= 2.1 && < 2.3 , transformers >= 0.3 && < 0.6+ if impl(ghc <8.0)+ build-depends:+ semigroups >=0.8 default-language: Haskell2010 ghc-options: -Wall -O2 -funbox-strict-fields @@ -88,7 +92,7 @@ , quadratic-irrational , numbers >= 3000.0 && < 3000.3 , QuickCheck >= 2.7 && < 3- , tasty >= 0.8 && < 1.2+ , tasty >= 0.8 , tasty-quickcheck >= 0.8 && < 0.11 default-language: Haskell2010 ghc-options: -Wall -O2 -funbox-strict-fields@@ -103,5 +107,6 @@ , filepath , mtl , numbers+ , arithmoi default-language: Haskell2010 ghc-options: -threaded -Wall
src/Numeric/QuadraticIrrational.hs view
@@ -1,5 +1,8 @@-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE ViewPatterns #-} -- | -- Module : Numeric.QuadraticIrrational@@ -57,7 +60,7 @@ , -- * Lenses _qi, _qi', _qiABD, _qiA, _qiB, _qiC, _qiD , -- * Numerical operations- qiZero, qiOne, qiIsZero+ qiIsZero , qiToFloat , qiAddI, qiSubI, qiMulI, qiDivI , qiAddR, qiSubR, qiMulR, qiDivR@@ -69,12 +72,15 @@ , module Numeric.QuadraticIrrational.CyclicList ) where -import Control.Monad.State (evalState, gets, modify)+import Control.Monad.Trans.State (evalState, gets, modify) import Data.Foldable (toList) import Data.List (foldl') import Data.Maybe (fromMaybe)+import Data.Proxy import Data.Ratio ((%), denominator, numerator)+import Data.Semigroup import qualified Data.Set as Set (empty, insert, member)+import Data.Type.Equality import Math.NumberTheory.Powers.Squares (integerSquareRoot) import Math.NumberTheory.Primes.Factorisation (factorise) import Text.Read (Lexeme (Ident), Read (readListPrec, readPrec),@@ -82,19 +88,23 @@ import Numeric.QuadraticIrrational.CyclicList import Numeric.QuadraticIrrational.Internal.Lens+import Numeric.QuadraticIrrational.Ring+import Numeric.QuadraticIrrational.SquareFree -- $setup -- >>> import Data.Number.CReal --- | @(a + b √c) \/ d@-data QI = QI !Integer- !Integer- !Integer- !Integer- deriving (Eq)+-- | @(a + b √c)@+data QI where+ QI :: KnownSquareFree c => QuadExt c Rational -> QI +instance Eq QI where+ (QI (q1 :: QuadExt c1 Rational)) == (QI (q2 :: QuadExt c2 Rational)) = case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of+ Nothing -> False+ Just Refl -> q1 == q2+ instance Show QI where- showsPrec p (QI a b c d) = showParen (p > 10)+ showsPrec p q = runQI q $ \a b c d -> showParen (p > 10) $ showString "qi " . showsPrec 11 a . showChar ' ' . showsPrec 11 b . showChar ' ' . showsPrec 11 c@@ -110,61 +120,6 @@ readListPrec = readListPrecDefault -instance Ord QI where- compare (QI a b c d) (QI a' b' c' d') = res- where- -- (a + b √c)/d ⋛ (a' + b' √c')/d'- -- (a + b √c) d' ⋛ (a' + b' √c') d- -- a d' + b d' √c ⋛ a' d + b' d √c'- -- a d' − a' d ⋛ b' d √c' − b d' √c- --- -- let i = a d' − a' d- -- j = b' d √c'- -- k = b d' √c- i = a * d' - a' * d- sqJ = sq b' * sq d * c'- sqK = sq b * sq d' * c-- -- i ⋛ j − k- --- -- sign (b' d √c') = sign b' because d ≥ 0 and c' ≥ 0- -- sign (b d' √c) = sign b because d' ≥ 0 and c ≥ 0- --- -- if j − k < 0 then (sign b') j² − (sign b) k² < 0- --- -- (sign i) |i| ⋛ sign ((sign b') j² − (sign b) k²) |j − k|- --- -- let snL = sign i- -- snR = sign ((sign b') j² − (sign b) k²)- snL = signum i- snR = signum (signum b' * sqJ - signum b * sqK)-- -- snL |i| ⋛ snR |j − k|- -- snL i² ⋛ snR (j − k)²- -- snL i² ⋛ snR (j² + k² − 2 j k)- -- snL i² − snR (j² + k²) ⋛ snR (−2) j k- -- snL i² − snR (j² + k²) ⋛ snR (−2) b b' d d' √c √c'- --- -- let q = snL i² − snR (j² + k²)- -- r = snR (−2) b b' d d' √c √c'- q = snL * sq i - snR * (sqJ + sqK)- sqR = 4 * sq b * sq b' * sq d * sq d' * c * c'-- -- q ⋛ r- --- -- sign (snR (−2) b b' d d' √c √c') = sign (snR (−2) b b')- --- -- let snL' = sign q- -- snR' = sign (snR (−2) b b')- snL' = signum q- snR' = signum (snR * (-2) * b * b')-- -- snL' |q| ⋛ snR' |r|- -- snL' q² ⋛ snR' r²- res = compare (snL' * sq q) (snR' * sqR)-- sq x = x*x- type QITuple = (Integer, Integer, Integer, Integer) -- | Given @a@, @b@, @c@ and @d@ such that @n = (a + b √c)\/d@, constuct a 'QI'@@ -173,18 +128,23 @@ -- >>> qi 3 4 5 6 -- qi 3 4 5 6 --+-- >>> qi 3 0 42 1+-- qi 3 0 42 1+--+-- >>> qi 1 1 (-1) 1+-- qi 1 1 (-1) 1+-- -- The fractions are reduced: -- -- >>> qi 30 40 5 60 -- qi 3 4 5 6------ If @b = 0@ then @c@ is zeroed and vice versa:+-- >>> qi (-30) (-40) (-5) (-60)+-- qi 3 4 (-5) 6 ----- >>> qi 3 0 42 1--- qi 3 0 0 1+-- If @c = 0@ then @b@ is zeroed and @c@ is put equal to 2: -- -- >>> qi 3 42 0 1--- qi 3 0 0 1+-- qi 3 0 2 1 -- -- The @b √c@ term is simplified: --@@ -194,33 +154,30 @@ -- If @c = 1@ (after simplification) then @b@ is moved to @a@: -- -- >>> qi 1 5 (2*2) 1--- qi 11 0 0 1+-- qi 11 0 2 1 qi :: Integer -- ^ a -> Integer -- ^ b -> Integer -- ^ c -> Integer -- ^ d -> QI-qi a b (nonNegative "qi" -> c) (nonZero "qi" -> d)- | b == 0 = reduceCons a 0 0 d- | c == 0 = reduceCons a 0 0 d- | c == 1 = reduceCons (a + b) 0 0 d- | otherwise = simplifyReduceCons a b c d+qi a _ 0 (nonZero "qi" -> d) = simplifyReduceCons a 0 2 d+qi a b c (nonZero "qi" -> d) = simplifyReduceCons a b c d {-# INLINE qi #-} -- Construct a 'QI' without simplifying @b √c@. Make sure it has already been -- simplified. qiNoSimpl :: Integer -> Integer -> Integer -> Integer -> QI-qiNoSimpl a b (nonNegative "qiNoSimpl" -> c) (nonZero "qiNoSimpl" -> d)- | b == 0 = reduceCons a 0 0 d- | c == 0 = reduceCons a 0 0 d- | c == 1 = reduceCons (a + b) 0 0 d+qiNoSimpl a b (nonZero "qiNoSimpl" -> c) (nonZero "qiNoSimpl" -> d)+ | b == 0 = reduceCons a 0 2 d+ | c == 0 = reduceCons a 0 2 d+ | c == 1 = reduceCons (a + b) 0 2 d | otherwise = reduceCons a b c d {-# INLINE qiNoSimpl #-} -- Simplify @b √c@ before constructing a 'QI'. simplifyReduceCons :: Integer -> Integer -> Integer -> Integer -> QI simplifyReduceCons a b (nonZero "simplifyReduceCons" -> c) d- | c' == 1 = reduceCons (a + b') 0 0 d+ | c' == 1 = reduceCons (a + b') 0 2 d | otherwise = reduceCons a b' c' d where ~(b', c') = separateSquareFactors b c {-# INLINE simplifyReduceCons #-}@@ -228,11 +185,11 @@ -- | Given @b@ and @c@ such that @n = b √c@, return a potentially simplified -- @(b, c)@. separateSquareFactors :: Integer -> Integer -> (Integer, Integer)-separateSquareFactors b (nonNegative "separateSquareFactors" -> c) =- case foldl' go (1,1) (factorise c) of- ~(bMul, c') -> (b*bMul, c')+separateSquareFactors b (nonZero "separateSquareFactors" -> c) =+ case foldl' go (1,1) (factorise (abs c)) of+ ~(bMul, c') -> (b*bMul, c' * signum c) where- go :: (Integer, Integer) -> (Integer, Int) -> (Integer, Integer)+ go :: Integral a => (Integer, Integer) -> (Integer, a) -> (Integer, Integer) go ~(i, j) ~(fac, pow) = i `seq` j `seq` fac `seq` pow `seq` if even pow@@ -241,9 +198,9 @@ -- Reduce the @a@, @b@, @d@ factors before constructing a 'QI'. reduceCons :: Integer -> Integer -> Integer -> Integer -> QI-reduceCons a b c (nonZero "reduceCons" -> d) =- QI (a `quot` q) (b `quot` q) c (d `quot` q)- where q = signum d * (a `gcd` b `gcd` d)+reduceCons a b c (nonZero "reduceCons" -> d) = case someSquareFreeVal c of+ Just (SomeSquareFree (Proxy :: Proxy c)) -> QI (QuadExt (a % d) (b % d) :: QuadExt c Rational)+ Nothing -> error $ "reduceCons: square root term is not square free and equals to " ++ show c {-# INLINE reduceCons #-} -- | Given @a@, @b@ and @c@ such that @n = a + b √c@, constuct a 'QI'@@ -251,11 +208,13 @@ -- -- >>> qi' 0.5 0.7 2 -- qi 5 7 2 10+-- >>> qi' 0.3 0.4 (-3)+-- qi 3 4 (-3) 10 qi' :: Rational -- ^ a -> Rational -- ^ b -> Integer -- ^ c -> QI-qi' a b (nonNegative "qi'" -> c) = n+qi' a b c = n where -- (aN/aD) + (bN/bD) √c = ((aN bD) + (bN aD) √c) / (aD bD) n = qi (aN*bD) (bN*aD) c (aD*bD)@@ -267,22 +226,35 @@ -- -- >>> runQI (qi 3 4 5 6) (\a b c d -> (a,b,c,d)) -- (3,4,5,6)+-- >>> runQI (qi 1 1 (-1) 1) (\a b c d -> (a,b,c,d))+-- (1,1,-1,1) runQI :: QI -> (Integer -> Integer -> Integer -> Integer -> a) -> a-runQI (QI a b c d) f = f a b c d+runQI (QI (QuadExt a b :: QuadExt c Rational)) f = f a' b' c d+ where+ c = squareFreeVal (Proxy :: Proxy c)+ d = lcm (denominator a) (denominator b)+ a' = numerator a * (d `quot` denominator a)+ b' = numerator b * (d `quot` denominator b) {-# INLINE runQI #-} -- | Given @n@ and @f@ such that @n = a + b √c@, run @f a b c@. -- -- >>> runQI' (qi' 0.5 0.7 2) (\a b c -> (a, b, c)) -- (1 % 2,7 % 10,2)+-- >>> runQI' (qi' 0.3 0.4 (-3)) (\a b c -> (a, b, c))+-- (3 % 10,2 % 5,-3) runQI' :: QI -> (Rational -> Rational -> Integer -> a) -> a-runQI' (QI a b c d) f = f (a % d) (b % d) c+runQI' (QI (QuadExt a b :: QuadExt c Rational)) f = f a b c+ where+ c = squareFreeVal (Proxy :: Proxy c) {-# INLINE runQI' #-} -- | Given @n@ such that @n = (a + b √c)\/d@, return @(a, b, c, d)@. -- -- >>> unQI (qi 3 4 5 6) -- (3,4,5,6)+-- >>> unQI (qi 1 1 (-1) 1)+-- (1,1,-1,1) unQI :: QI -> (Integer, Integer, Integer, Integer) unQI n = runQI n (,,,) {-# INLINE unQI #-}@@ -291,6 +263,8 @@ -- -- >>> unQI' (qi' 0.5 0.7 2) -- (1 % 2,7 % 10,2)+-- >>> unQI' (qi' 0.3 0.4 (-3))+-- (3 % 10,2 % 5,-3) unQI' :: QI -> (Rational, Rational, Integer) unQI' n = runQI' n (,,) {-# INLINE unQI' #-}@@ -299,9 +273,13 @@ -- -- >>> view _qi (qi 3 4 5 6) -- (3,4,5,6)+-- >>> view _qi (qi 1 1 (-1) 1)+-- (1,1,-1,1) -- -- >>> over _qi (\(a,b,c,d) -> (a+10, b+10, c+10, d+10)) (qi 3 4 5 6) -- qi 13 14 15 16+-- >>> over _qi (\(a,b,c,d) -> (a+10, b+10, c+10, d+10)) (qi 1 1 (-1) 1)+-- qi 4 0 2 1 _qi :: Lens' QI (Integer, Integer, Integer, Integer) _qi f n = (\ ~(a',b',c',d') -> qi a' b' c' d') <$> f (unQI n) {-# INLINE _qi #-}@@ -310,9 +288,13 @@ -- -- >>> view _qi' (qi' 0.5 0.7 2) -- (1 % 2,7 % 10,2)+-- >>> view _qi' (qi' 0.3 0.4 (-3))+-- (3 % 10,2 % 5,-3) -- -- >>> over _qi' (\(a,b,c) -> (a/5, b/6, c*3)) (qi 3 4 5 6) -- qi 9 10 15 90+-- >>> over _qi' (\(a,b,c) -> (a/5, b/6, c*3)) (qi 1 1 (-1) 1)+-- qi 6 5 (-3) 30 _qi' :: Lens' QI (Rational, Rational, Integer) _qi' f n = (\ ~(a',b',c') -> qi' a' b' c') <$> f (unQI' n) {-# INLINE _qi' #-}@@ -322,9 +304,13 @@ -- -- >>> view _qiABD (qi 3 4 5 6) -- (3,4,6)+-- >>> view _qiABD (qi 1 1 (-1) 1)+-- (1,1,1) -- -- >>> over _qiABD (\(a,b,d) -> (a+10, b+10, d+10)) (qi 3 4 5 6) -- qi 13 14 5 16+-- >>> over _qiABD (\(a,b,d) -> (a+10, b+10, d+10)) (qi 1 1 (-1) 1)+-- qi 1 1 (-1) 1 _qiABD :: Lens' QI (Integer, Integer, Integer) _qiABD f (unQI -> ~(a,b,c,d)) = (\ ~(a',b',d') -> qiNoSimpl a' b' c d') <$> f (a,b,d)@@ -335,9 +321,13 @@ -- -- >>> view _qiA (qi 3 4 5 6) -- 3+-- >>> view _qiA (qi 1 1 (-1) 1)+-- 1 -- -- >>> over _qiA (+ 10) (qi 3 4 5 6) -- qi 13 4 5 6+-- >>> over _qiA (+ 10) (qi 1 1 (-1) 1)+-- qi 11 1 (-1) 1 _qiA :: Lens' QI Integer _qiA = _qiABD . go where go f ~(a,b,d) = (\a' -> (a',b,d)) <$> f a@@ -347,9 +337,13 @@ -- -- >>> view _qiB (qi 3 4 5 6) -- 4+-- >>> view _qiB (qi 1 1 (-1) 1)+-- 1 -- -- >>> over _qiB (+ 10) (qi 3 4 5 6) -- qi 3 14 5 6+-- >>> over _qiB (+ 10) (qi 1 1 (-1) 1)+-- qi 1 11 (-1) 1 _qiB :: Lens' QI Integer _qiB = _qiABD . go where go f ~(a,b,d) = (\b' -> (a,b',d)) <$> f b@@ -359,9 +353,13 @@ -- -- >>> view _qiC (qi 3 4 5 6) -- 5+-- >>> view _qiC (qi 1 1 (-1) 1)+-- -1 -- -- >>> over _qiC (+ 10) (qi 3 4 5 6) -- qi 3 4 15 6+-- >>> over _qiC (+ 10) (qi 1 1 (-1) 1)+-- qi 4 0 2 1 _qiC :: Lens' QI Integer _qiC = _qi . go where go f ~(a,b,c,d) = (\c' -> (a,b,c',d)) <$> f c@@ -371,42 +369,32 @@ -- -- >>> view _qiD (qi 3 4 5 6) -- 6+-- >>> view _qiD (qi 1 1 (-1) 1)+-- 1 -- -- >>> over _qiD (+ 10) (qi 3 4 5 6) -- qi 3 4 5 16+-- >>> over _qiD (+ 10) (qi 1 1 (-1) 1)+-- qi 1 1 (-1) 11 _qiD :: Lens' QI Integer _qiD = _qiABD . go where go f ~(a,b,d) = (\d' -> (a,b,d')) <$> f d --- | The constant zero.------ >>> qiZero--- qi 0 0 0 1-qiZero :: QI-qiZero = qi 0 0 0 1-{-# INLINE qiZero #-}---- | The constant one.------ >>> qiOne--- qi 1 0 0 1-qiOne :: QI-qiOne = qi 1 0 0 1-{-# INLINE qiOne #-}- -- | Check if the value is zero. ----- >>> map qiIsZero [qiZero, qiOne, qiSubR (qi 7 0 0 2) 3.5]--- [True,False,True]+-- >>> map qiIsZero [qi 0 0 0 1, qi 1 0 0 1, qiSubR (qi 7 0 0 2) 3.5, qi 0 9 0 1, qi 0 0 (-1) 5]+-- [True,False,True,True,True] qiIsZero :: QI -> Bool--- If b = 0 then c = 0 and vice versa, guaranteed by the constructor.-qiIsZero (unQI -> ~(a,b,_,_)) = a == 0 && b == 0+qiIsZero (QI q) = q == 0 {-# INLINE qiIsZero #-} -- | Convert a 'QI' number into a 'Floating' one.+-- Returns NaN for imaginary values. -- -- >>> qiToFloat (qi 3 4 5 6) == ((3 + 4 * sqrt 5)/6 :: Double) -- True+-- >>> qiToFloat (qi 2 3 (-5) 7)+-- NaN qiToFloat :: Floating a => QI -> a qiToFloat (unQI -> ~(a,b,c,d)) = (fromInteger a + fromInteger b * sqrt (fromInteger c)) / fromInteger d@@ -416,219 +404,213 @@ -- -- >>> qi 3 4 5 6 `qiAddI` 1 -- qi 9 4 5 6+-- >>> qi 1 1 (-1) 1 `qiAddI` 1+-- qi 2 1 (-1) 1 qiAddI :: QI -> Integer -> QI-qiAddI n x = over _qiABD go n- where go ~(a,b,d) = a `seq` b `seq` d `seq` x `seq` (a + d*x, b, d)+qiAddI (QI q) x = QI (q + fromInteger x) {-# INLINE qiAddI #-} -- | Add a 'Rational' to a 'QI'. -- -- >>> qi 3 4 5 6 `qiAddR` 1.2 -- qi 51 20 5 30+-- >>> qi 1 1 (-1) 1 `qiAddR` 1.2+-- qi 11 5 (-1) 5 qiAddR :: QI -> Rational -> QI-qiAddR n x = over _qiABD go n- where- -- n = (a + b √c)/d + xN/xD- -- n = ((a + b √c) xD)/(d xD) + (d xN)/(d xD)- -- n = ((a xD + d xN) + b xD √c)/(d xD)- go ~(a,b,d) =- a `seq` b `seq` d `seq` xN `seq` xD `seq` (a*xD + d*xN, b*xD, d*xD)- (xN, xD) = (numerator x, denominator x)+qiAddR (QI q) x = QI (q + fromRational x) {-# INLINE qiAddR #-} -- | Subtract an 'Integer' from a 'QI'. -- -- >>> qi 3 4 5 6 `qiSubI` 1 -- qi (-3) 4 5 6+-- >>> qi 1 1 (-1) 1 `qiSubI` 1+-- qi 0 1 (-1) 1 qiSubI :: QI -> Integer -> QI-qiSubI n x = qiAddI n (negate x)+qiSubI (QI q) x = QI (q - fromInteger x) {-# INLINE qiSubI #-} -- | Subtract a 'Rational' from a 'QI'. -- -- >>> qi 3 4 5 6 `qiSubR` 1.2 -- qi (-21) 20 5 30+-- >>> qi 1 1 (-1) 1 `qiSubR` 1.2+-- qi (-1) 5 (-1) 5 qiSubR :: QI -> Rational -> QI-qiSubR n x = qiAddR n (negate x)+qiSubR (QI q) x = QI (q - fromRational x) {-# INLINE qiSubR #-} -- | Multiply a 'QI' by an 'Integer'. -- -- >>> qi 3 4 5 6 `qiMulI` 2 -- qi 3 4 5 3+-- >>> qi 1 1 (-1) 1 `qiMulI` 2+-- qi 2 2 (-1) 1 qiMulI :: QI -> Integer -> QI-qiMulI n x = over _qiABD go n- where go ~(a,b,d) = a `seq` b `seq` d `seq` x `seq` (a*x, b*x, d)+qiMulI (QI q) x = QI (q * fromInteger x) {-# INLINE qiMulI #-} -- | Multiply a 'QI' by a 'Rational'. -- -- >>> qi 3 4 5 6 `qiMulR` 0.5 -- qi 3 4 5 12+-- >>> qi 1 1 (-1) 1 `qiMulR` 0.5+-- qi 1 1 (-1) 2 qiMulR :: QI -> Rational -> QI-qiMulR n x = over _qiABD go n- where- -- n = (a + b √c)/d xN/xD- -- n = (a xN + b xN √c)/(d xD)- go ~(a,b,d) = a `seq` b `seq` d `seq` xN `seq` xD `seq` (a*xN, b*xN, d*xD)- (xN, xD) = (numerator x, denominator x)+qiMulR (QI q) x = QI (q * fromRational x) {-# INLINE qiMulR #-} --- | Divice a 'QI' by an 'Integer'.+-- | Divide a 'QI' by an 'Integer'. -- -- >>> qi 3 4 5 6 `qiDivI` 2 -- qi 3 4 5 12+-- >>> qi 1 1 (-1) 1 `qiDivI` 2+-- qi 1 1 (-1) 2 qiDivI :: QI -> Integer -> QI-qiDivI n (nonZero "qiDivI" -> x) = over _qiABD go n- where go ~(a,b,d) = a `seq` b `seq` d `seq` x `seq` (a, b, d*x)+qiDivI (QI q) (nonZero "qiDivI" -> x) = QI (q * fromRational (1 % x)) {-# INLINE qiDivI #-} --- | Divice a 'QI' by a 'Rational'.+-- | Divide a 'QI' by a 'Rational'. -- -- >>> qi 3 4 5 6 `qiDivR` 0.5 -- qi 3 4 5 3+-- >>> qi 1 1 (-1) 1 `qiDivR` 0.5+-- qi 2 2 (-1) 1 qiDivR :: QI -> Rational -> QI-qiDivR n (nonZero "qiDivR" -> x) = qiMulR n (recip x)+qiDivR (QI q) (nonZero "qiDivR" -> x) = QI (q * fromRational (recip x)) {-# INLINE qiDivR #-} -- | Negate a 'QI'. -- -- >>> qiNegate (qi 3 4 5 6) -- qi (-3) (-4) 5 6+-- >>> qiNegate (qi 1 1 (-1) 1)+-- qi (-1) (-1) (-1) 1 qiNegate :: QI -> QI-qiNegate n = over _qiABD go n- where go ~(a,b,d) = a `seq` b `seq` d `seq` (negate a, negate b, d)+qiNegate (QI q) = QI (negate q) {-# INLINE qiNegate #-} -- | Compute the reciprocal of a 'QI'. -- -- >>> qiRecip (qi 5 0 0 2)--- Just (qi 2 0 0 5)+-- Just (qi 2 0 2 5) -- -- >>> qiRecip (qi 0 1 5 2) -- Just (qi 0 2 5 5) ----- >>> qiRecip qiZero+-- >>> qiRecip (qi 1 1 (-1) 1)+-- Just (qi 1 (-1) (-1) 2)+--+-- >>> qiRecip (qi 0 0 0 1) -- Nothing qiRecip :: QI -> Maybe QI-qiRecip n@(unQI -> ~(a,b,c,d))- -- 1/((a + b √c)/d) =- -- d/(a + b √c) =- -- d (a − b √c) / ((a + b √c) (a − b √c)) =- -- d (a − b √c) / (a² − b² c) =- -- (a d − b d √c) / (a² − b² c)- | qiIsZero n = Nothing- | denom == 0 = error ("qiRecip: Failed for " ++ show n)- | otherwise = Just (set _qiABD (a * d, negate (b * d), denom) n)- where denom = (a*a - b*b * c)+qiRecip (QI 0) = Nothing+qiRecip (QI q) = Just (QI (recip q)) -- | Add two 'QI's if the square root terms are the same or zeros. ----- >>> qi 3 4 5 6 `qiAdd` qiOne+-- >>> qi 3 4 5 6 `qiAdd` qi 1 0 5 1 -- Just (qi 9 4 5 6) -- -- >>> qi 3 4 5 6 `qiAdd` qi 3 4 5 6 -- Just (qi 3 4 5 3) --+-- >>> qi 3 4 (-5) 6 `qiAdd` qi 3 4 (-5) 6+-- Just (qi 3 4 (-5) 3)+-- -- >>> qi 0 1 5 1 `qiAdd` qi 0 1 6 1 -- Nothing qiAdd :: QI -> QI -> Maybe QI-qiAdd n@(unQI -> ~(a,b,c,d)) n'@(unQI -> ~(a',b',c',d'))- -- n = (a + b √c)/d + (a' + b' √c')/d'- -- n = ((a + b √c) d' + (a' + b' √c') d)/(d d')- -- if c = c' then n = ((a d' + a' d) + (b d' + b' d) √c)/(d d')- | c == 0 = Just (set _qiABD (a*d' + a'*d, b'*d, d*d') n')- | c' == 0 = Just (set _qiABD (a*d' + a'*d, b*d' , d*d') n)- | c == c' = Just (set _qiABD (a*d' + a'*d, b*d' + b'*d, d*d') n)- | otherwise = Nothing+qiAdd (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =+ case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of+ Nothing -> Nothing+ Just Refl -> Just (QI (q1 + q2)) -- | Subtract two 'QI's if the square root terms are the same or zeros. ----- >>> qi 3 4 5 6 `qiSub` qiOne+-- >>> qi 3 4 5 6 `qiSub` qi 1 0 5 1 -- Just (qi (-3) 4 5 6) -- -- >>> qi 3 4 5 6 `qiSub` qi 3 4 5 6--- Just (qi 0 0 0 1)+-- Just (qi 0 0 5 1) --+-- >>> qi 3 4 (-5) 6 `qiSub` qi 3 4 (-5) 6+-- Just (qi 0 0 (-5) 1)+-- -- >>> qi 0 1 5 1 `qiSub` qi 0 1 6 1 -- Nothing qiSub :: QI -> QI -> Maybe QI-qiSub n n' = qiAdd n (qiNegate n')+qiSub (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =+ case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of+ Nothing -> Nothing+ Just Refl -> Just (QI (q1 - q2)) -- | Multiply two 'QI's if the square root terms are the same or zeros. ----- >>> qi 3 4 5 6 `qiMul` qiZero--- Just (qi 0 0 0 1)+-- >>> qi 3 4 5 6 `qiMul` qi 0 0 5 1+-- Just (qi 0 0 5 1) ----- >>> qi 3 4 5 6 `qiMul` qiOne+-- >>> qi 3 4 5 6 `qiMul` qi 1 0 5 1 -- Just (qi 3 4 5 6) -- -- >>> qi 3 4 5 6 `qiMul` qi 3 4 5 6 -- Just (qi 89 24 5 36) --+-- >>> qi 3 4 (-5) 6 `qiMul` qi 3 4 (-5) 6+-- Just (qi (-71) 24 (-5) 36)+-- -- >>> qi 0 1 5 1 `qiMul` qi 0 1 6 1 -- Nothing qiMul :: QI -> QI -> Maybe QI-qiMul n@(unQI -> ~(a,b,c,d)) n'@(unQI -> ~(a',b',c',d'))- -- n = (a + b √c)/d (a' + b' √c')/d'- -- n = (a a' + a b' √c' + a' b √c + b b' √c √c')/(d d')- -- if c = 0 then n = (a a' + a b' √c')/(d d')- -- if c' = 0 then n = (a a' + a' b √c)/(d d')- -- if c = c' then n = ((a a' + b b' c) + (a b' + a' b) √c)/(d d')- | c == 0 = Just (set _qiABD (a*a' , a*b' , d*d') n')- | c' == 0 = Just (set _qiABD (a*a' , a'*b, d*d') n)- | c == c' = Just (set _qiABD (a*a' + b*b'*c, a*b' + a'*b, d*d') n)- | otherwise = Nothing+qiMul (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =+ case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of+ Nothing -> Nothing+ Just Refl -> Just (QI (q1 * q2)) -- | Divide two 'QI's if the square root terms are the same or zeros. ----- >>> qi 3 4 5 6 `qiDiv` qiZero+-- >>> qi 3 4 5 6 `qiDiv` qi 0 0 5 1 -- Nothing ----- >>> qi 3 4 5 6 `qiDiv` qiOne+-- >>> qi 3 4 5 6 `qiDiv` qi 1 0 5 1 -- Just (qi 3 4 5 6) -- -- >>> qi 3 4 5 6 `qiDiv` qi 3 4 5 6--- Just (qi 1 0 0 1)+-- Just (qi 1 0 5 1) -- -- >>> qi 3 4 5 6 `qiDiv` qi 0 1 5 1 -- Just (qi 20 3 5 30) --+-- >>> qi 3 4 (-5) 6 `qiDiv` qi 0 1 (-5) 1+-- Just (qi 20 (-3) (-5) 30)+-- -- >>> qi 0 1 5 1 `qiDiv` qi 0 1 6 1 -- Nothing qiDiv :: QI -> QI -> Maybe QI-qiDiv n n' = qiMul n =<< qiRecip n'+qiDiv _ (QI 0) = Nothing+qiDiv (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =+ case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of+ Nothing -> Nothing+ Just Refl -> Just (QI (q1 / q2)) -- | Exponentiate a 'QI' to an 'Integer' power. -- -- >>> qi 3 4 5 6 `qiPow` 0--- qi 1 0 0 1+-- qi 1 0 5 1 -- -- >>> qi 3 4 5 6 `qiPow` 1 -- qi 3 4 5 6 -- -- >>> qi 3 4 5 6 `qiPow` 2 -- qi 89 24 5 36+--+-- >>> qi 1 1 (-1) 1 `qiPow` 4+-- qi (-4) 0 (-1) 1 qiPow :: QI -> Integer -> QI-qiPow num (nonNegative "qiPow" -> pow) = go num pow- where- go _ 0 = qiOne- go n 1 = n- go n p- | even p = go (sudoQIMul n n) (p `div` 2)- | otherwise = go' (sudoQIMul n n) ((p-1) `div` 2) n-- -- Like go but multiplied with n'.- go' _ 0 n' = n'- go' n 1 n' = sudoQIMul n n'- go' n p n'- | even p = go' (sudoQIMul n n) (p `div` 2) n'- | otherwise = go' (sudoQIMul n n) ((p-1) `div` 2) (sudoQIMul n n')-- -- Multiplying a QI with its own power will always succeed.- sudoQIMul n n' = case qiMul n n' of ~(Just m) -> m+qiPow (QI q) pow = QI (getProduct (stimes pow (Product q))) -- | Compute the floor of a 'QI'.+-- Throws an error for imaginary values. -- -- >>> qiFloor (qi 10 0 0 2) -- 5@@ -638,6 +620,10 @@ -- -- >>> qiFloor (qi 10 2 5 2) -- 7+--+-- >>> qiFloor (qi 1 1 (-1) 1)+-- *** Exception: ...+-- ... qiFloor :: QI -> Integer qiFloor (unQI -> ~(a,b,c,d)) = -- n = (a + b √c)/d@@ -654,7 +640,7 @@ -- @[2; 2] = 2 + 1\/2 = 5\/2@. -- -- >>> continuedFractionToQI (2,CycList [2] [])--- qi 5 0 0 2+-- qi 5 0 2 2 -- -- The golden ratio is @[1; 1, 1, …]@. --@@ -666,7 +652,7 @@ continuedFractionToQI :: (Integer, CycList Integer) -> QI continuedFractionToQI (i0_, is_) = qiAddI (go is_) i0_ where- go (CycList as bs) = goNonCyc as (if null bs then qiZero else goCyc bs)+ go (CycList as bs) = goNonCyc as (if null bs then qi 0 0 0 1 else goCyc bs) goNonCyc ((pos -> i):is) final = sudoQIRecip (qiAddI (goNonCyc is final) i) goNonCyc [] final = final@@ -705,6 +691,7 @@ pos = positive "continuedFractionToQI" -- | Convert a 'QI' into a (possibly periodic) simple continued fraction.+-- Throws an error for imaginary values. -- -- @5\/2 = 2 + 1\/2 = [2; 2]@. --@@ -718,6 +705,10 @@ -- -- >>> qiToContinuedFraction (qi 987601513930253257378987883 1 14116473325908285531353005 81983584717737887813195873886) -- (0,CycList [83,78,65,75,69] [32,66,65,68,71,69,82])+--+-- >>> qiToContinuedFraction (qi 1 1 (-1) 1)+-- *** Exception: ...+-- ... qiToContinuedFraction :: QI -> (Integer, CycList Integer) qiToContinuedFraction num@@ -772,10 +763,6 @@ low = integerSquareRoot n high | low*low == n = low | otherwise = low + 1--nonNegative :: (Num a, Ord a, Show a) => String -> a -> a-nonNegative name = validate name "non-negative" (>= 0)-{-# INLINE nonNegative #-} positive :: (Num a, Ord a, Show a) => String -> a -> a positive name = validate name "positive" (> 0)
src/Numeric/QuadraticIrrational/CyclicList.hs view
@@ -13,6 +13,7 @@ ( CycList(..) ) where +import Data.Foldable (toList) import Data.Monoid ((<>)) -- | A container for a possibly cyclic list.@@ -29,6 +30,10 @@ deriving (Eq, Ord, Read, Show, Functor) instance Foldable CycList where- foldMap f (CycList as bs) = foldMap f as <> if null bs then mempty else cycleAppend (foldMap f bs)+ toList (CycList as []) = as+ toList (CycList as bs) = as <> cycle bs++ foldMap f (CycList as []) = foldMap f as+ foldMap f (CycList as bs) = foldMap f as <> cycleAppend (foldMap f bs) where cycleAppend x = x <> cycleAppend x
+ src/Numeric/QuadraticIrrational/Ring.hs view
@@ -0,0 +1,48 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}++module Numeric.QuadraticIrrational.Ring+ ( QuadExt(..)+ ) where++import Data.Proxy+import Numeric.QuadraticIrrational.SquareFree++data QuadExt (c :: SquareFree) t = QuadExt t t+ deriving (Eq)++instance (KnownSquareFree c, Show t) => Show (QuadExt c t) where+ showsPrec p (QuadExt a b) = showParen (p > 10)+ $ showString "QuadExt " . showsPrec 11 a . showChar '+' . showsPrec 11 b . showChar '√' . showsPrec 11 c+ where+ c = squareFreeVal (Proxy :: Proxy c)++instance (KnownSquareFree c, Num t) => Num (QuadExt c t) where+ QuadExt a1 b1 + QuadExt a2 b2 = QuadExt (a1 + a2) (b1 + b2)+ QuadExt a1 b1 - QuadExt a2 b2 = QuadExt (a1 - a2) (b1 - b2)+ negate (QuadExt a b) = QuadExt (negate a) (negate b)+ abs q = q+ signum _ = 1+ fromInteger n = QuadExt (fromInteger n) 0++ QuadExt a1 b1 * QuadExt a2 b2 = QuadExt a b+ where+ a = a1 * a2 + fromInteger c * b1 * b2+ b = a1 * b2 + a2 * b1+ c = squareFreeVal (Proxy :: Proxy c)++instance (KnownSquareFree c, Fractional t) => Fractional (QuadExt c t) where+ fromRational r = QuadExt (fromRational r) 0++ recip (QuadExt a b) = QuadExt (a / d) (negate b / d)+ where+ d = a * a - fromInteger c * b * b+ c = squareFreeVal (Proxy :: Proxy c)++ QuadExt a1 b1 / QuadExt a2 b2 = QuadExt (a / d) (b / d)+ where+ QuadExt a b = QuadExt a1 b1 * QuadExt a2 (negate b2) :: QuadExt c t+ d = a2 * a2 - fromInteger c * b2 * b2+ c = squareFreeVal (Proxy :: Proxy c)
+ src/Numeric/QuadraticIrrational/SquareFree.hs view
@@ -0,0 +1,48 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}++module Numeric.QuadraticIrrational.SquareFree+ ( SquareFree()+ , KnownSquareFree(..)+ , SomeSquareFree(..)+ , someSquareFreeVal+ , sameSquareFree+ ) where++import Data.Proxy+import Data.Type.Equality+import GHC.TypeLits+import Math.NumberTheory.ArithmeticFunctions+import Unsafe.Coerce++data SquareFree = Negative Nat | Positive Nat++class KnownSquareFree (c :: SquareFree) where+ squareFreeVal :: Proxy c -> Integer++instance KnownNat n => KnownSquareFree ('Negative n) where+ squareFreeVal (Proxy :: Proxy ('Negative n)) = negate $ natVal (Proxy :: Proxy n)++instance KnownNat n => KnownSquareFree ('Positive n) where+ squareFreeVal (Proxy :: Proxy ('Positive n)) = natVal (Proxy :: Proxy n)++data SomeSquareFree where+ SomeSquareFree :: KnownSquareFree c => Proxy c -> SomeSquareFree++someSquareFreeVal :: Integer -> Maybe SomeSquareFree+someSquareFreeVal 0 = Nothing+someSquareFreeVal n = case moebius n of+ MoebiusZ -> Nothing+ _ -> case someNatVal n of+ Just (SomeNat (Proxy :: Proxy n)) -> Just (SomeSquareFree (Proxy :: Proxy ('Positive n)))+ Nothing -> case someNatVal (negate n) of+ Just (SomeNat (Proxy :: Proxy n)) -> Just (SomeSquareFree (Proxy :: Proxy ('Negative n)))+ Nothing -> Nothing++sameSquareFree :: (KnownSquareFree a, KnownSquareFree b) => Proxy a -> Proxy b -> Maybe (a :~: b)+sameSquareFree x y+ | squareFreeVal x == squareFreeVal y = Just (unsafeCoerce Refl)+ | otherwise = Nothing
tests/QuadraticIrrational.hs view
@@ -5,9 +5,7 @@ import Data.Number.CReal (CReal) import Test.Tasty (TestTree, testGroup)-import Test.Tasty.QuickCheck (Arbitrary (arbitrary, shrink),- NonNegative (NonNegative), NonZero (NonZero), Property, Testable, (===),- (==>), conjoin, counterexample, testProperty)+import Test.Tasty.QuickCheck import Numeric.QuadraticIrrational import Numeric.QuadraticIrrational.Internal.Lens@@ -87,16 +85,9 @@ approxEq' (qiToFloat (qi a b c d)) (approxQI a b c d) , testProperty "compare equals" $ \a ->- conjoin [ a === a, compare a a === EQ ]+ conjoin [ a === a ] `const` (a :: QI) - , testProperty "compare" $ \a b ->- let a' = qiToFloat a :: RefFloat- b' = qiToFloat b :: RefFloat- in conjoin [ (a == b) === (a' == b')- , compare a b === compare a' b'- ]- , testProperty "qiAddI" $ \n x -> approxEq' (qiToFloat (qiAddI n x)) (qiToFloat n + fromInteger x) @@ -163,7 +154,7 @@ let cf@(_, CycList _ xs) = qiToContinuedFraction n -- Limit the length of the periodic part for speed. in (length xs <= 100) ==>- n === continuedFractionToQI cf+ (qiToFloat n :: Double) === qiToFloat (continuedFractionToQI cf) , testProperty "continuedFractionApproximate" $ \n -> let cf = qiToContinuedFraction n@@ -173,16 +164,16 @@ ] withCompatiblePair :: Testable p- => (QI -> QI -> p) -> QI -> QI -> Bool -> Bool -> Property-withCompatiblePair f n0_ n1_ c0Zero c1Zero =+ => (QI -> QI -> p) -> QI -> QI -> Property+withCompatiblePair f n0_ n1_ = counterexample ("n0 = " ++ show n0) . counterexample ("n1 = " ++ show n1) $ f n0 n1 where n0 = runQI n0_ $ \a b c d ->- qi a b (if c0Zero then 0 else c) d+ qi a b c d n1 = runQI n0_ $ \_ _ c _ -> runQI n1_ $ \a b _ d ->- qi a b (if c1Zero then 0 else c) d+ qi a b c d approxQI :: Integer -> Integer -> Integer -> Integer -> RefFloat approxQI a b c d =
tests/doctests.hs view
@@ -2,9 +2,8 @@ module Main (main) where -import Control.Applicative (empty) import Control.Monad (guard)-import Control.Monad.List (ListT (ListT), liftIO, runListT)+import Data.Foldable (toList) import Data.List (isSuffixOf) import System.Directory (doesDirectoryExist, doesFileExist, getDirectoryContents)@@ -17,17 +16,22 @@ doctest ("-isrc" : sources) findSources :: FilePath -> IO [FilePath]-findSources dir = runListT (goDir dir)+findSources = goDir where- goItem :: FilePath -> FilePath -> ListT IO FilePath- goItem _ ('.':_) = empty+ goItem :: FilePath -> FilePath -> IO [FilePath]+ goItem _ ('.':_) = pure [] goItem parent name = do let path = parent </> name- isDir <- liftIO (doesDirectoryExist path)- isFile <- liftIO (doesFileExist path)- if | isDir -> goDir path- | isFile -> goFile path- | otherwise -> empty+ isDir <- doesDirectoryExist path+ isFile <- doesFileExist path+ if | isDir -> goDir path+ | isFile -> pure $ toList $ goFile path+ | otherwise -> pure [] - goDir path = goItem path =<< ListT (getDirectoryContents path)+ goDir :: FilePath -> IO [FilePath]+ goDir path = do+ items <- getDirectoryContents path+ concat <$> traverse (goItem path) items++ goFile :: FilePath -> Maybe FilePath goFile path = path <$ guard (".hs" `isSuffixOf` path)